Tutorial 16 – Laminar Flow in Microfluidic Device

Posted on July 17, 2017 - Tutorials

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 $k = 1e-7 m^2/s^2$. These assumptions result in the following mathematical model:

    $$\vec u \nabla \phi = \nabla \cdot (k \nabla \phi)$$

with $\phi = 0$ at inlet 15, $\phi = 1$ at inlet 16 and homogeonuous Neumann condition $\vec n \cdot \nabla \phi = 0$ elsewhere.

Code

% 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,[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

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)

Tags