Modeling with Feel++
Advection Models
JSON file
Let’s describe how the json file is build.
Properties
Model  Name 

Advection 

Stokes 

Boundary Conditions
Boundary conditions are set in the json files
"BoundaryConditions":
{
"field":
{
"<condition type>":
{
"marker":
{
"expr": "value"
}
}
}
}
where all component of vectorial field can be access with <field>_x
for the x component where <field>
is the name of the vectorial field.
The condition types are explain in the two following tables.
Name  Options  Type 

Dirichlet 
faces, edges and componentwise 
"Dirichlet" 
Neumann 
scalar, vectorial 
"Neumann_scalar" or "Neumann_vectorial" 
Action
Here is a example code, that use this model on a ring domain.
Feel++ code
Here is the code
First at all, we define our model type with
typedef FeelModels::Advection< Simplex<FEELPP_DIM,1>, Lagrange<OrderAdvection, Scalar,Continuous,PointSetFekete> > model_type;
This definition allows us to create our advection model object FM like this
auto Adv = model_type::New("advection");
The method New
retrieve all data from the configuration and json files, as well build a mesh if need.
With this object, we can initialize our model parameters, such as boundaries conditions. Data on our model and on the numeric solver are then save and print on the terminal. This is made by
Adv>init(); Adv>printAndSaveInfo();
Now that our model is completed, we can solve the associated problem. To begin the resolution
Adv>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.
Adv>solve(); Adv>exportResults();
If it’s not, our model is time reliant, and a loop on time is necessary. Our model is then solved and the results are export at each time step.
for ( ; !Adv>timeStepBase()>isFinished(); Adv>updateTimeStep() ) { Adv>solve(); Adv>exportResults(); }
Code
{% include "../Examples/advection_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
# advection [advection] geofile=$cfgdir/ring2d.geo filename=$cfgdir/ring2d.json [exporter] directory=simul/advect2d/ring/data format=ensightgold
It can also be resolution dependent parameters such as mesh elements size, methods used to define our problem and solvers.
# time [advection.bdf] order=2 [ts] timeinitial=0.0 timestep=1 timefinal=1 steady=true [advection.gmsh] hsize=0.03 # backend advection and projection pcfactormatsolverpackagetype=mumps pctype=lu #kspmonitor=1 kspconvergedreason=true kspmaxit=1000 #snesmonitor=1 snesconvergedreason=true snesmaxitreuse=3 sneskspmaxitreuse=20
In this case, we choose LU as the preconditioner method, with a mesh size equal to 0.03. As for time discretization, we use a BDF at order 2.
Code
{% include "../Examples/ring2d.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": "Ring2d",
"ShortName": "Ring2d",
"Model": "Advection",
In this case, we have only boundary conditions to define. Here, we impose homogeneous Dirichlet conditions.
"BoundaryConditions":
{
"advection":
{
"Dirichlet":
{
"Bottom":
{
"expr":"0"
},
"Left":
{
"expr":"0"
},
"InnerCircle":
{
"expr":"0"
},
"OuterCircle":
{
"expr":"0"
}
}
}
}
Code
{% include "../Examples/ring2d.json" %}
Compilation/Execution
Once you’ve a build dir, you just have to realise the command make
at
{buildir}/applications/models/advection
This will generate an executable named feelpp_application_advection_2d
. To execute it, you need to give the path of the cfg file associated to your case, with configfile
.
For example
./feelpp_application_advection_2d configfile={sourcedir}/applications/models/advection/ring/ring2d.cfg
is how to execute the case ahead.
The result files are then stored by default in
feel/simul/advect2d/{domain_shape}/data/{processor_used}
If we return once again at our example, the result files are in
feel/simul/advect2d/ring/data/np_1
CFD Toolbox
Feel++ Book Contributors https://github.com/feelpp/book.feelpp.org/graphs/contributors
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"
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"
}
}
2.1. Generalised Newtonian fluid
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. Boundary Conditions
We start by a listing of boundary conditions supported in fluid mechanics model.
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.2. Dirichlet on pressure
p & = g \\ \boldsymbol{u}_f \times {\boldsymbol{ n }} & = \boldsymbol{0}
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.4. Slip
\boldsymbol{u}_f \cdot \boldsymbol{ n } = 0
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.6. Outlet flow
Option  Value  Default value  Description 

