This tutorial is intended for the full version of the toolbox.
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');