Feel++ Book Contributors https://github.com/feelpp/book.feelpp.org/graphs/contributors
Benchmarks are available at CFD Benchmarks. 
Feel++ Book Contributors https://github.com/feelpp/book.feelpp.org/graphs/contributors
Feel++ Book Contributors https://github.com/feelpp/book.feelpp.org/graphs/contributors
Benchmarks are available at CFD Benchmarks. 
Notation  Quantity  Unit 

\(\rho_f\) 
fluid density 
\(kg \cdot m^{3}\) 
\(\boldsymbol{u}_f\) 
fluid velocity 
\(m \cdot s^{1}\) 
\(\boldsymbol{\sigma}_f\) 
fluid stress tensor 
\(N \cdot m^{2}\) 
\(\boldsymbol{f}^t_f\) 
source term 
\(kg \cdot m^{3} \cdot s^{1}\) 
\(p_f\) 
pressure fields 
\(kg \cdot m^{1} \cdot s^{2}\) 
\(\mu_f\) 
dynamic viscosity 
\(kg \cdot m^{1} \cdot s^{1}\) 
\(\bar{U}\) 
characteristic inflow velocity 
\(m \cdot s^{1}\) 
\(\nu\) 
kinematic viscosity 
\(m^2 \cdot s^{1}\) 
\(L\) 
characteristic length 
\(m\) 
NavierStokes model is used to model incompressible Newtonian fluid. It can be described by these conservative laws :
we complete this set of equations with the fluid constitutive law
\boldsymbol{\sigma}_{f} = p_f \boldsymbol{I} + 2\mu_f D(\boldsymbol{u}_{f})
with strain tensor \(D(\boldsymbol{u}_{f})\) defined by :
D(\boldsymbol{u}_{f}) = \frac{1}{2} (\nabla_\mathrm{x} \boldsymbol{u}_f + (\nabla_\mathrm{x} \boldsymbol{u}_f)^T)
An alternative model is the Stokes model. It is valid in the case of small Reynolds number. It corresponds to the same formulation than NavierStokes equations but without the convective term \(\left( \boldsymbol{u}_{f} \cdot \nabla_{\mathrm{x}} \right) \boldsymbol{u}_{f}\) .
Let’s define a bounded domain \(\Omega \subset \mathbb{R}^p\) (\(p=2,3\)) decomposed into two subdomains \(\Omega_1\) and \(\Omega_2\). We denote \(\Gamma\) the interface between the two partitions. The goal of the level set method is to track implicitly the interface \(\Gamma(t)\) moving at a velocity \(\mathbf{u}\). The level set method has been described in [osher] and its main ingredient is a continuous scalar function \(\phi\) (the /level set/ function) defined on the whole domain. This function is chosen to be positive in \(\Omega_1\), negative in \(\Omega_2\) and zero on \(\Gamma\). The motion of the interface is then described by the advection of the level set function with a divergence free velocity field \(\mathbf{u}\):
A convenient choice for \(\phi\) is a signed distance function to the interface. Indeed, the property \(\nabla \phi = 1\) of distance functions eases the numerical solution and gives a convenient support for delta and Heaviside functions).
Nevertheless, it is known that the advection equation \eqref{eq:advection} does not conserve the property \(\nabla \phi=1\). Thus, when \(\nabla \phi\) is far from \(1\) we have to reset \(\phi\) as a distance function without moving the interface. To do so we can either use an HamiltonJacobi method or the fast marching method (see \cite{Winkelmann2007} for details about the fast marching method).
In twofluid flow simulations, we need to define some quantities related to the interface such as the density, the viscosity, or some interface forces. To this end, we introduce the smoothed Heaviside and delta functions:
where \(\varepsilon\) is a parameter defining a ``numerical thickness'' of the interface. A typical value of \(\varepsilon\) is \(1.5 h\) where \(h\) is the mesh size of elements crossed by the isovalue \(0\) of the level set function.
The Heaviside function is used to define parameters having different values on each subdomains. For example, we define the density of twofluid flow as \(\rho = \rho_2 + (\rho_1\rho_2) H_{\varepsilon}(\phi)\) (we use a similar expression for the viscosity \(\nu\)). Regarding the delta function, it is used to define quantities on the interface. In particular, in the variational formulations, we replace integrals over the interface \(\Gamma\) by integrals over the entire domain \(\Omega\) using the smoothed delta function. If \(\phi\) is a signed distance function, we have : \(\int_{\Gamma} 1 \simeq \int_{\Omega} \delta_{\varepsilon}(\phi)\). If \(\phi\) is not close enough to a distance function, then \(\int_{\Gamma} 1 \simeq \int_{\Omega} \nabla \phi \delta_{\varepsilon}(\phi)\) which still tends to the measure of \(\Gamma\) as \(\varepsilon\) vanishes. However, if \(\phi\) is not a distance function, the support of \(\delta_{\varepsilon}\) can have a different size on each side of the interface. More precisely, the support of \(\delta_{\varepsilon}\) is narrowed on the side where \(\nabla \phi>1\) and enlarged on regions where \(\nabla \phi<1\). It has been shown in [cottet] that replacing \(\phi\) by \(\frac{\phi}{\nabla \phi}\) has the property that \(\nabla \frac{\phi}{\nabla \phi} \simeq 1\) near the interface and has the same isovalue \(0\) as \(\phi\). Thus, replacing \(\phi\) by \(\frac{\phi}{\nabla \phi}\) as support of the delta function does not move the interface. Moreover, the spread interface has the same size on each part of the levelset \(\phi=0\). It reads \(\int_{\Gamma} 1 \simeq \int_{\Omega} \delta_{\varepsilon}(\frac{\phi}{ \nabla \phi})\). The same technique is used for the Heaviside function.
We use Feel++ to discretize and solve the problem. Equation \eqref{eq:advection} is solved using a stabilized finite element method. We have implemented several stabilization methods such as Streamline Upwind Diffusion (SUPG), Galerkin Least Square (GLS) and Subgrid Scale (SGS). A general review of these methods is available in [Franca92]. Other available methods include the Continuous Interior Penalty method (CIP) are described in \cite{Burman2006}. The variational formulation at the semidiscrete level for the stabilized equation \eqref{eq:advection} reads, find \(\phi_h \in {\mathbb R}_h^k\) such that \(\forall \psi_h \in {\mathbb R}_h^k\) :
where \(S(\cdot, \cdot)\) stands for the stabilization bilinear form (see section \ref{sec:membrinext} for description of \({\mathbb R}_h^k\) and \(\mathbf{u}_h\)). In our case, we use a BDF2 scheme which needs the solution at the two previous time step to compute the one at present time. For the first time step computation or after a reinitialization we use an Euler scheme.
The twofluid flow problem can be expressed by
Where \(\boldsymbol{f}_\phi\) is the force obtain by projection of the density of interfacial forces on the domain \(\Omega\)
The CFD Toolbox supports both the Stokes and the incompressible NavierStokes equations.
The fluid mechanics model (NavierStokes
or Stokes
) can be selected in json file:
"Model": "NavierStokes"
The next step is to define the fluid material by setting its properties namely the density \(\rho_f\) and viscosity \(\mu_f\). In next table, we find the correspondance between the mathematical names and the json names.
Parameter  Symbol 