model 
free,windkessel 
free 
select an outlet modeling 
3.6.1. Free outflow
\boldsymbol{\sigma}_{f} \boldsymbol{ n } = \boldsymbol{0}
3.6.2. 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.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.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 
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.
5. Post Processing
"PostProcess":
{
"Fields":["field1","field2",...],
"Measures":
{
"<measure type>":
{
"label":
{
"<range type>":"value",
"fields":["field1","field3"]
}
}
}
}
5.1. 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
5.2. 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"
}
}
5.3. 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"],
}
6. Stabilization methods
6.1. GLS family
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 
6.2. CIP family
Documentation pending 
7. Run simulation
programme avalaible :

feelpp_toolbox_fluid_2d

feelpp_toolbox_fluid_3d
mpirun np 4 feelpp_toolbox_fluid_2d configfile=<myfile.cfg>
Equations
The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows
8. Linear elasticity
9. Hyperelasticity
9.1. SaintVenantKirchhoff
9.2. NeoHookean
9.2.1. Isochoric part : \(\boldsymbol{\Sigma}_s^\text{iso}\)
Name  \(\mathcal{W}_S(J_s)\)  \(\boldsymbol{\Sigma}_s^{\text{iso}}\) 

NeoHookean 
\(\mu_s J^{2/3}(\boldsymbol{I}  \frac{1}{3} \text{tr}(\boldsymbol{C}) \ \boldsymbol{C}^{1}) \) 
9.2.2. Volumetric part : \(\boldsymbol{\Sigma}_s^\text{vol}\)
Name  \(\mathcal{W}_S(J_s)\)  \(\boldsymbol{\Sigma}_s^\text{vol}\) 

classic 
\(\frac{\kappa}{2} \left( J_s  1 \right)^2\) 

