Tutorial 17 – Iterative Solvers

Introduction

In this tutorial we will present iterative solver in microfuidic application. For more detail about microfluid please refer to tutorial 16.

Code

% Clear MATLAB Command Window and Workspace
clc;
clear;



% Read and display mesh
[p,e,t] = importMeshGmsh('microfluidicH3D_coarse.msh');
%[p,e,t] = importMeshGmsh('microfluidicH3D.msh');
displayMesh(p,t);

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


% Define fluid
nu = 0.1;


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


% Define inlet velocity profile
vel = [1 0 0];


% Options of the iterative solver
opt.alpha_p = 0.02;
opt.pres_iter = 8;
opt.bicg_iter = 48;
opt.bicg_eps = 1e-10;

% Iterate nonlinear terms
for iter = 1:maxiter
    % Compute shear rate and viscosity for non-Newtonian fluids
    % gamma = shearRate(p,t,u);
    % nu = nuinf + (nu0-nuinf)*(1+(lambda*gamma).^2).^((n-1)/2);

    % Assemble matrix and right hand side
    [NS, F] = assembleNavierStokesMatrix(p,e,t,nu,u(indices.indu),u(indices.indv),u(indices.indw),'nosupg');
    % [NS, F] = assembleNavierStokesMatrix(p,e,t,nu,u(indices.indu),u(indices.indv),u(indices.indw),'supgDoublyAsymptotic');

    % Apply boundary conditions
    [NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,[77 78],'inlet', [1 0 0]);
    [NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,[80 81],'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
    tic
    %[u, convergence] = iterativeSolver(NS, F, u, indices, convergence, 'direct');
    %[u, convergence] = iterativeSolver(NS, F, u, indices, convergence, 'SIMPLE');
    [u, convergence] = iterativeSolver(NS, F, u, indices, convergence, 'SIMPLE',opt);
    toc
end


% Plot solution
figure(3);
displaySolution(p,t,u(indices.indu),'x-velocity');

% Solve scalar transport problem
k = 0.001;
%k = 0.1;


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


% Apply boundary conditions
[D,F] = imposeScalarBoundaryCondition(p,e,D,F,77,'value',0);
[D,F] = imposeScalarBoundaryCondition(p,e,D,F,78,'value',1);


% Solve system of equations
phi = D\F;
% Display solution

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

figure(5);
pressure = generatePressureData(u,p,t);
displaySolution(p,t,pressure,'Pressure');