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; % 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)