simo1985 
\(\frac{\kappa}{2} \left( ln(J_s) \right)\) 
10. Axisymmetric reduced model
We interest us here to a 1D reduced model, named generalized string.
The axisymmetric form, which will interest us here, is a tube of length \(L\) and radius \(R_0\). It is oriented following the axis \(z\) and \(r\) represent the radial axis. The reduced domain, named \(\Omega_s^*\) is represented by the dotted line. So the domain, where radial displacement \(\eta_s\) is calculated, is \(\Omega_s^*=\lbrack0,L\rbrack\).
We introduce then \(\Omega_s^{'*}\), where we also need to estimate a radial displacement as before. The unique variance is this displacement direction.
The mathematical problem associated to this reduce model can be describe as
where \(\eta_s\), the radial displacement that satisfy this equation, \(k\) is the Timoshenko’s correction factor, and \(\gamma_v\) is a viscoelasticity parameter. The material is defined by its density \(\rho_s^*\), its Young’s modulus \(E_s\), its Poisson’s ratio \(\nu_s\) and its shear modulus \(G_s\)
At the end, we take \( \eta_s=0\text{ on }\partial\Omega_s^*\) as a boundary condition, which will fixed the wall to its extremities.
CSM Toolbox
11. Models
The solid mechanics model can be selected in json file :
"Model":"HyperElasticity"
Model  Name in json 

Linear Elasticity 

Hyper Elasticity 

When materials are closed to incompressibility formulation in displacement/pressure are available.
Model  Name  Volumic law 

SaintVenantKirchhoff 

classic, simo1985 
NeoHookean 

classic, simo1985 
option: mechanicalproperties.compressible.volumic_law
12. Materials
The Lamé coefficients are deducing from the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material we work on and can be express
"Materials":
{
"<marker>":
{
"name":"solid",
"E":"1.4e6",
"nu":"0.4",
"rho":"1e3"
}
}
where E
stands for the Young’s modulus in Pa, nu
the Poisson’s ratio (
dimensionless ) and rho
the density in \(kg\cdot m^{3}\).
13. Boundary Conditions
Name  Options  Type 

Dirichlet 
faces, edges and componentwise 
"Dirichlet" 
Neumann 
scalar, vectorial 
"Neumann_scalar" or "Neumann_vectorial" 
Pressure follower , 
Nonlinear boundary condition set in deformed domain 
TODO 
Robin 
TODO 
TODO 
14. Body forces
Name  Options  Type 

Expression 
Vectorial 
"VolumicForces" 
15. Post Process
15.1. Exports for visualisation
The fields allowed to be exported in the Fields
section are:

displacement

velocity

acceleration

stress or normalstress

pressure

materialproperties

pid

fsi

VonMises

Tresca

principalstresses

all
15.2. Measures

Points

Maximum

Minimum

VolumeVariation
15.2.1. Points
Same syntax as FluidMechanics with available Fields :

displacement

velocity

acceleration

pressure

principalstress0

principalstress1

principalstress2

sigma_xx, sigma_xy, …
15.2.2. Maximum/Minimum
The Maximum and minimum can be evaluated and save on .csv file. User need to define (i) <Type> ("Maximum" or "Minimum"), (ii) "<tag>" representing this data in the .csv file, (iii) "<marker>" representing the name of marked entities and (iv) the field where extremum is computed.
"<Type>":
{
"<tag>":
{
"markers":"marker>",
"fields":["displacement","velocity"]
}
}
15.2.3. VolumeVariation
"VolumeVariation":<marker>
16. Run simulations
programme avalaible :

feelpp_toolbox_solid_2d

feelpp_toolbox_solid_3d
mpirun np 4 feelpp_toolbox_solid_2d configfile=<myfile.cfg>
Fluid Structure Interaction Models
18. Fluid structure coupling conditions
In order to have a correct fluidstructure model, we need to add to the solid model and the fluid model equations some coupling conditions :
\(\boldsymbol{(1)}, \boldsymbol{(2)}, \boldsymbol{(3)}\) are the fluidstruture coupling conditions, respectively velocities continuity, constraint continuity and geometric continuity.
18.1. Fluid structure coupling conditions with 1D reduced model
For the coupling conditions, between the 2D fluid and 1D structure, we need to modify the original ones \(\boldsymbol{(1)},\boldsymbol{(2)}, \boldsymbol{(3)}\) by
18.2. Variables, symbols and units
Notation 
Quantity 
Unit 
\(\boldsymbol{u}_f\) 
fluid velocity 
\(m.s^{1}\) 
\(\boldsymbol{\sigma}_f\) 
fluid stress tensor 
\(N.m^{2}\) 
\(\boldsymbol{\eta}_s\) 
displacement 
\(m\) 
\(\boldsymbol{F}_s\) 
deformation gradient 
dimensionless 
\(\boldsymbol{\Sigma}_s\) 
second PiolaKirchhoff tensor 
\(N.m^{2}\) 
\(\mathcal{A}_f^t\) 
Arbitrary Lagrangian Eulerian ( ALE ) map 
dimensionless 
and
ThermoElectric Toolbox
This section focuses on defining and solving of a thermoelectric model, which couples heat transfer and electrostatic models. The current flow inside an electrical conductor leads to the heating of the material, due to Joule effect. The temperature in the conductor can be deduced from the standard heat equation whose source term reads from the electrical potential to mimic the Joule effect.
In the following sections, the domain is denoted by \(\Omega \subset \mathbb{R}^d\), \(d=2,3\) and its border by \(\Gamma\).
19. Electrostatic problem
The electrostatic model consists in a diffusion equation. The electrical potential \(V\) satisfies
where \(\sigma\) is the electrical conductivity.
The electrical conductivity \(\sigma\) often depends on the temperature. The thermoelectric model becomes nonlinear in this case. 
The current flow within the conductor is modeled by a difference in electrical potential between the current input \(\Gamma_{in} \subset \Gamma\) and the current output \(\Gamma_{out} \subset \Gamma\).

\(V=0 ~~\text{on}~~ \Gamma_{in}\)

\(V=V_{out} ~~\text{on}~~ \Gamma_{out}\).
We consider in this model that all other boundaries are electrically insulating. This conditions reads from the current density \(\mathbf{j} =  \sigma \nabla V\) as

\(\mathbf{j}\cdot\mathbf{n} = \sigma\nabla V\cdot\mathbf{n} = 0 ~~\text{on}~~ \Gamma \setminus (\Gamma_{in} \cup \Gamma_{out})\).
The local form of Joule’s law defines the heat source corresponding to Joule effect as the dot product of the current density \(\mathbf{j} =  \sigma \nabla V\) by the electric field \(E=\nabla V\). This product \(P = \sigma \nabla V \cdot \nabla V\) will be the source term of the heat equation defined in the next section.
20. Heat transfer problem
The temperature \(T\) satisfies the standard heat equation with the previously introduced Joule effect as source term,
where \(k\) is the thermal conductivity.
As for the electrical conductivity \(\sigma\), the thermal conductivity \(k\) can also depends on temperature. 
The thermal flux can eventually be controled by a water cooling process which limit the heating. The thermal exchange between the electrical conductor and the cooling water is governed by a convection phenomenon involving a heat transfert coefficient \(h\) and the water temperature \(T_c\). This cooling process is handled by a Robin condition on the cooled regions \(\Gamma_{c} \subset \Gamma\)

\(k\nabla T\cdot\mathbf{n} = h(TT_c) ~~\text{on}~~ \Gamma_{C}\)
We consider that the thermal exchanges are limited to the cooled region, which amounts to apply an homogeneous Neumann condition on all other boundaries.

\(k\nabla T\cdot \mathbf{n} = 0 ~~\text{on}~~ \Gamma \setminus \Gamma_{C}\)
21. ThermoElectric problem
In the thermoelectric model, we have 4 parameters

\(\sigma\) the electrical conductivity

\(k\) the thermal conductivity

\(T_c\) the water cooling temperature

\(h\) the heat transfer coefficient
These parameters can be scalars of fields. 
We gather in the following table parameter ranges, nominal values as well as units for \(\sigma, k, T_c, h\).
Parameters 
Ranges 
Nominal value 
Units 
\(\sigma\) 
\([52.10^{6};58.10^{6}\)] 
\(53\cdot 10^{6}\) 
\(S \cdot m^{1}\) 
\(k\) 
\([360;380\)] 
\(370\) 
\(W\cdot m^{1} \cdot K^{1}$\) 
\(T_c\) 
\([293;310\)] 
\(300\) 
K 
\(h\) 
\([70000;90000\)] 
\(850000\) 
\(W \cdot m^{2} \cdot K^{1}\) 
In the linear case, we first solve for \(V\) and then for \(T\) using \(V\) to compute the Joule effect that generates heat inside \(\Omega\).
22. Running the Thermoelectric application
Using Docker, you can run Feel++ model application and in particular the thermoelectric model using the following command
$ docker run it v $HOME/feel:/feel feelpp/feelpptoolboxes:latest
Then type the following command in docker environment to run the model
$ cd Testcases/models/thermoelectric/test
$ mpirun np 4 /usr/local/bin/feelpp_toolbox_thermoelectric_3d configfile model.cfg
Unresolved directive in README.adoc  include::Electromagnetism/README.adoc[]