\(\mu_f\) 

\(\rho_f\) 

A Materials
section is introduced in json file in order to configure the fluid properties. For each mesh marker, we can define the material properties associated.
"Materials":
{
"<marker>"
{
"name":"water",
"rho":"1.0e3",
"mu":"1.0"
}
}
The non Newtonian properties are defined in cfg file in fluid section.
The viscosity law is activated by:
option  values 

viscosity.law 
newtonian, power_law, walburnschneck_law, carreau_law, carreauyasuda_law 
Then, each model are configured with the options reported in the following table:
Viscosity law  options  unit 

power_law 
power_law.k power_law.n 
dimensionless dimensionless 
walburnschneck_law 
hematocrit TPMA walburnschneck_law.C1 walburnschneck_law.C2 walburnschneck_law.C3 walburnschneck_law.C4 
Percentage g/l dimensionless dimensionless dimensionless l/g 
carreau_law 
viscosity.zero_shear viscosity.infinite_shear carreau_law.lambda carreau_law.n 
\(kg.m^{1}.s^{1}\) dimensionless dimensionless 
carreauyasuda_law 
viscosity.zero_shear viscosity.infinite_shear carreauyasuda_law.lambda carreauyasuda_law.n carreauyasuda_law.a 
\(kg/(m \times s)\) \(kg/(m \times s)\) dimensionless dimensionless dimensionless 
We start by a listing of boundary conditions supported in fluid mechanics model.
A Dirichlet condition on velocity field reads:
\boldsymbol{u}_f = \boldsymbol{g} \quad \text{ on } \Gamma
or only a component of vector \(\boldsymbol{u}_f =(u_f^1,u_f^2 ,u_f^3 )\)
u_f^i = g \quad \text{ on } \Gamma
Several methods are available to enforce the boundary condition:
elimination
Nitsche
Lagrange multiplier
p & = g \\ \boldsymbol{u}_f \times {\boldsymbol{ n }} & = \boldsymbol{0}
Name  Expression 

