This tutorial is intended for the full version of the toolbox.
Code
% Clear MATLAB Command Window and Workspace clear; clc; % Import and display mesh [p,e,t] = importMeshGmsh('carotidBifurcation.msh'); displayMesh(p,t);
% Convert mesh to second order [p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t); % Define fluid nu = 3.8; % mm^2/s % Init solution or solve a stationary problem [u, convergence] = initSolution(p,t,[0 0 0],0); % Define timestepping settings dt = 0.125; nstep = 20; % Assemble mass matrix M = assembleMassMatrix(p,t,1/dt); movieObj = initMovie('carotidFlow.avi',12); for step = 1:nstep time = step*dt; % Compute shear rate and viscosity for non-Newtonian fluids if needed gamma = shearRate(p,t,u); % nu = nuinf + (nu0-nuinf)*(1+(lambda*gamma).^2).^((n-1)/2); % Assemble Navier-Stokes matrix [NS, F] = assembleNavierStokesMatrix(p,e,t,nu,u(indices.indu),u(indices.indv),u(indices.indw),'supgDoublyAsymptotic'); % Add unsteady terms [NS,F] = addUnsteadyTerms(NS,F,M,u,indices); % Impose boundary conditions vz = 400*max(sin(time),0); [NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,256,'inlet', [0 0 vz]); [NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,259,'wall'); u = NS\F; % Preview flow field and shoot movie figure(3) displaySolution(p,t,gamma,'Shear rate [1/s]',[0 1600]); drawnow; writeVideo(movieObj,getframe(gcf)); % Report time step number disp(step); end
close(movieObj);