Microfluidic simulation in QuickerSim CFD Toolbox for MATLAB

Today, we want to provide you an example source code which shows you how to do your own microfluidic simulation in the free version of our QuickerSim CFD Toolbox. You will see, how to simulate a full Navier-Stokes fluid flow problem in laminar regime and add a scalar transport equation to compute concentration field of a species. Here, the example comes only with pictures and pure source code. Full description is included in pdf tutorial files for this case, after you download our toolbox from the homepage: http://www.quickersim.com/cfd-toolbox-for-matlab/index

% Solution of the flow problem
clc;
clear;

&nbsp;

% Read mesh
[p,e,t] = importMeshGmsh('microfluidicH.msh');
p = p/1000; % p array stores mesh coordinates - rescale to [mm]
displayMesh2D(p,t);

% Convert mesh to second order
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);

% Define fluid kinematic viscosity of water
nu = 1e-6; % water viscosity [m^2/s]

% Init solution and convergence criteria
[u, convergence] = initSolution(p,t,[0.05 0],0);
maxres = 1e-9;
maxiter = 25;

% Define inlet velocity profile
vel = [0.05 0]; % x-velocity of 50 mm/s = 0.05 m/s

% Iterate nonlinear terms
for iter = 1:maxiter
% Assemble matrix and right hand side
[NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu,u(indices.indu),u(indices.indv),'nosupg');

% Apply boundary conditions
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,15,'inlet', vel);
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,16,'inlet', vel);
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,18,'outflow',0.01); % non-zero pressure at the top outlet
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[19],'wall');

% Compute and plot residuals
[stop, convergence] = computeResiduals(NS,F,u,size(p),convergence,maxres);
plotResiduals(convergence,2);

% Break if solution converged
if(stop)
break;
end

% Solve equations
u = NS\F;
end

% Plot velocity field
figure(3);
vmag = sqrt(u(indices.indu).^2+u(indices.indv).^2);
displaySolution2D(p,t,vmag,'Velocity magnitude');

% Solution of the Scalar Transport Problem

% Define scalar diffusivity
k = 1e-7;

% Assemble problem matrix
[D,F] = assembleDiffusionMatrix2D(p,t,k);
D = D + assembleScalarConvectionMatrix2D(p,t,k,u(indices.indu),u(indices.indv),'supgDoublyAsymptotic');

% Apply boundary conditions
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,15,'value',0);
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,16,'value',1);

% Solve system of equations
phi = D\F;

% Simple Postprocessing

figure(4);
displaySolution2D(p,t,phi,'Concentration');

figure(5);
pressure = generatePressureData(u,p,t);
rho = 1000;
displaySolution2D(p,t,rho*pressure,'Pressure [Pa]');

% We can also monitor average pressure at each of the inlets with the
% following code
[~,p_in_15] = boundaryIntegral2D(p,e,rho*pressure,15);
[~,p_in_16] = boundaryIntegral2D(p,e,rho*pressure,16);

disp(p_in_15)
disp(p_in_16)

Tags

Search