Fluid dynamics: Simulation Setup

We provide here an example of data and main files which are used to setup a simulation.

Datafile options

First set some constants to be used later on:

H = 0.41;
U = 2.25;

Define the 3 components of the volumetric forces (i.e. the right-hand side term $\mathbf{f}$ )

% Source term
data.force{1} = @(x, y, z, t, param)(0.*x.*y);
data.force{2} = @(x, y, z, t, param)(0.*x.*y);
data.force{3} = @(x, y, z, t, param)(0.*x.*y);

Set the flags for the boundary conditions: on boundaries [2, 10, 6] we impose Dirichlet (essential) conditions on all the 3 components of the velocity. At the outflow (boundary 3), we impose a Neumann (natural) condition

% flags
data.flag_dirichlet{1} = [2 10 6];
data.flag_neumann{1}   = [3];

data.flag_dirichlet{2} = [2 10 6];
data.flag_neumann{2}   = [3];

data.flag_dirichlet{3} = [2 10 6];
data.flag_neumann{3}   = [3];

Specify the functions $\mathbf{g}$ and $\mathbf{h}$ for Dirichlet and Neumann boundary conditions

% Dirichlet
data.bcDir{1} = @(x, y, z, t, param)( sin( pi*t/8 ) * 16*U*y.*z.*(H-y).*(H-z)/H^4.*(x==0)+ 0.*x.*y); 
data.bcDir{2} = @(x, y, z, t, param)(0.*x.*y); 
data.bcDir{3} = @(x, y, z, t, param)(0.*x.*y); 

% Neumann
data.bcNeu{1} = @(x, y, z, t, param)(0.*x.*y);
data.bcNeu{2} = @(x, y, z, t, param)(0.*x.*y);
data.bcNeu{3} = @(x, y, z, t, param)(0.*x.*y);

Set the inital condition on the velocity (zero in this case)

% initial condition
data.u0{1}    = @(x, y, z, t, param)(0.*x.*y);
data.u0{2}    = @(x, y, z, t, param)(0.*x.*y);
data.u0{3}    = @(x, y, z, t, param)(0.*x.*y);

Set physical parameters $\mu$ and $\rho$

% Model parameters
data.dynamic_viscosity   = 1e-3;
data.density             = 1;

Set tolerance and maximum number of iterations of the nonlinear solver (Newton)

% Nonlinear solver
data.NonlinearSolver.tol         = 1e-6; 
data.NonlinearSolver.maxit       = 30;

Set stabilization to be used (only ‘SUPG’ available for the time being)

% Stabilization
data.Stabilization = 'SUPG';

Set parameters for the linear solver and preconditioner. In this case, if a parallel pool is open we use gmres with a one-level Additive Schwarz preconditioner; otherwise we use a direct solver

% detect parallel pool
poolobj = gcp('nocreate');
if isempty(poolobj)
    poolsize = 0;
else
    poolsize = poolobj.NumWorkers;
end

if poolsize > 0
    % Linear Solver
    data.LinearSolver.type              = 'gmres'; % MUMPS, backslash, gmres
    data.LinearSolver.tol               = 1e-8;
    data.LinearSolver.maxit             = 500;
    data.LinearSolver.gmres_verbosity   = 5;
    
    % Preconditioner
    data.Preconditioner.type              = 'AdditiveSchwarz'; % AdditiveSchwarz, None, ILU
    data.Preconditioner.local_solver      = 'matlab_lu'; % matlab_lu, MUMPS
    data.Preconditioner.overlap_level     = 2;
    data.Preconditioner.mumps_reordering  = 7;
    data.Preconditioner.num_subdomains    = poolsize; %poolsize, number of subdomains

else
    % Linear Solver
    data.LinearSolver.type              = 'backslash'; % MUMPS, backslash, gmres
    data.LinearSolver.mumps_reordering  = 7;

    % Preconditioner
    data.Preconditioner.type         = 'None'; % AdditiveSchwarz, None, ILU
end

Set time parameters: BDF order, initial and final time, timestep and convective term treatment

% time 
data.time.BDF_order  = 2;
data.time.t0         = 0;
data.time.dt         = 0.005;
data.time.tf         = 2;
data.time.nonlinearity  = 'implicit';% 'semi-implicit'

Set output options. In this case we require to compute and export to a txt file the aerodynamic forces (i.e. drag and lift coefficients) computed on the boundary with flag 6

%% Output options
data.Output.DragLift.computeDragLift = 1;
data.Output.DragLift.factor          = 2/(0.1*1^2*0.41);
data.Output.DragLift.flag            = 6;
data.Output.DragLift.filename        = 'Results/AerodynamicForces_Re100.txt';

It also possible to export to a txt file the flow rate computed on selected boundaries (e.g. those with flag 2,3,4,5,6,7,200 as below)

%% Output options
data.Output.FlowRates.computeFlowRates = 1;
data.Output.FlowRates.flag             = [2:7 200];
data.Output.FlowRates.filename         = 'Results/FlowRates.txt';

Resistance boundary conditions

A family of RC type outflow boundary conditions for cardiovascular simulations are available. Please note that they are supported only in 3D in combination with the implicit solver. They are formulated as follows:

