Tutorial 18 – Blood Flow Simulation

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);