CFD Model
1. Notations and units
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\) 
2. Equations
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}\) .
3. CFD Toolbox
3.1. Models
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"
3.2. Materials
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"
}
}
3.2.1. Generalised Newtonian fluid
A non newtonian fluid is characterized by a non constant viscosity, which is a function of strain rate \(\boldsymbol{D}\left(\boldsymbol{u}_{f}\right)\).
We start by introducing a metric of the rate of deformation, denoted by \(\dot{\gamma}\):
\[ \dot{\gamma} \equiv \sqrt{2 \ tr \left( \boldsymbol{D}\left(\boldsymbol{u}_{f}\right)^{2} \right) } \] 
We represent the viscosity \(\mu_f\) as a function of \(\dot{\gamma}\). The deviatoric stress tensor \(\boldsymbol{\tau}\) is obtained thanks to generalised Newtonian model, which takes the following form:
\[ \boldsymbol{\tau} = 2 \mu_f \left(\dot{\gamma} \right) \boldsymbol{D}\left(\boldsymbol{u}_{f}\right) \] 
The simplest example of a generalised Newtonian model is the powerlaw fluid, which has a viscosity given by:
\[ \mu_f \left(\dot{\gamma} \right) = k \dot{\gamma}^{n1} \] 
where \(k\) and \(n < 1\) are two parameters related to fluid properties.
Blood flow viscosity
In the context of blood flow modeling, an extension of the power model was proposed by Walburn and Schneck. The parameters \(k\) and \(n\) are related to the hematocrit \(Ht\) and Total Proteins Minus Albumin (TPMA) as follows

and \(C_i, i=1,..,4\) are parameters to fit with experimental data.
Another family of generalised Newtonian model can be defined from a function \(\Phi\) express by:
\[ \Phi\left( \dot{\gamma}, \mu_{\infty},\mu_{0} \right) = \frac{\mu\left(\dot{\gamma}\right)  \mu_{\infty}}{\mu_{0}\mu_{\infty}} \] 
where \(\mu_0\) and \(\mu_{\infty}\) are the asymptotic viscosities at zero and infinite shear rate.
Viscosity law  \(\Phi\left( \dot{\gamma}, \mu_{\infty},\mu_{0} \right) \) 

Carreau 
\(\left(1+\left(\lambda\dot{\gamma}\right)^{2}\right)^{(n1)/2}\) 
CarreauYasuda 
\(\left(1+\left(\lambda\dot{\gamma}\right)^{a}\right)^{(n1)/a}\) 
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 
3.3. Boundary Conditions
We start by a listing of boundary conditions supported in fluid mechanics model.
3.3.1. Dirichlet on velocity
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
3.3.2. Dirichlet on pressure
\[\begin{split} p & = g \\ \boldsymbol{u}_f \times {\boldsymbol{ n }} & = \boldsymbol{0} \end{split}\] 
3.3.3. Neumann
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}\) 
3.3.4. Slip
\[ \boldsymbol{u}_f \cdot \boldsymbol{ n } = 0 \] 
3.3.5. Inlet
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:
 Constant profile
\[ \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} \] 
 Parabolic profile
\[ \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:
 velocity_max

\(\max( g ) = \alpha \)
 flow_rate

\(\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 
3.3.6. Outlet flow
Option  Value  Default value  Description 

model 
free,windkessel 
free 
select an outlet modeling 
Free outflow
\[ \boldsymbol{\sigma}_{f} \boldsymbol{ n } = \boldsymbol{0} \] 
Windkessel model
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 
3.3.7. Implementation of boundary conditions in json
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>",
// ...
}
}
}
}
3.3.8. Options summary
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 
3.4. Body forces
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.
3.5. Post Processing
"PostProcess":
{
"Fields":["field1","field2",...],
"Measures":
{
"<measure type>":
{
"label":
{
"<range type>":"value",
"fields":["field1","field3"]
}
}
}
}
Exports for vizualisation
The fields allowed to be exported in the Fields
section are:

velocity

pressure

displacement

vorticity

stress or normalstress

wallshearstress

density

viscosity

pid

alemesh
Measures

Points

Force