$\begin{equation} \label{eq:RC_BC} p + RC \frac{dp}{dt} = RQ + p_0 \qquad \text{on } \Gamma_o \end{equation}$

where $p$ is the average pressure on the outflow boundary $\Gamma_o$, $R$ is the resistance, $C$ is the capacitance, $Q = \int_{\Gamma_o} \mathbf{v} \cdot \mathbf{n} \,d\Gamma$ is the flow rate, and $p_0$ sets the physiological values of intramural blood pressure. From the algorithmic standpoint, equation \eqref{eq:RC_BC} is discretized in time using a semi-implicit approach:

$\begin{equation} p^{n+1} + RC \frac{p^{n+1} - p^n}{\Delta t} = RQ^n + p_0 \qquad \text{on } \Gamma_o, \end{equation}$

As a result, this condition is imposed by setting $\boldsymbol{\sigma}^{n+1}\mathbf{n} = -p^{n+1}\mathbf{n}$ on $\Gamma_o$, with

$\begin{equation} \label{eq:RC_BC3} p^{n+1} = \frac{\Delta t}{\Delta t + RC} \left( Q^n R + p^n \frac{R C}{\Delta t} + p_0 \right). \end{equation}$

If $C=0$, \eqref{eq:RC_BC3} reduces to the following (explicit) resistance BC

$p^{n+1} = Q^n R + p_0.$

In the datafile, this condition can be imposed by setting the following fields:

data.flag_resistance     = [2 4:7];
data.OutFlow_Resistance  = [178.39 1120.66 1008.25 1311.84 674.94];
data.OutFlow_RefPressure = @(t) 80 * 1333;% corresponds to 80 mmHg
data.OutFlow_Capacity    = [7.817 1.234 1.373 1.050 2.037]*1e-4;

In Fluid-Structure Interaction simulations, if the values of the resistance coefficients are not available, they can be estimated using a first order Taylor approximation of the formula for absorbing BCs:

$R \approx \frac{\sqrt{\rho \beta}}{\sqrt{2}A_f^{5/4}} \quad \text{with} \quad \beta = \frac{\sqrt{\pi} h_s R_f E}{1-\nu^2},$

where $E$ and $\nu$ are the arterial wall Young modulus and Poisson coefficient, respectively, $A_f$ is the area of the fluid outflow section, $R_f$ its radius, and $h_s$ is such that the wall thickness equals $h_s R_f$. In the datafile, it is sufficient to set the following fields:

data.OutFlow_ComputeApproximateResistance = true;
data.OutFlow_Poisson                      = 0.45;
data.OutFlow_Young                        = 1e+07;
data.OutFlow_WallThicknessPercentage      = 0.12;% wall thickness is 12% of the lumen radius

Main file

Create a subfolder Results where to save the txt with the values of the aerodynamic forces and a subfolder Figures where to save the vtk files

[~,~,~]  = mkdir('Results');
[~,~,~]  = mkdir('Figures');

Set the spatial dimension of the problem $d = 2,3$ and load the mesh fluid_L1.msh saved in the mesh subfolder

%% load P1 mesh 
dim      =  3; % 3D simulation
[vertices, boundaries, elements] = msh_to_Mmesh('mesh/fluid_L1', dim);

Set the finite element spaces for velocity and pressure: either $\mathbb{P}_2$-$\mathbb{P}_1$, $\mathbb{P}_1^{\text{bubble}}$-$\mathbb{P}_1$ or $\mathbb{P}_1$-$\mathbb{P}_1$

%% specify FEM approximation (velocity and pressure)
fem        = {'P2', 'P1'};% {'B1', 'P1'}, {'P1', 'P1'};

Launch the Navier-Stokes solver. At each timestep, the solution is exported to a vtk file with rootname SolP2_

%% solve  
[U, MESH, DATA] = NSt_Solver(dim, elements, vertices, boundaries, fem, 'datafileName', [], 'Figures/SolP2_');

Load the drag and lift coefficients and do some postprocessing

%% Postprocessing

RefDrag = 6.185331;
RefLift = 9.40135e-3;

[t, Drag, Lift] = importAerodynamicForces( DATA.Output.DragLift.filename );

fprintf('\nredbKIT RESULT:\n Drag = %1.3f, Lift = %1.5e\n',  Drag(end), Lift(end))
fprintf('\nREFERENCE VALUES:\n Drag = %1.3f, Lift = %1.5e\n',   RefDrag,  RefLift)

handle = figure;
subplot(2,1,1)
plot(t,Drag,'-b','LineWidth',2)
hold on
plot(t,Drag*0+RefDrag,'--k','LineWidth',2)
legend('Drag Coefficient','Reference Value')
xlabel('Time [s]')
grid on

subplot(2,1,2)
plot(t,Lift,'-r','LineWidth',2)
hold on
plot(t,Lift*0+RefLift,'--k','LineWidth',2)
legend('Lift Coefficient','Reference Value')
xlabel('Time [s]')
grid on

saveas(handle,'Results/Re20_DragLift','epsc');
saveas(handle,'Results/Re20_DragLift','fig');