Neumann_scalar 
\(\boldsymbol{\sigma}_{f} \boldsymbol{n} = g \ \boldsymbol{n} \) 
Neumann_vectorial 
\(\boldsymbol{\sigma}_{f} \boldsymbol{n} = \boldsymbol{g} \) 
Neumann_tensor2 
\(\boldsymbol{\sigma}_{f} \boldsymbol{n} = g \ \boldsymbol{n}\) 
\boldsymbol{u}_f \cdot \boldsymbol{ n } = 0
The boundary condition at inlets allow to define a velocity profile on a set of marked faces \(\Gamma_{\mathrm{inlet}}\) in fluid mesh:
\boldsymbol{u}_f =  g \ \boldsymbol{ n } \quad \text{ on } \Gamma_{\mathrm{inlet}}
The function \(g\) is computed from flow velocity profiles:
\text{Find } g \in C^0(\Gamma) \text{ such that } \\ \begin{eqnarray} g &=& \beta \quad &\text{ in } \Gamma \setminus \partial\Gamma \\ g&=&0 \quad &\text{ on } \partial\Gamma \end{eqnarray}
\text{Find } g \in H^2(\Gamma) \text{ such that : } \\ \begin{eqnarray} \Delta g &=& \beta \quad &\text{ in } \Gamma \\ g&=&0 \quad &\text{ on } \partial\Gamma \end{eqnarray}
where \(\beta\) is a constant determined by adding a constraint to the inflow:
\(\max( g ) = \alpha \)
\(\int_\Gamma ( g \ \boldsymbol{n} ) \cdot \boldsymbol{n} = \alpha\)
Option  Value  Default value  Description 

shape 

select a shape profile for inflow 

constraint 

give a constraint wich controle velocity 

expr 
string 
symbolic expression of constraint value 
Option  Value  Default value  Description 

