Introduction and Problem Definition
This tutorial is intended to show, how a scalar transport problem can be solved in QuickerSim CFD Toolbox for MATLAB. We will be dealing here with a laminar flow in a microfluidic application at Reynolds number of 50. Domain geometry with ids of Physical Lines, where we will be specyfying boundary conditions is given below.
The flow will be simulated with incompressible Navier-Stokes equations. Inlet velocity at each of the inlets (id = 15 and 16) will be equal to 50 mm/s. Pressure value at outlet 17 will be equal to zero and on outlet 18, average pressure value will be 10 Pa (please remember that definition of pressure in the toolbox is physical pressure divided by density – which means that for water, we will need to apply average pressure of 0.01 m^2/s^2).
Futhermore, we define a scalar transport problem with zero concentration at inlet 15, concentration equal to unity at inlet 16 and diffusivity of . These assumptions result in the following mathematical model:
with at inlet 15, at inlet 16 and homogeonuous Neumann condition elsewhere.
% Clear MATLAB Command Window and Workspace clc; clear; % Read, scale (from mm to meters and SI units) and display mesh [p,e,t] = importMeshGmsh('microfluidicH.msh'); p = p/1000; % p array stores mesh coordinates displayMesh2D(p,t);
% Convert mesh to second order [p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t); % Define fluid kinematic viscosity nu = 1e-6; % Init solution and convergence criteria [u, convergence] = initSolution(p,t,[0.05 0],0); maxres = 1e-8; 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); [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,,'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
Now, we present the code which solves passive scalar transport problem. Please note that we exploit the flag ‘supgDoublyAsymptotic’ to turn on stabilization of convective terms. Otherwise, if the mesh is not fine enough strong concentration gradients could be smoothed out due to excessive numerical diffusion or unstable oscillations in the solution field could be present. Refer to Tutorial 4 on SUPG stabilization for details.
% 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; % Display solution % We display concentration field to track species concentration and the % pressure field to evaluate visually pressure losses. 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)