FlowRate

Pressure

VelocityDivergence
Points
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"
}
}
Flow rate
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"
}
}
Export user functions
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"],
}
3.6. Action
Let’s finish with a simple example in order to show how this works. We will interest us to a fluid flow into a cavity in 3D.
Feel++ code
Here is the code
First at all, we define our model type with
typedef FeelModels::FluidMechanics< Simplex<FEELPP_DIM,1>, Lagrange<OrderVelocity,Vectorial,Continuous,PointSetFekete>, Lagrange<OrderPressure,Scalar,Continuous,PointSetFekete> > model_type;
We choose here a \(\mathbb{P}_2\) space for the velocity order and \(\mathbb{P}_1\) space for the pressure order. This definition allows us to create our fluid model object FM like this
auto FM = model_type::New("fluid");
The method New
retrieves all data from the configuration and json files and build a mesh if needed.
With this object, we can initialize our model parameters, such as velocity or boundaries conditions. Data on our model and on the numeric solver are then save and print on the terminal. This is made by
FM>init(); FM>printAndSaveInfo();
Now that our model is completed, we can solve the associated problem. To begin the resolution
FM>isStationary()
determine if our model is stationary or not.
If it is, then we need to solve our system only one time and export the obtained results.
FM>solve(); FM>exportResults();
If it’s not, our model is time reliant, and a loop on time is necessary. Our model is then solve and the results are export at each time step.
for ( ; !FM>timeStepBase()>isFinished(); FM>updateTimeStep() ) { FM>solve(); FM>exportResults(); }
Code
{% include "../Examples/fluid_model.cpp" %}
Configuration file
The config file is used to define options, linked to our case, we would have the possibility to change at will. It can be, for example, files paths as follows
[fluid] geofile=$cfgdir/cavity3d.geo filename=$cfgdir/cavity3d.json [exporter] directory=applications/models/fluid/cavity3d/$fluid_tag
It can also be resolution dependent parameters such as mesh elements size, methods used to define our problem and solvers.
[fluid] solver=Oseen #Picard,Newton linearsystemcstupdate=false jacobianlinearupdate=false snesmonitor=true snesmaxit=100 snesmaxitreuse=100 sneskspmaxit=1000 sneskspmaxitreuse=100 pctype=lu #gasm,lu,fieldsplit,ilu
In this case, we use Oseen to define our problem, we set the update of linear system constant and jacobian linear as "no update", we discretize values associated to SNES ( Scalable Nonlinear Equations Solvers ), and finally we choose LU as the preconditioner method.
Code
{% include "../Examples/cavity3d.cfg" %}
Json file
First at all, we define some general information like the name ( and short name ) and the model we would like to use
"Name": "Fluid Mechanics",
"ShortName":"Fluid",
"Model":"NavierStokes",
Then we define the material properties. In our case, the fluid, define by rho` its density in \(kg\cdot m^{3}\) and mu
its dynamic viscosity in \(kg\cdot (m \cdot s)^{1}\), is the only material we have to define.
"Materials":
{
"Fluid":{
"name":"myFluidMat",
"rho":"1.0",
"mu":"0.01"
}
},
The boundary conditions are the next aspect we define. Here, we impose on the velocity \(u_f\) Dirichlet conditions at two specific places : lid
and wall
.
"BoundaryConditions":
{
"velocity":
{
"Dirichlet":
{
"lid":
{
"expr":"{ 1,0,0}:x:y:z"
},
"wall":
{
"expr":"{0,0,0}"
}
}
}
}
The post process aspect is the last one to define. We choose the fields we want to export ( velocity, pressure and pid ). Furthermore, we want to measure forces on wall
and the pressure at point \(A\).
"PostProcess":
{
"Fields":["velocity","pressure","pid"],
"Measures":
{
"Forces":"wall",
"Points":
{
"pointA":
{
"coord":"{0.5,0.5,0.5}",
"fields":"pressure"
}
}
}
}
}
Code
{% include "../Examples/cavity3d.json" %}
Compilation/Execution
Once you’ve a build dir, you just have to realise the command make
at
{buildir}/applications/models/fluid
This will generate executables named feelpp_application_fluid_*
. To execute it, you need to give the path of the cfg file associated to your case, with configfile
.
For example
./feelpp_application_fluid_3d configfile={sourcedir}/applications/models/fluid/cavity/cavity3d.cfg
is how to execute the case ahead.
The result files are then stored by default in
feel/applications/models/fluid/{case_name}/ {velocity_space}{pression_space}{Geometric_order}/{processor_used}
If we return once again at our example, the result files are in
feel/applications/models/fluid/cavity3d/P2P1G1/np_1
Examples
4. Turek & Hron CFD Benchmark
4.1. Introduction
We implement the benchmark proposed by Turek and Hron, on the behavior of drag and lift forces of a flow around an object composed by a pole and a bar, see Figure [imggeometry1].
The software and the numerical results were initially obtained from Vincent Chabannes.
This benchmark is linked to the Turek Hron CSM and Turek Hron FSI benchmarks. 
4.2. Problem Description
We consider a 2D model representative of a laminar incompressible flow around an obstacle. The flow domain, named \(\Omega_f\), is contained into the rectangle \( \lbrack 0,2.5 \rbrack \times \lbrack 0,0.41 \rbrack \). It is characterised, in particular, by its dynamic viscosity \(\mu_f\) and by its density \(\rho_f\). In this case, the fluid material we used is glycerine.
In order to describe the flow, the incompressible NavierStokes model is chosen for this case, define by the conservation of momentum equation and the conservation of mass equation. At them, we add the material constitutive equation, that help us to define \(\boldsymbol{\sigma}_f\)
The goal of this benchmark is to study the behavior of lift forces \(F_L\) and drag forces \(F_D\), with three different fluid dynamics applied on the obstacle, i.e on \(\Gamma_{obst}\), we made rigid by setting specific structure parameters. To simulate these cases, different mean inflow velocities, and thus different Reynolds numbers, will be used.
4.2.1. Boundary conditions
We set

on \(\Gamma_{in}\), an inflow Dirichlet condition : \( \boldsymbol{u}_f=(v_{in},0) \)

on \(\Gamma_{wall}\) and \(\Gamma_{obst}\), a homogeneous Dirichlet condition : \( \boldsymbol{u}_f=\boldsymbol{0} \)

on \(\Gamma_{out}\), a Neumann condition : \( \boldsymbol{\sigma}_f\boldsymbol{ n }_f=\boldsymbol{0} \)
4.2.2. Initial conditions
We use a parabolic velocity profile, in order to describe the flow inlet by \( \Gamma_{in} \), which can be express by
\[ v_{cst} = 1.5 \bar{U} \frac{4}{0.1681} y \left(0.41y\right) \] 
where \(\bar{U}\) is the mean inflow velocity.
However, we want to impose a progressive increase of this velocity profile. That’s why we define
\[ v_{in} = \left\{ \begin{aligned} & v_{cst} \frac{1\cos\left( \frac{\pi}{2} t \right) }{2} \quad & \text{ if } t < 2 \\ & v_{cst} \quad & \text{ otherwise } \end{aligned} \right. \] 
With t the time.
Moreover, in this case, there is no source term, so \(f_f\equiv 0\).
4.3. Inputs
The following table displays the various fixed and variables parameters of this testcase.
Name  Description  Nominal Value  Units  

\(l\) 
elastic structure length 
\(0.35\) 
\(m\) 

\(h\) 
elastic structure height 
\(0.02\) 
\(m\) 

\(r\) 
cylinder radius 
\(0.05\) 
\(m\) 

\(C\) 
cylinder center coordinates 
\((0.2,0.2)\) 
\(m\) 

\(\nu_f\) 
kinematic viscosity 
\(1\times 10^{3}\) 
\(m^2/s\) 

\(\mu_f\) 
dynamic viscosity 
\(1\) 
\(kg/(m \times s)\) 

\(\rho_f\) 
density 
\(1000\) 
\(kg/m^3\) 

\(f_f\) 
source term 
0 
\(kg/(m^3 \times s)\) 

\(\bar{U}\) 
characteristic inflow velocity 

\(m/s\) 
4.4. Outputs
As defined above, the goal of this benchmark is to measure the drag and lift forces, \(F_D\) and \(F_L\), to control the fluid solver behavior. They can be obtain from
\[ (F_D,F_L)=\int_{\Gamma_{obst}}\boldsymbol{\sigma}_f \boldsymbol{ n }_f \] 
where \(\boldsymbol{n}_f\) the outer unit normal vector from \(\partial \Omega_f\).
4.5. Discretization
To realize these tests, we made the choice to used \(P_N\)\(P_{N1}\) TaylorHood finite elements, described by Chabannes, to discretize space. With the time discretization, we use BDF, for Backward Differentation Formulation, schemes at different orders \(\)q\(\).
4.5.1. Solvers
Here are the different solvers ( linear and nonlinear ) used during results acquisition.
type 
gmres 
relative tolerance 
1e13 
max iteration 
1000 
reuse preconditioner 
false 
relative tolerance 
1e8 
steps tolerance 
1e8 
max iteration 
CFD1/CFD2 : 100  CFD3 : 50 
max iteration with reuse 
CFD1/CFD2 : 100  CFD3 : 50 
reuse jacobian 
false 
reuse jacobian rebuild at first Newton step 
true 
relative tolerance 
1e5 
max iteration 
1000 
max iteration with reuse 
CFD1/CFD2 : 100  CFD3 : 1000 
reuse preconditioner 
false 
reuse preconditioner rebuild at first Newton step 
false 
type 
lu 
package 
mumps 
4.6. Running the model
The configuration files are in toolboxes/fluid/TurekHron
. The different cases are implemented in the corresponding .cfg
files e.g. cfd1.cfg
, cfd2.cfg
and cfd3.cfg
.
The command line in feelpptoolboxes docker reads
$ mpirun np 4 /usr/local/bin/feelpp_toolbox_fluid_2d configfile cfd1.cfg
The result files are then stored by default in
feel/applications/models/fluid/TurekHron/"case_name"/"velocity_space""pression_space""Geometric_order"/"processor_used"
For example, for CFD2 case executed on \(\)12\(\) processors, with a \(\)P_2\(\) velocity approximation space, a \(\)P_1\(\) pressure approximation space and a geometric order of \(\)1\(\), the path is
feel/toolboxes/fluid/TurekHron/cfd2/P2P1G1/np_12
4.7. Results
Here are results from the different cases studied in this benchmark.
4.7.1. CFD1
\(\mathbf{N_{geo}}\)  \(\mathbf{N_{elt}}\)  \(\mathbf{N_{dof}}\)  Drag  Lift 

Reference Turek and Hron 
14.29 
1.119 

1 
9874 
45533 (\(P_2/P_1\)) 
14.217 
1.116 
1 
38094 
173608 (\(P_2/P_1\)) 
14.253 
1.120 
1 
59586 
270867 (\(P_2/P_1\)) 
14.262 
1.119 
2 
7026 
78758 (\(P_3/P_2\)) 
14.263 
1.121 
2 
59650 
660518 (\(P_3/P_2\)) 
14.278 
1.119 
3 
7026 
146057 (\(P_4/P_3\)) 
14.270 
1.120 
3 
59650 
1228831 (\(P_4/P_3\)) 
14.280 
1.119 
All the files used for this case can be found in this rep [geo file, config file, json file]
4.7.2. CFD2
\(\mathbf{N_{geo}}\)  \(\mathbf{N_{elt}}\)  \(\mathbf{N_{dof}}\)  Drag  Lift 

Reference Turek and Hron 
136.7 
10.53 

1 
7020 
32510 (\(P_2/P_1\)) 
135.33 
10.364 
1 
38094 
173608 (\(P_2/P_1\)) 
136.39 
10.537 
1 
59586 
270867 (\(P_2/P_1\)) 
136.49 
10.531 
2 
7026 
78758 (\(P_3/P_2\)) 
136.67 
10.548 
2 
59650 
660518 (\(P_3/P_2\)) 
136.66 
10.532 
3 
7026 
146057 (\(P_4/P_3\)) 
136.65 
10.539 
3 
59650 
1228831 (\(P_4/P_3\)) 
136.66 
10.533 
All the files used for this case can be found in this rep [geo file, config file, json file]
4.7.3. CFD3
As CFD3 is timedependent ( from BDF use ), results will be expressed as
\[ mean ± amplitude [frequency] \] 
where

mean is the average of the min and max values at the last period of oscillations.
\[ mean=\frac{1}{2}(max+min) \] 

amplitude is the difference of the max and the min at the last oscillation.
\[ amplitude=\frac{1}{2}(maxmin) \] 

frequency can be obtain by Fourier analysis on periodic data and retrieve the lowest frequency or by the following formula, if we know the period time T.
\[ frequency=\frac{1}{T} \] 
\(\mathbf{\Delta t}\)  \(\mathbf{N_{geo}}\)  \(\mathbf{N_{elt}}\)  \(\mathbf{N_{dof}}\)  \(\mathbf{N_{bdf}}\)  Drag  Lift 

0.005 
Reference Turek and Hron 
439.45 ± 5.6183[4.3956] 
−11.893 ± 437.81[4.3956] 
0.01 
1 
8042 
37514 (\(P_2/P_1\)) 
2 
437.47 ± 5.3750[4.3457] 
9.786 ± 437.54[4.3457] 
2 
2334 
26706 (\(P_3/P_2\)) 
2 
439.27 ± 5.1620[4.3457] 
8.887 ± 429.06[4.3457] 

2 
7970 
89790 (\(P_2/P_2\)) 
2 
439.56 ± 5.2335[4.3457] 
11.719 ± 425.81[4.3457] 
0.005 
1 
3509 
39843\(\)(P_3/P_2)\(\) 
2 
438.24 ± 5.5375[4.3945] 
11.024 ± 433.90[4.3945] 
1 
8042 
90582 (\(P_3/P_2\)) 
2 
439.25 ± 5.6130[4.3945] 
10.988 ± 437.70[4.3945] 

2 
2334 
26706 (\(P_3/P_2\)) 
2 
439.49 ± 5.5985[4.3945] 
10.534 ± 441.02[4.3945] 

2 
7970 
89790 (\(P_3/P_2\)) 
2 
439.71 ± 5.6410[4.3945] 
11.375 ± 438.37[4.3945] 

3 
3499 
73440 (\(P_4/P_3\)) 
3 
439.93 ± 5.8072[4.3945] 
14.511 ± 440.96[4.3945] 

4 
2314 
78168 (\(P_5/P_4\)) 
2 
439.66 ± 5.6412[4.3945] 
11.329 ± 438.93[4.3945] 
0.002 
2 
7942 
89482 (stem:[P_3/P_2)\(\) 
2 
439.81 ± 5.7370[4.3945] 
13.730 ± 439.30[4.3945] 
3 
2340 
49389 (\(P_4/P_3\)) 
2 
440.03 ± 5.7321[4.3945] 
13.250 ± 439.64[4.3945] 

3 
2334 
49266 (\(P_4/P_3\)) 
3 
440.06 ± 5.7773[4.3945] 
14.092 ± 440.07[4.3945] 
All the files used for this case can be found in this rep [geo file, config file, json file].
4.8. Geometrical Order
Add a section on geometrical order. 
4.9. Conclusion
The reference results, Turek and Hron, have been obtained with a time step \(\Delta t=0.05\). When we compare our results, with the same step and \(\mathrm{BDF}_2\), we observe that they are in accordance with the reference results.
With a larger \(\Delta t\), a discrepancy is observed, in particular for the drag force. It can also be seen at the same time step, with a higher order \(\mathrm{BDF}_n\) ( e.g. \(\mathrm{BDF}_3\) ). This suggests that the couple \(\)\Delta t=0.05\(\) and \(\mathrm{BDF}_2\) isn’t enough accurate.
4.10. Bibliography

[TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow, Lecture Notes in Computational Science and Engineering, 2006.

[Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Universitée de Grenoble, 2013.
5. Multifluid flows
We introduce here the multifluid flows benchmarks and the associated methdology.
5.1. Introduction
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).
5.2. Interface related quantities
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.
5.3. Numerical implementation and coupling with the fluid solver
We use the finite element C library <<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.
5.4. Twofluid flow
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\)
5.5. 2D Drops Benchmark
This benchmark has been proposed and realised by Hysing. It allows us to verify our level set code, our NavierStokes solver and how they couple together.
Computer codes, used for the acquisition of results, are from Vincent Doyeux, with the use of Chabannes's NavierStokes code.
5.5.1. Problem Description
We want to simulate the rising of a 2D bubble in a Newtonian fluid. The bubble, made of a specific fluid, is placed into a second one, with a higher density. Like this, the bubble, due to its lowest density and by the action of gravity, rises.
The equations used to define fluid bubble rising in an other are the NavierStokes for the fluid and the advection one for the level set method. As for the bubble rising, two forces are defined :

The gravity force : \(\boldsymbol{f}_g=\rho_\phi\boldsymbol{g}\)

The surface tension force : \(\boldsymbol{f}_{st}=\int_\Gamma\sigma\kappa\boldsymbol{ n } \)
We denote \( \Omega\times\rbrack0,3\rbrack \) the interest domain with \( \Omega=(0,1)\times(0,2) \). \(\Omega\) can be decompose into \(\Omega_1\), the domain outside the bubble and \(\Omega_2\) the domain inside the bubble and \(\Gamma\) the interface between these two.
Durig this benchmark, we will study two different cases : the first one with a ellipsoidal bubble and the second one with a squirted bubble.
Boundary conditions

On the lateral walls, we imposed slip conditions

On the horizontal walls, no slip conditions are imposed : \(\boldsymbol{u}=0 \)
Initial conditions
In order to let the bubble rise, its density must be inferior to the density of the exterior fluid, so \(\rho_1>\rho_2\)
5.5.2. Inputs
The following table displays the various fixed and variables parameters of this testcase.
Name  Description  Nominal Value  Units 

\(\boldsymbol{g}\) 
gravity acceleration 
\((0,0.98)\) 
\(m/s^2\) 
\(l\) 
length domain 
\(1\) 
\(m\) 
\(h\) 
height domain 
\(2\) 
\(m\) 
\(r\) 
bubble radius 
\(0.25\) 
\(m\) 
\(B_c\) 
bubble center 
\((0.5,0.5)\) 
\(m\) 
5.5.3. Outputs
In the first place, the quantities we want to measure are \(X_c\) the position of the center of the mass of the bubble, the velocity of the center of the mass \(U_c\) and the circularity \(c\), define as the ratio between the perimeter of a circle and the perimeter of the bubble. They can be expressed by
After that, we interest us to quantitative points for comparison as \(c_{min}\), the minimum of the circularity and \(t_{c_{min}}\), the time needed to obtain this minimum, \(u_{c_{max}}\) and \(t_{u_{c_{max}}}\) the maximum velocity and the time to attain it, or \(y_c(t=3)\) the position of the bubble at the final time step. We add a second maximum velocity \(u_{max}\) and \(u_{c_{max_2}}\) and its time \(t_{u_{c_{max_2}}}\) for the second test on the squirted bubble.
5.5.4. Discretization
This is the parameters associate to the two cases, which interest us here.
Case 
\(\rho_1\) 
\(\rho_2\) 
\(\mu_1\) 
\(\mu_2\) 
\(\sigma\) 
Re 
\(E_0\) 
ellipsoidal bubble (1) 
1000 
100 
10 
1 
24.5 
35 
10 
squirted bubble (2) 
1000 
1 
10 
0.1 
1.96 
35 
125 
5.5.6. Results
Test 2
We describe the different quantitative results for the two studied cases.
h 
\(c_{min}\) 
\(t_{c_{min}}\) 
\(u_{c_{max}}\) 
\(t_{u_{c_{max}}}\) 
\(y_c(t=3)\) 
lower bound 
0.9011 
1.8750 
0.2417 
0.9213 
1.0799 
upper bound 
0.9013 
1.9041 
0.2421 
0.9313 
1.0817 
0.02 
0.8981 
1.925 
0.2400 
0.9280 
1.0787 
0.01 
0.8999 
1.9 
0.2410 
0.9252 
1.0812 
0.00875 
0.89998 
1.9 
0.2410 
0.9259 
1.0814 
0.0075 
0.9001 
1.9 
0.2412 
0.9251 
1.0812 
0.00625 
0.8981 
1.9 
0.2412 
0.9248 
1.0815 
h 
\(c_{min}\) 
\(t_{c_{min}}\) 
\(u_{c_{max_1}}\) 
\(t_{u_{c_{max_1}}}\) 
\(u_{c_{max_2}}\) 
\(t_{u_{c_{max_2}}}\) 
\(y_c(t=3)\) 
lower bound 
0.4647 
2.4004 
0.2502 
0.7281 
0.2393 
1.9844 
1.1249 
upper bound 
0.5869 
3.0000 
0.2524 
0.7332 
0.2440 
2.0705 
1.1380 
0.02 
0.4744 
2.995 
0.2464 
0.7529 
0.2207 
1.8319 
1.0810 
0.01 
0.4642 
2.995 
0.2493 
0.7559 
0.2315 
1.8522 
1.1012 
0.00875 
0.4629 
2.995 
0.2494 
0.7565 
0.2324 
1.8622 
1.1047 
0.0075 
0.4646 
2.995 
0.2495 
0.7574 
0.2333 
1.8739 
1.1111 
0.00625 
0.4616 
2.995 
0.2496 
0.7574 
0.2341 
1.8828 
1.1186 
5.5.7. Bibliography

[Hysing] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska, Quantitative benchmark computations of twodimensional bubble dynamics, International Journal for Numerical Methods in Fluids, 2009.

[Chabannes] V. Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles. PhD thesis, Université de Grenoble, 2013.

[Doyeux] V. Doyeux, Modélisation et simulation de systèmes multifluides, Application aux écoulements sanguins, PhD thesis, Université de Grenoble, 2014.
5.6. 3D Drop benchmark
The previous section described the strategy we used to track the interface. We couple it now to the Navier Stokes equation solver described in \cite{chabannes11:_high}. In the current section, we present a 3D extension of the 2D benchmark introduced in \cite{Hysing2009} and realised using Feel++ in \cite{Doyeux2012}.
5.6.1. Benchmark problem
The benchmark objective is to simulate the rise of a 3D bubble in a Newtonian fluid. The equations solved are the incompressible Navier Stokes equations for the fluid and the advection for the level set:
where \(\rho\) is the density of the fluid, \(\nu\) its viscosity, and \(\mathbf{g} \approx (0, 0.98)^T\) is the gravity acceleration.
The computational domain is \(\Omega \times \rbrack0, T\rbrack \) where \(\Omega\) is a cylinder which has a radius \(R\) and a heigth \(H\) so that \(R=0.5\) and \(H=2\) and \(T=3\). We denote \(\Omega_1\) the domain outside the bubble \( \Omega_1= \{\mathbf{x}  \phi(\mathbf{x})>0 \} \), \(\Omega_2\) the domain inside the bubble \( \Omega_2 = \{\mathbf{x}  \phi(\mathbf{x})<0 \} stem:[ and stem:[\Gamma\) the interface \( \Gamma = \{\mathbf{x}  \phi(\mathbf{x})=0 \} \). On the lateral walls and on the bottom walls, noslip boundary conditions are imposed, i.e. \(\mathbf{u} = 0\) and \(\mathbf{t} \cdot (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) \cdot \mathbf{n}=0\) where \(\mathbf{n}\) is the unit normal to the interface and \(\mathbf{t}\) the unit tangent. Neumann condition is imposed on the top wall i.e. \(\dfrac{\partial \mathbf{u}}{\partial \mathbf{n}}=\mathbf{0}\). The initial bubble is circular with a radius \(r_0 = 0.25\) and centered on the point \((0.5, 0.5, 0.)\). A surface tension force \(\mathbf{f}_{st}\) is applied on \(\Gamma\), it reads : \(\mathbf{f}_{st} = \int_{\Gamma} \sigma \kappa \mathbf{n} \simeq \int_{\Omega} \sigma \kappa \mathbf{n} \delta_{\varepsilon}(\phi)\) where \(\sigma\) stands for the surface tension between the twofluids and \(\kappa = \nabla \cdot (\frac{\nabla \mathbf{\phi}}{\nabla \phi})\) is the curvature of the interface. Note that the normal vector \(\mathbf{n}\) is defined here as \(\mathbf{n}=\frac{\nabla \phi}{\nabla \phi}\).
We denote with indices \(1\) and \(2\) the quantities relative to the fluid in respectively in \(\Omega_1\) and \(\Omega_2\). The parameters of the benchmark are \(\rho_1\), \(\rho_2\), \(\nu_1\), \(\nu_2\) and \(\sigma\) and we define two dimensionless numbers: first, the Reynolds number which is the ratio between inertial and viscous terms and is defined as : \(Re = \dfrac{\rho_1 \sqrt{\mathbf{g} (2r_0)^3}}{\nu_1}\); second, the E\"otv\"os number which represents the ratio between the gravity force and the surface tension \(E_0 = \dfrac{4 \rho_1 \mathbf{g} r_0^2}{\sigma}\). The table below reports the values of the parameters used for two different test cases proposed in~\cite{Hysing2009}.
Tests 
\(\rho_1\) 
\(\rho_2\) 
\(\nu_1\) 
\(\nu_2\) 
\(\sigma\) 
Re 
\(E_0\) 
Test 1 (ellipsoidal bubble) 
1000 
100 
10 
1 
24.5 
35 
10 
Test 2 (skirted bubble) 
1000 
1 
10 
0.1 
1.96 
35 
125 
The quantities measured in \cite{Hysing2009} are \(\mathbf{X_c}\) the center of mass of the bubble, \(\mathbf{U_c}\) its velocity and the circularity. For the 3D case we extend the circularity to the sphericity defined as the ratio between the surface of a sphere which has the same volume and the surface of the bubble which reads \(\Psi(t) = \dfrac{4\pi\left(\dfrac{3}{4\pi} \int_{\Omega_2} 1 \right)^{\frac{2}{3}}}{\int_{\Gamma} 1}\).
5.6.2. Simulations parameters
The simulations have been performed on the supercomputer SUPERMUC using 160 or 320 processors. The number of processors was chosen depending on the memory needed for the simulations. The table Numerical parameters used for the test 1 simulations summarize for the test 1 the different simulation properties and the table Mesh caracteristics give the carachteristics of each mesh.
h 
Number of processors 
\(\Delta t\) 
Time per iteration (s) 
Total Time (h) 
0.025 
360 
0.0125 
18.7 
1.25 
0.02 
360 
0.01 
36.1 
3.0 
0.0175 
180 
0.00875 
93.5 
8.9 
0.015 
180 
0.0075 
163.1 
18.4 
0.0125 
180 
0.00625 
339.7 
45.3 
h 
Tetrahedra 
Points 
Order 1 
Order 2 
0.025 
73010 
14846 
67770 
1522578 
0.02 
121919 
23291 
128969 
2928813 
0.0175 
154646 
30338 
187526 
4468382 
0.015 
217344 
41353 
292548 
6714918 
0.0125 
333527 
59597 
494484 
11416557 
The NavierStokes equations are linearized using the Newton’s method and we used a KSP method to solve the linear system. We use an Additive Schwarz Method for the preconditioning (GASM) and a LU method as a sub preconditionner. We run the simulations looking for solutions in finite element spaces spanned by Lagrange polynomials of order \((2,1,1)\) for respectively the velocity, the pressure and the level set.
5.6.3. Results Test 1: Ellipsoidal bubble
Accordind to the 2D results we expect that the drop would became ellipsoid. The figure~\ref{subfig:elli_sh} shows the shape of the drop at the final time step. The contour is quite similar to the one we obtained in the two dimensions simulations. The shapes are similar and seems to converge when the mesh size is decreasing. The drop reaches a stationary circularity and its topology does not change. The velocity increases until it attains a constant value. Figure~\ref{subfig:elli_uc} shows the results we obtained with different mesh sizes.