model 
free,windkessel 
free 
select an outlet modeling 
\boldsymbol{\sigma}_{f} \boldsymbol{ n } = \boldsymbol{0}
We use a 3element Windkessel model for modeling an outflow boundary condition. We define \(P_l\) a pressure and \(Q_l\) the flow rate. The outflow model is discribed by the following system of differential equations:
\left\{ \begin{aligned} C_{d,l} \frac{\partial \pi_l}{\partial t} + \frac{\pi_l}{R_{d,l}} = Q_l \\ P_l = R_{p,l} Q_l + \pi_l \end{aligned} \right.
Coefficients \(R_{p,l}\) and \(R_{d,l}\) represent respectively the proximal and distal resistance. The constant \(C_{d,l}\) is the capacitance of blood vessel. The unknowns \(P_l\) and \(\pi_l\) are called proximal pressure and distal pressure. Then we define the coupling between this outflow model and the fluid model by these two relationships:
\begin{align} Q_l &= \int_{\Gamma_l} \boldsymbol{u}_f \cdot \boldsymbol{ n }_f \\ \boldsymbol{\sigma}_f \boldsymbol{ n }_f &= P_l \boldsymbol{ n }_f \end{align}
Option  Value  Description 

windkessel_coupling 
explicit, implicit 
coupling type with the NavierStokes equation 
windkessel_Rd 
real 
distal resistance 
windkessel_Rp 
real 
proximal resistance 
windkessel_Cd 
real 
capacitance 
Boundary conditions are set in the json files in the category BoundaryConditions
.
Then <field>
and <bc_type>
are chosen from type of boundary condition.
The parameter <marker>
corresponds to mesh marker where the boundary condition is applied.
Finally, we define some specific options inside a marker.
"BoundaryConditions":
{
"<field>":
{
"<bc_type>":
{
"<marker>":
{
"<option1>":"<value1>",
"<option2>":"<value2>",
// ...
}
}
}
}
Field  Name  Option  Entity 

velocity 
Dirichlet 
expr type number alemesh_bc 
faces, edges, points 
velocity_x velocity_y velocity_z 
Dirichlet 
expr type number alemesh_bc 
faces, edges, points 
velocity 
Neumann_scalar 
expr number alemesh_bc 
faces 
velocity 
Neumann_vectorial 
expr number alemesh_bc 
faces 
velocity 
Neumann_tensor2 
expr number alemesh_bc 
faces 
velocity 
slip 
alemesh_bc 
faces 
pressure 
Dirichlet 
expr number alemesh_bc 
faces 
fluid 
outlet 
number alemesh_bc model windkessel_coupling windkessel_Rd windkessel_Rp windkessel_Cd 
faces 
fluid 
inlet 
expr shape constraint number alemesh_bc 
faces 
Body forces are also defined in BoundaryConditions
category in json file.
"VolumicForces":
{
"<marker>":
{
"expr":"{0,0,gravityCst*7850}:gravityCst"
}
}
The marker corresponds to mesh elements marked with this tag. If the marker is an empty string, it corresponds to all elements of the mesh.
"PostProcess":
{
"Fields":["field1","field2",...],
"Measures":
{
"<measure type>":
{
"label":
{
"<range type>":"value",
"fields":["field1","field3"]
}
}
}
}
The fields allowed to be exported in the Fields
section are:
velocity
pressure
displacement
vorticity
stress or normalstress
wallshearstress
density
viscosity
pid
alemesh
Points
Force
FlowRate
Pressure
VelocityDivergence
In order to evaluate velocity or pressure at specific points and save the results in .csv file, the user must define:
"<tag>" representing this data in the .csv file
the coordinate of point
the fields evaluated ("pressure" or "velocity")
"Points":
{
"<tag>":
{
"coord":"{0.6,0.2,0}",
"fields":"pressure"
},
"<tag>":
{
"coord":"{0.15,0.2,0}",
"fields":"velocity"
}
}
The flow rate can be evaluated and save on .csv file. The user must define:
"<tag>" representing this data in the .csv file
"<face_marker>" representing the name of marked face
the fluid direction ("interior_normal" or "exterior_normal") of the flow rate.
"FlowRate":
{
"<tag>":
{
"markers":"<face_marker>",
"direction":"interior_normal"
},
"<tag>":
{
"markers":"<face_marker>",
"direction":"exterior_normal"
}
}
A function defined by a symbolic expression can be represented as a finite element field thanks to nodal projection. This function can be exported.
"Functions":
{
"toto":{ "expr":"x*y:x:y"},
"toto2":{ "expr":"0.5*ubar*x*y:x:y:ubar"},
"totoV":{ "expr":"{2*x,y}:x:y"}
},
"PostProcess":
{
"Fields":["velocity","pressure","pid","totoV","toto","toto2"],
}
Galerkin leatSquare (GLS) stabilization methods can be activated from the cfg file by adding stabilizationgls=1
in the fluid
prefix.
Others options available are enumerated in the next table and must be given with the prefix fluid.stabilizationgls
.
Option  Value  Default value  Description 

type 


type of stabilization 
parameter.method 


method used to compute the stabilization parameter 
parameter.hsize.method 


method used for evalute an element mesh size 
parameter.eigenvalue.penallambdaK 
real 
0. 
add a mass matrix scaled by this factor in the eigen value problem for the stabilization parameter 
convectiondiffusion.location.expressions 
string 
if given, the stabilization is apply only on mesh location which verify 
Documentation pending 
programme avalaible :
feelpp_toolbox_fluid_2d
feelpp_toolbox_fluid_3d
mpirun np 4 feelpp_toolbox_fluid_2d configfile=<myfile.cfg>