CFD Toolbox

1.1. Models

The CFD Toolbox supports both the Stokes and the incompressible Navier-Stokes equations.

The fluid mechanics model (Navier-Stokes or Stokes) can be selected in json file:

Listing : select fluid model
"Model": "Navier-Stokes"

1.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.

Table 1. Correspondance between fluid parameters and symbols in JSon files
Parameter Symbol

\(\mu_f\)

mu

\(\rho_f\)

rho

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.

Listing : Materials section
"Materials":
{
    "<marker>"
    {
        "name":"water",
        "rho":"1.0e3",
        "mu":"1.0"
    }
}

1.2.1. Generalised Newtonian fluid

The non Newtonian properties are defined in cfg file in fluid section.

The viscosity law is activated by:

Table 2. Viscosity law
option values

viscosity.law

newtonian, power_law, walburn-schneck_law, carreau_law, carreau-yasuda_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

walburn-schneck_law

hematocrit

TPMA

walburn-schneck_law.C1

walburn-schneck_law.C2

walburn-schneck_law.C3

walburn-schneck_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

carreau-yasuda_law

viscosity.zero_shear

viscosity.infinite_shear

carreau-yasuda_law.lambda

carreau-yasuda_law.n

carreau-yasuda_law.a

\(kg/(m \times s)\)

\(kg/(m \times s)\)

dimensionless

dimensionless

dimensionless

1.3. Boundary Conditions

We start by a listing of boundary conditions supported in fluid mechanics model.

1.3.1. Dirichlet on velocity

A Dirichlet condition on velocity field reads:

Dirichlet condition

\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

1.3.2. Dirichlet on pressure

p & = g \\ \boldsymbol{u}_f \times {\boldsymbol{ n }} & = \boldsymbol{0}

1.3.3. Neumann

Table 3. Neumann options
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}\)

1.3.4. Slip

\boldsymbol{u}_f \cdot \boldsymbol{ n } = 0

1.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\)

Table 4. Inlet flow options
Option Value Default value Description

shape

constant,parabolic

select a shape profile for inflow

constraint

velocity_max,flow_rate

give a constraint wich controle velocity

expr

string

symbolic expression of constraint value

1.3.6. Outlet flow

Table 5. Outlet flow options
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 3-element 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}

Table 6. Windkessel options
Option Value Description

windkessel_coupling

explicit, implicit

coupling type with the Navier-Stokes equation

windkessel_Rd

real

distal resistance

windkessel_Rp

real

proximal resistance

windkessel_Cd

real

capacitance

1.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.

Listing : boundary conditions in json
"BoundaryConditions":
{
    "<field>":
    {
        "<bc_type>":
        {
            "<marker>":
            {
                "<option1>":"<value1>",
                "<option2>":"<value2>",
                // ...
            }
        }
    }
}

1.3.8. Options summary

Table 7. Boundary conditions
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

1.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.

1.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 normal-stress

  • wall-shear-stress

  • 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"
    }
}
Forces

compute lift and drag

"Forces":["fsi-wall","fluid-cylinder"]
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"],
}

1.6. Stabilization methods

1.6.1. GLS family

Galerkin leat-Square (GLS) stabilization methods can be activated from the cfg file by adding stabilization-gls=1 in the fluid prefix. Others options available are enumerated in the next table and must be given with the prefix fluid.stabilization-gls.

Table 8. GLS stabilzation options in CFG
Option Value Default value Description

type

gls,supg,gls-no-pspg, supg-pspg, pspg

gls

type of stabilization

parameter.method

eigenvalue,doubly-asymptotic-approximation

eigenvalue

method used to compute the stabilization parameter

parameter.hsize.method

hmin,h,meas

hmin

method used for evalute an element mesh size

parameter.eigenvalue.penal-lambdaK

real

0.

add a mass matrix scaled by this factor in the eigen value problem for the stabilization parameter

convection-diffusion.location.expressions

string

if given, the stabilization is apply only on mesh location which verify expr>0

1.6.2. CIP family

Documentation pending

1.7. Run simulation

programme avalaible :

  • feelpp_toolbox_fluid_2d

  • feelpp_toolbox_fluid_3d

mpirun -np 4 feelpp_toolbox_fluid_2d --config-file=<myfile.cfg>

1. Equations

The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows

\[ \rho^*_{s} \frac{\partial^2 \boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot \left(\boldsymbol{F}_s \boldsymbol{\Sigma}_s\right) = \boldsymbol{f}^t_s\]

1.1. Linear elasticity

\[\begin{align} \boldsymbol{F}_s &= \text{Identity} \\ \boldsymbol{\Sigma}_s &=\lambda_s tr( \boldsymbol{\epsilon}_s)\boldsymbol{I} + 2\mu_s\boldsymbol{\epsilon}_s \end{align}\]

1.2. Hyper-elasticity

1.2.1. Saint-Venant-Kirchhoff

\[\boldsymbol{\Sigma}_s=\lambda_s tr( \boldsymbol{E}_s)\boldsymbol{I} + 2\mu_s\boldsymbol{E}_s\]

1.2.2. Neo-Hookean

\[\boldsymbol{\Sigma}_s= \mu_s J^{-2/3}(\boldsymbol{I} - \frac{1}{3} \text{tr}(\boldsymbol{C}) \ \boldsymbol{C}^{-1})\]
\[\boldsymbol{\Sigma}_s^ = \boldsymbol{\Sigma}_s^\text{iso} + \boldsymbol{\Sigma}_s^\text{vol}\]
Isochoric part : \(\boldsymbol{\Sigma}_s^\text{iso}\)
Table 9. Isochoric law
Name \(\mathcal{W}_S(J_s)\) \(\boldsymbol{\Sigma}_s^{\text{iso}}\)

Neo-Hookean

\(\mu_s J^{-2/3}(\boldsymbol{I} - \frac{1}{3} \text{tr}(\boldsymbol{C}) \ \boldsymbol{C}^{-1}) \)

Volumetric part : \(\boldsymbol{\Sigma}_s^\text{vol}\)
Table 10. Volumetric law
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)\)

1.3. 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.

Reduced Model Geometry
Figure 1 : Geometry of the reduce model

The mathematical problem associated to this reduce model can be describe as

\[ \rho^*_s h \frac{\partial^2 \eta_s}{\partial t^2} - k G_s h \frac{\partial^2 \eta_s}{\partial x^2} + \frac{E_s h}{1-\nu_s^2} \frac{\eta_s}{R_0^2} - \gamma_v \frac{\partial^3 \eta}{\partial x^2 \partial t} = f_s.\]

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.

2. CSM Toolbox

2.1. Models

The solid mechanics model can be selected in json file :

Listing : select solid model
"Model":"Hyper-Elasticity"
Table 11. Table of Models for model option
Model Name in json

Linear Elasticity

Elasticity

Hyper Elasticity

Hyper-Elasticity

When materials are closed to incompressibility formulation in displacement/pressure are available.

Table 12. Table of Models for material_law with hyper elasticity model
Model Name Volumic law

Saint-Venant-Kirchhoff

SaintVenantKirchhoff

classic, simo1985

NeoHookean

NeoHookean

classic, simo1985

option: mechanicalproperties.compressible.volumic_law

2.2. 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

\[\lambda_s = \frac{E_s\nu_s}{(1+\nu_s)(1-2\nu_s)} \hspace{0.5 cm} , \hspace{0.5 cm} \mu_s = \frac{E_s}{2(1+\nu_s)}\]
Materials section
"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}\).

2.3. Boundary Conditions

Table 13. Boundary conditions
Name Options Type

Dirichlet

faces, edges and component-wise

"Dirichlet"

Neumann

scalar, vectorial

"Neumann_scalar" or "Neumann_vectorial"

Pressure follower ,

Nonlinear boundary condition set in deformed domain

TODO

Robin

TODO

TODO

2.4. Body forces

Table 14. Volumic forces
Name Options Type

Expression

Vectorial

"VolumicForces"

2.5. Post Process

2.5.1. Exports for visualisation

The fields allowed to be exported in the Fields section are:

  • displacement

  • velocity

  • acceleration

  • stress or normal-stress

  • pressure

  • material-properties

  • pid

  • fsi

  • Von-Mises

  • Tresca

  • principal-stresses

  • all

2.5.2. Measures

  • Points

  • Maximum

  • Minimum

  • VolumeVariation

Points

Same syntax as FluidMechanics with available Fields :

  • displacement

  • velocity

  • acceleration

  • pressure

  • principal-stress-0

  • principal-stress-1

  • principal-stress-2

  • sigma_xx, sigma_xy, …​

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"]
    }
}
VolumeVariation
"VolumeVariation":<marker>

2.6. Run simulations

programme avalaible :

  • feelpp_toolbox_solid_2d

  • feelpp_toolbox_solid_3d

mpirun -np 4 feelpp_toolbox_solid_2d --config-file=<myfile.cfg>

3. Fluid Structure Interaction Models

The Fluid Structure models are formed from the combination of a Solid model and a Fluid model.

3.2. Fluid structure coupling conditions

In order to have a correct fluid-structure model, we need to add to the solid model and the fluid model equations some coupling conditions :

\[ \frac{\partial \boldsymbol{\eta_{s}} }{\partial t} - \boldsymbol{u}_f \circ \mathcal{A}_{f}^t = \boldsymbol{0} , \quad \text{ on } \Gamma_{fsi}^* \times \left[t_i,t_f \right] \quad \boldsymbol{(1)}\]
\[ \boldsymbol{F}_{s} \boldsymbol{\Sigma}_{s} \boldsymbol{n}^*_s + J_{\mathcal{A}_{f}^t} \hat{\boldsymbol{\sigma}}_f \boldsymbol{F}_{\mathcal{A}_{f}^t}^{-T} \boldsymbol{n}^*_f = \boldsymbol{0} , \quad \text{ on } \Gamma_{fsi}^* \times \left[t_i,t_f \right] \quad \boldsymbol{(2)}\]
\[ \boldsymbol{\varphi}_s^t - \mathcal{A}_{f}^t = \boldsymbol{0} , \quad \text{ on } \Gamma_{fsi}^* \times \left[t_i,t_f \right] \quad \boldsymbol{(3)}\]

\(\boldsymbol{(1)}, \boldsymbol{(2)}, \boldsymbol{(3)}\) are the fluid-struture coupling conditions, respectively velocities continuity, constraint continuity and geometric continuity.

3.2.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

\[\dot{\eta}_s \boldsymbol{e}_r - \boldsymbol{u}_f = \boldsymbol{0} \quad \boldsymbol{(1.2)}\]
\[f_s + \left(J_{\mathcal{A}_f^t} \boldsymbol{F}_{\mathcal{A}_f^t}^{-T} \hat{\boldsymbol{\sigma}}_f \boldsymbol{n}^*_f\right) \cdot \boldsymbol{e}_r = 0 \quad \boldsymbol{(2.2)}\]
\[\boldsymbol{\varphi}_s^t - \mathcal{A}_f^t = \boldsymbol{0} \quad \boldsymbol{(3.2)}\]

3.2.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 Piola-Kirchhoff tensor

\(N.m^{-2}\)

\(\mathcal{A}_f^t\)

Arbitrary Lagrangian Eulerian ( ALE ) map

dimensionless

and

\[\boldsymbol{F}_{\mathcal{A}_f^t} = \boldsymbol{\mathrm{x}}^* + \nabla \mathcal{A}_f^t\]
\[J_{\mathcal{A}_f^t} = det(\boldsymbol{F}_{\mathcal{A}_f^t})\]

4. Advection Models

Table of Contents

4.1. JSON file

Let’s describe how the json file is build.

4.1.1. Properties

Table 15. Table of Models for model option
Model Name

Advection

Advection

Stokes

Stokes

4.1.2. Boundary Conditions

Boundary conditions are set in the json files

Example for boundary conditions
"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.

Table 16. Boundary conditions
Name Options Type

Dirichlet

faces, edges and component-wise

"Dirichlet"

Neumann

scalar, vectorial

"Neumann_scalar" or "Neumann_vectorial"

4.2. Action

Here is a example code, that use this model on a ring domain.

4.2.1. 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" %}

4.2.2. 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]
time-initial=0.0
time-step=1
time-final=1
steady=true

[advection.gmsh]
hsize=0.03

# backend advection and projection
pc-factor-mat-solver-package-type=mumps
pc-type=lu

#ksp-monitor=1
ksp-converged-reason=true
ksp-maxit=1000
#snes-monitor=1
snes-converged-reason=true
snes-maxit-reuse=3
snes-ksp-maxit-reuse=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" %}

4.2.3. 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" %}

4.2.4. 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 --config-file.

For example

./feelpp_application_advection_2d --config-file={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

5. Thermo-Electric Toolbox

Table of Contents

This section focuses on defining and solving of a thermo-electric 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\).

5.1. Electrostatic problem

The electrostatic model consists in a diffusion equation. The electrical potential \(V\) satisfies

\[\nabla\cdot(-\sigma\nabla V) = 0 ~~\text{in}~~ \Omega\]

where \(\sigma\) is the electrical conductivity.

The electrical conductivity \(\sigma\) often depends on the temperature. The thermo-electric model becomes non-linear 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.

5.2. Heat transfer problem

The temperature \(T\) satisfies the standard heat equation with the previously introduced Joule effect as source term,

\[\nabla\cdot(-k\nabla T) = \sigma\nabla V\cdot \nabla V ~~\text{in}~~ \Omega\]

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(T-T_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}\)

5.3. Thermo-Electric problem

In the thermo-electric 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\).

Table 17. Table of parameters for the simulation

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\).

5.4. Running the Thermo-electric application

Using Docker, you can run Feel++ model application and in particular the thermo-electric model using the following command

Docker command
$ docker run -it -v $HOME/feel:/feel feelpp/feelpp-toolboxes:latest

Then type the following command in docker environment to run the model

Running Thermo-Electric model
$ cd Testcases/models/thermoelectric/test
$ mpirun -np 4 /usr/local/bin/feelpp_toolbox_thermoelectric_3d --config-file model.cfg

6. Learning Feel++

:leveloffset:+1 = Generic Partial Differential Equations

:leveloffset:+1

7. The Laplacian

7.1. Problem statement

We are interested in this section in the conforming finite element approximation of the following problem:

Laplacian problem

Look for \(u\) such that

\[\left\{\begin{split} -\Delta u &= f \text{ in } \Omega\\ u &= g \text{ on } \partial \Omega_D\\ \frac{\partial u}{\partial n} &=h \text{ on } \partial \Omega_N\\ \frac{\partial u}{\partial n} + u &=l \text{ on } \partial \Omega_R \end{split}\right.\]
\(\partial \Omega_D\), \(\partial \Omega_N\) and \(\partial \Omega_R\) can be empty sets. In the case \(\partial \Omega_D =\partial \Omega_R = \emptyset\), then the solution is known up to a constant.

In the implementation presented later, \(\partial \Omega_D =\partial \Omega_N = \partial \Omega_R = \emptyset\), then we set Dirichlet boundary conditions all over the boundary. The problem then reads like a standard laplacian with inhomogeneous Dirichlet boundary conditions:

Laplacian Problem with inhomogeneous Dirichlet conditions

Look for \(u\) such that

Inhomogeneous Dirichlet Laplacian problem
\[-\Delta u = f\ \text{ in } \Omega,\quad u = g \text{ on } \partial \Omega\]

7.2. Variational formulation

We assume that \(f, h, l \in L^2(\Omega)\). The weak formulation of the problem then reads:

Laplacian problem variational formulation

Look for \(u \in H^1_{g,\Gamma_D}(\Omega)\) such that

Variational formulation
\[\displaystyle\int_\Omega \nabla u \cdot \nabla v +\int_{\Gamma_R} u v = \displaystyle \int_\Omega f\ v+ \int_{\Gamma_N} g\ v + \int_{\Gamma_R} l\ v,\quad \forall v \in H^1_{0,\Gamma_D}(\Omega)\]

7.3. Conforming Approximation

We now turn to the finite element approximation using Lagrange finite element. We assume \(\Omega\) to be a segment in 1D, a polygon in 2D or a polyhedron in 3D. We denote \(V_\delta \subset H^1(\Omega)\) an approximation space such that \(V_{g,\delta} \equiv P^k_{c,\delta}\cap H^1_{g,\Gamma_D}(\Omega)\).

The weak formulation reads:

Laplacian problem weak formulation

Look for \(u_\delta \in V_\delta \) such that

\[\displaystyle\int_{\Omega_\delta} \nabla u_{\delta} \cdot \nabla v_\delta +\int_{\Gamma_{R,\delta}} u_\delta\ v_\delta = \displaystyle \int_{\Omega_\delta} f\ v_\delta+ \int_{\Gamma_{N,\delta}} g\ v_\delta + \int_{\Gamma_{R,\delta}} l\ v_\delta,\quad \forall v_\delta \in V_{0,\delta}\]
from now on, we omit \(\delta\) to lighten the notations. Be careful that it appears both the geometrical and approximation level.

7.4. Feel++ Implementation

In Feel++, \(V_{g,\delta}\) is not built but rather \(P^k_{c,\delta}\).

The Dirichlet boundary conditions can be treated using different techniques and we use from now on the elimination technique.

We start with the mesh

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
the keyword auto enables type inference, for more details see Wikipedia C++11 page.

Next the discretization setting by first defining Vh=Pch<k>(mesh) \(\equiv P^k_{c,h}\), then elements of Vh and expressions f, n and g given by command line options or configuration file.

    auto Vh = Pch<2>( mesh );
    auto u = Vh->element("u");
    auto mu = doption(_name="mu");
    auto f = expr( soption(_name="functions.f"), "f" );
    auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
    auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
    auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
    auto g = expr( soption(_name="functions.g"), "g" );
    auto v = Vh->element( g, "g" );

at the following line

    auto v = Vh->element( g, "g" );

v is set to the expression g, which means more precisely that v is the interpolant of g in Vh.

the variational formulation is implemented below, we define the bilinear form a and linear form l and we set strongly the Dirichlet boundary conditions with the keyword on using elimination. If we don’t find Dirichlet, Neumann or Robin in the list of physical markers in the mesh data structure then we impose Dirichlet boundary conditions all over the boundary.

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=f*id(v));
    l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
    l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
    toc("l");

    tic();
    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=mu*gradt(u)*trans(grad(v)) );
    a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
    a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
    //! if no markers Robin Neumann or Dirichlet are present in the mesh then
    //! impose Dirichlet boundary conditions over the entire boundary
    if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
        a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
    toc("a");

    tic();
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);
    toc("a.solve");

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

    tic();
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

    toc("Exporter");
    return 0;

}

We have the following correspondance:

Element sets Domain

elements(mesh)

\(\Omega\)

boundaryfaces(mesh)

\(\partial \Omega\)

markedfaces(mesh,"Dirichlet")

\(\Gamma_D\)

markedfaces(mesh,"Neumann")

\(\Gamma_R\)

markedfaces(mesh,"Robin")

\(\Gamma_R\)

next we solve the algebraic problem

Listing: solve algebraic system
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);

next we compute the \(L^2\) norm of \(u_\delta-g\), it could serve as an \(L^2\) error if \(g\) was manufactured to be the exact solution of the Laplacian problem.

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

and finally we export the results, by default it is in the ensight gold format and the files can be read with Paraview and Ensight. We save both \(u\) and \(g\).

Listing: export Laplacian results
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

7.5. Testcases

The Feel++ Implementation comes with testcases in 2D and 3D.

7.5.1. circle

circle is a 2D testcase where \(\Omega\) is a disk whose boundary has been split such that \(\partial \Omega=\partial \Omega_D \cup \partial \Omega_N \cup \partial \Omega_R\).

Here are some results we can observe after use the following command

cd Testcases/quickstart/circle
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file circle.cfg

This give us some data such as solution of our problem or the mesh used in the application.

ucircle

meshCircle

Solution \(u_\delta\)

Mesh

7.5.2. feelpp2d and feelpp3d

This testcase solves the Laplacian problem in \(\Omega\) an quadrangle or hexadra containing the letters of Feel++

feelpp2d

After running the following command

cd Testcases/quickstart/feelpp2d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file feelpp2d.cfg

we obtain the result \(u_\delta\) and also the mesh

ufeelpp2d

/images/Laplacian/TestCases/Feelpp2d/meshfeelpp2d.png[]

Solution \(u_\delta\)

Mesh

feelpp3d

We can launch this application with the current line

cd Testcases/quickstart/feelpp3d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_3d --config-file feelpp3d.cfg

When it’s finish, we can extract some informations

ufeelpp3d

meshfeelpp3d

Solution \(u_\delta\)

Mesh

8. Levelset

Having the possibility to determine where two regions meeting can be really useful in some scientific domains. That’s why the levelset method is born.

8.1. Levelset introduction

8.1.1. Levelset function

By using a scalar function \phi, define on all regions as the null value is obtained when it’s placed on an interface of two domains.

We denote \Omega_1 and \Omega_2 two domains with \Gamma the interface betwen them. Then \phi can be define as

\phi(\boldsymbol{x}) = \left\{ \begin{array}{cccc} \text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in &\Omega_1 \\ 0, & \boldsymbol{x}& \in &\Gamma\\ -\text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in & \Omega_2 \end{array} \right.

with \text{dist}(\boldsymbol{x}, \Gamma ) = \underset{\boldsymbol{y} \; \in \; \Gamma}{\min}( |\boldsymbol{x} - \boldsymbol{y}| ).

This function \phi had also the following property |\nabla\phi|=1.

Moreover, the unit normal vector \boldsymbol{n} outgoing from the interface and the curvature \mathcal{\kappa} can be obtained from the levelset function.

\boldsymbol{n}=\frac{\nabla\phi}{|\nabla\phi} \\ \mathcal{\kappa}=\nabla \cdot \boldsymbol{n}= \nabla \cdot \frac{\nabla\phi}{|\nabla\phi|}

Now we have exposed the levelset function, we need to define how the levelset will evolve and will spread into all the space. To do this, we use the following advection equation :

\partial_t\phi+\boldsymbol{u}\cdot\nabla\phi=0

where \boldsymbol{u} is an incompressible velocity field.

8.1.2. Heaviside and Dirac functions

We define also the regularized Heaviside function H_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \dfrac{1}{2} \left(1+\dfrac{\phi}{\varepsilon}+\dfrac{\sin\left(\dfrac{\pi \phi}{\varepsilon}\right)}{\pi}\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 1, & \phi \geq \varepsilon \end{array} \right.

and the regularized Dirac function \delta_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \displaystyle\dfrac{1}{2 \varepsilon} \left(1+\cos\left(\dfrac{\pi \phi}{\varepsilon}\right)\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 0, & \phi \geq \varepsilon. \end{array} \right.

The first one gives a different value to each side of the interface ( here 0 in and 1 out ). The second one allow us to define quantities, with value different from 0 at the interface. A typical value of \varepsilon in literature is 1.5h where h is the mesh step size.

It should be noted that these functions allow us to determine respectively the volume and the surface of the interface by V^+_{\varepsilon} = \int_{\Omega} H_\varepsilon \\ S^{\Gamma}_{\varepsilon} = \int_{\Omega} \delta_\varepsilon

9. Solid rotation of a slotted disk

We describe the benchmark proposed by Zalesak.

Computer codes, used for the acquisition of results, are from Vincent Doyeux.

9.1. Problem Description

In order to test our interface propagation method, i.e. the levelset method \phi, we will study the rotation of a slotted disk into a square domain. The geometry can be represented as

Slotted Disk Geometry
Figure 1 : Initial Geometry.

We denote \Omega, the square domain [0,1]\times[0,1]. The center of the slotted disk is placed at (0.5,0.75).
To model the rotation, we will apply an angular velocity, centered in (0.5,0.5), as the disk is back to its initial position after t_f=628.

During this test, we observe three different errors to measure the quality of our method. With these values, two kinds of convergence will be studied : the time convergence, with different time step on an imposed grid and the space one, where the space discretization and the time step are linked by a relation. Several stabilization methods are used such as CIP ( Continuous Interior Penalty ) or SUPG ( Streamline-Upwing/Petrov-Galerkin ).

9.1.1. Boundary conditions

We set a Neumann boundary condition on the boundary of the domain.

9.1.2. Initial conditions

The velocity is imposed as \boldsymbol{u}=\left( \frac{\pi}{314} (50-y),\frac{\pi}{314} (x-50) \right)

Here is the velocity look in the square domain

Velocity Geometry
Figure 2 : Imposed velocity in \Omega.

9.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 18. Fixed and Variable Input Parameters

Name

Description

Nominal Value

Units

r

disk radius

0.15

m

l

slot base

0.05

m

h

slot height

0.25

m

t_f

slotted disk rotation period

628

s

9.3. Outputs

We observe during this benchmark three different errors.
First at all, the mass error, define by e_{\text{m}} = \frac{ \left| m_{\phi_f} - m_{\phi_0}\right| }{m_{\phi_0}} = \frac{ \left| \displaystyle \int_{\Omega} \chi( \phi_f < 0 ) - \displaystyle \int_{\Omega} \chi( \phi_0 < 0 ) \right| }{ \displaystyle \int_{\Omega} \chi( \phi_0 < 0 )}

where \chi is the characteristic function. \chi( f( \phi ) ) = \left\{ \begin{array}{rcl} 1 & \text{ if } & f( \phi ) \neq 0 \\ 0 & \text{ if } & f( \phi ) = 0 \end{array} \right.

However, mass is gain and loose at different emplacements on the mesh, and at the same time, with the level set method.

Secondly, the sign change error e_{\text{sc}} = \sqrt{ \int_\Omega \left( (1-H_0) - (1-H_f) \right)^2 }

with H_0=H_\epsilon(\phi_0) and H_f=H_\epsilon(\phi_f), H_\epsilon the smoothed Heaviside function of thickness 2ε.

This error is better to define the interface displacement. In fact, we can determine where \phi_0\phi_f<0, in other words where the interface has moved.

Finally, we define the classical L^2 error at the interface, as e_{L^2} = \sqrt{ \frac{1}{\displaystyle \int_\Omega \chi( \delta(\phi_0) > 0 ) } \int_\Omega (\phi_0 - \phi_f)^2 \chi( \delta(\phi_0) > 0 ) }.

9.4. Discretization

9.4.1. Time convergence

For this case, we set a fixed grid with mesh step size h=0.04, and so 72314 degree of freedom on a \mathbb{P}^1.

Then, after the disk made one round, we measure the errors obtained from two different discretizations ( BDF2 and Euler ) and compared them.

We repeat this with several time step dt\in \{2.14, 1, 0.5, 0.25, 0.20\}.

Only one stabilization method is used : SUPG

9.4.2. Space convergence

We define the following relation, between time step and mesh step size : dt=C\frac{h}{U_{max}}

where C<1 constant and U_{max} the maximum velocity of \Omega.

From the definition of our velocity, U_{max} is reached at the farthest point from the center of \Omega. In this case, we have U_{max}=0.007, and we set C=0.8.

We use the BDF2 method for time discretization. As in time convergence, we wait one round of the disk to measure the errors and we repeat this test for these values of h: 0.32, 0.16, 0.08, 0.04.

We compare the results from different stabilization methods : CIP, SUPG, GLS ( Galerkin-Least-Squares ) and SGS ( Sub-Grid Scale ).

9.6. Results

9.6.1. Time convergence

Table 19. Time convergence with Euler scheme

dt

e_{L^2}

e_{sc}

e_m

2.14

0.0348851

0.202025

0.202025

1.00

0.0187567

0.147635

0.147635

0.5

0.0098661

0.10847

0.10847

0.25

0.008791

0.0782569

0.0782569

0.20

0.00803373

0.0670677

0.0670677

Table 20. Time convergence with BDF2 scheme

dt

e_{L^2}

e_{sc}

e_m

2.14

0.0118025

0.0906617

0.0492775

1.00

0.00436957

0.0445275

0.0163494

0.5

0.00173637

0.0216359

0.0100621

0.25

0.001003

0.0125971

0.00354644

0.20

0.000949343

0.0117449

0.00317368

Time mass error
Figure 2 : Mass error of Zalesac benchmark
Time sign change
Figure 3 : Sign change error of Zalesac benchmark
Time L2 error
Figure 4 : L^2 error of Zalesac benchmark
Time shape
Figure 5 : Slotted disk shape after a round

9.6.2. Space convergence

stab

h

e_{L^2}

e_{sc}

e_m

CIP

0.32

0.0074

0.072

0.00029

0.16

0.0046

0.055

0.00202

0.08

0.0025

0.033

0.00049

0.04

0.0023

0.020

0.00110

SUPG

0.32

0.012

0.065

0.01632

0.16

0.008

0.049

0.07052

0.08

0.004

0.030

0.00073

0.04

0.001

0.018

0.00831

GLS

0.32

0.013

0.066

0.02499

0.16

0.008

0.051

0.05180

0.08

0.004

0.031

0.00805

0.04

0.001

0.019

0.00672

SGS

0.32

0.012

0.065

0.01103

0.16

0.008

0.050

0.07570

0.08

0.004

0.030

0.00084

0.04

0.001

0.018

0.00850

Space sign change
Figure 6 : Sign change error of Zalesac benchmark
Space L2 error
Figure 7 : L^2 error of Zalesac benchmark
Space Shape
Figure 8 : Slotted disk shape after a round

9.6.3. Conclusion

Let’s begin with time convergence results. Tables shows us that sign change error is better to define the quality of the chosen scheme than the mass error. In fact, the loss of mass somewhere can be nullify by a gain of mass elsewhere. Sign change error shows half an order gain from Euler scheme to BDF2 one, as L^2 errors show us a gain of one order. For the slotted disk shape, BDF2 uses the two previous iterations to obtain the current result, while Euler only need the previous iteration. This explain why we can see an asymmetrical tendency in the first one.

As for space convergence, the different stabilization methods we used give us the same convergence rate equals to 0.6, with close error values, for the sign change error. For the L^2 error case, it’s not as evident as the previous one. Aside the CIP stabilization method with a 0.6 convergence rate, the others show us a convergence rate of 0.9.

9.7. Bibliography

References for this benchmark
  • [Zalesak] Steven T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, Journal of Computational Physics, 1979.

  • [Doyeux] Vincent Doyeux, Modelisation et simulation de systemes multi-fluides. Application aux ecoulements sanguins., Physique Numérique [physics.comp-ph], Université de Grenoble, 2014

:leveloffset:-1

10. Computational Fluid Dynamics Benchmarks

11. Turek & Hron CFD Benchmark

11.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 Geometry of the reduce model.

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.

11.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.

TurekHron Geometry
Figure 1. Geometry of the Turek & Hron CFD Benchmark

In order to describe the flow, the incompressible Navier-Stokes 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.

11.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} \)

11.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.41-y\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\).

11.3. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 21. Fixed and Variable Input Parameters
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

CFD1 CFD2 CFD3

\(0.2\)

\(1\)

\(2\)

\(m/s\)

11.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\).

11.5. Discretization

To realize these tests, we made the choice to used \(P_N\)-\(P_{N-1}\) Taylor-Hood finite elements, described by Chabannes, to discretize space. With the time discretization, we use BDF, for Backward Differentation Formulation, schemes at different orders \(q\).

11.5.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

Table 22. KSP configuration

type

gmres

relative tolerance

1e-13

max iteration

1000

reuse preconditioner

false

Table 23. SNES configuration

relative tolerance

1e-8

steps tolerance

1e-8

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

Table 24. KSP in SNES configuration

relative tolerance

1e-5

max iteration

1000

max iteration with reuse

CFD1/CFD2 : 100 | CFD3 : 1000

reuse preconditioner

false

reuse preconditioner rebuild at first Newton step

false

Table 25. Preconditioner configuration

type

lu

package

mumps

11.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 feelpp-toolboxes docker reads

Command line to execute CFD1 testcase
$ mpirun -np 4 /usr/local/bin/feelpp_toolbox_fluid_2d --config-file cfd1.cfg

The result files are then stored by default in

Results Directory
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

11.7. Results

Here are results from the different cases studied in this benchmark.

11.7.1. CFD1

Table 26. Results for 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]

11.7.2. CFD2

Table 27. Results for 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]

11.7.3. CFD3

As CFD3 is time-dependent ( 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}(max-min)

  • 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}

Table 28. Results for CFD3
\(\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 (\(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].

TurekHron CFD3 results
Figure 2. Lift and drag forces

11.8. Geometrical Order

Add a section on geometrical order.

11.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.

11.10. Bibliography

References for this benchmark
  • [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure 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.

Unresolved directive in /var/lib/buildkite-agent/builds/sd-87660-1/feelpp/www-dot-feelpp-dot-org/pages/man/04-learning/CFD/README.adoc - include::MultiFluid/README.adoc[]

12. 2D Drops Benchmark

This benchmark has been proposed and realised by Hysing. It allows us to verify our level set code, our Navier-Stokes solver and how they couple together.

Computer codes, used for the acquisition of results, are from Vincent Doyeux, with the use of Chabannes's Navier-Stokes code.

12.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 Navier-Stokes 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.

2D Bubble Geo
Figure 3. Geometry used in 2D Bubble Benchmark

Durig this benchmark, we will study two different cases : the first one with a ellipsoidal bubble and the second one with a squirted bubble.

12.1.1. Boundary conditions

  • On the lateral walls, we imposed slip conditions

\[\begin{eqnarray} \boldsymbol{u}\cdot\boldsymbol{n}&=&0 \\ t\cdot(\nabla\boldsymbol{u}+^t\nabla\boldsymbol{u})\cdot \boldsymbol{n}&=&0 \end{eqnarray}\]
  • On the horizontal walls, no slip conditions are imposed : \(\boldsymbol{u}=0 \)

12.1.2. 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\)

12.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 29. Fixed and Variable Input Parameters
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\)

12.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

\[\boldsymbol{X}_c = \dfrac{ \displaystyle \int_{\Omega_2} \boldsymbol{x}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{ \displaystyle \int_\Omega \boldsymbol{x} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }\]
\[\boldsymbol{U}_c = \dfrac{\displaystyle \int_{\Omega_2} \boldsymbol{u}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{\displaystyle \int_\Omega \boldsymbol{u} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }\]
\[c = \dfrac{\left(4 \pi \displaystyle \int_{\Omega_2} 1 \right)^{\frac{1}{2}}}{ \displaystyle \int_{\Gamma} 1} = \dfrac{ \left(4 \pi \displaystyle \int_{\Omega} (1 - H_\varepsilon(\phi)) \right) ^{\frac{1}{2}}}{ \displaystyle \int_{\Omega} \delta_\varepsilon(\phi)}\]

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.

12.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

12.6. Results

12.6.2. Test 2

We describe the different quantitative results for the two studied cases.

Table 30. Results comparison between benchmark values and our results for the ellipsoidal case

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

Table 31. Results comparison between benchmark values and our results for the squirted case

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

12.7. Bibliography

References for this benchmark
  • [Hysing] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska, Quantitative benchmark computations of two-dimensional 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 multi-fluides, Application aux écoulements sanguins, PhD thesis, Université de Grenoble, 2014.

13. 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}.

13.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:

\[\begin{array}[lll] \rho\rho(\phi(\mathbf{x}) ) \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) + \nabla p - \nabla \cdot \left( \nu(\phi(\mathbf{x})) (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) \right) &=& \rho ( \phi(\mathbf{x}) ) \mathbf{g}, \\ \nabla \cdot \mathbf{u} &=& 0, \\ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi &=& 0, \end{array}\]

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, no-slip 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 two-fluids 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}.

Table 32. Numerical parameters taken for the benchmarks.

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}\).

13.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: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation. summarize for the test 1 the different simulation properties and the table Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom give the carachteristics of each mesh.

Table 33. Numerical parameters used for the test 1 simulations: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation.

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

Table 34. Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom

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 Navier-Stokes 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.

13.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.

13.4. Bibliography

.

References for this benchmark
  • [cottet]

  • [Feelpp] C. Prud’homme et al.

  • [osher] Osher1988, book_Sethian, book_Osher

  • [Franca1992] Franca 1992

14. Computational Solid Mechanics

Computational solid mechanics ( or CSM ) is a part of mechanics that studies solid deformations or motions under applied forces or other parameters.

14.1. Second Newton’s law

The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows

\[\rho^*_s \frac{\partial^2\boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot (\boldsymbol{F}_s\boldsymbol{\Sigma}_s) = \boldsymbol{f}^t_s\]

It’s define here into a Lagrangian frame.

14.1.1. Variables, symbols and units

Notation

Quantity

Unit

\(\rho_s^*\)

strucure density

\(kg/m^3\)

\(\boldsymbol{\eta}_s\)

displacement

\(m\)

\(\boldsymbol{F}_s\)

deformation gradient

\(\boldsymbol{\Sigma}_s\)

second Piola-Kirchhoff tensor

\(N/m^2\)

\(f_s^t\)

body force

\(N/m^2\)

14.2. Lamé coefficients

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

\[\lambda_s = \frac{E_s\nu_s}{(1+\nu_s)(1-2\nu_s)} \hspace{0.5 cm} , \hspace{0.5 cm} \mu_s = \frac{E_s}{2(1+\nu_s)}\]

14.3. 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.

Reduced Model Geometry
Figure 1 : Geometry of the reduce model

The mathematical problem associated to this reduce model can be describe as

\[ \rho^*_s h \frac{\partial^2 \eta_s}{\partial t^2} - k G_s h \frac{\partial^2 \eta_s}{\partial x^2} + \frac{E_s h}{1-\nu_s^2} \frac{\eta_s}{R_0^2} - \gamma_v \frac{\partial^3 \eta}{\partial x^2 \partial t} = f_s.\]

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.

inlude::FSI/README.adoc[]

15. Thermal Building Environment

16. ISO 10211:2007 Thermal bridges in building construction

16.1. Introduction

ISO 10211:2007 sets out the specifications for a three-dimensional and a two-dimensional geometrical model of a thermal bridge for the numerical calculation of:

  1. heat flows, in order to assess the overall heat loss from a building or part of it;

  2. minimum surface temperatures, in order to assess the risk of surface condensation.

These specifications include the geometrical boundaries and subdivisions of the model, the thermal boundary conditions, and the thermal values and relationships to be used.

ISO 10211:2007 is based upon the following assumptions:

  1. all physical properties are independent of temperature;

  2. there are no heat sources within the building element.

ISO 10211:2007 can also be used for the derivation of linear and point thermal transmittances and of surface temperature factors. More information here.

16.2. Implementation

Only the 2D specifications have been implemented.

17. Running the testcase

$ mpirun -np 4 /usr/local/bin/feelpp_thermodyn_2d --config-file thermo2dCase2.cfg

inlude::HeatFluid/README.adoc[]

18. Optimization problems

This section presents some mathematical optimization problems. The following examples source codes are located in "doc/manual/opt/".

:leveloffset:-1

19. Generic Partial Differential Equations

:leveloffset:+1

20. The Laplacian

20.1. Problem statement

We are interested in this section in the conforming finite element approximation of the following problem:

Laplacian problem

Look for \(u\) such that

\[\left\{\begin{split} -\Delta u &= f \text{ in } \Omega\\ u &= g \text{ on } \partial \Omega_D\\ \frac{\partial u}{\partial n} &=h \text{ on } \partial \Omega_N\\ \frac{\partial u}{\partial n} + u &=l \text{ on } \partial \Omega_R \end{split}\right.\]
\(\partial \Omega_D\), \(\partial \Omega_N\) and \(\partial \Omega_R\) can be empty sets. In the case \(\partial \Omega_D =\partial \Omega_R = \emptyset\), then the solution is known up to a constant.

In the implementation presented later, \(\partial \Omega_D =\partial \Omega_N = \partial \Omega_R = \emptyset\), then we set Dirichlet boundary conditions all over the boundary. The problem then reads like a standard laplacian with inhomogeneous Dirichlet boundary conditions:

Laplacian Problem with inhomogeneous Dirichlet conditions

Look for \(u\) such that

Inhomogeneous Dirichlet Laplacian problem
\[-\Delta u = f\ \text{ in } \Omega,\quad u = g \text{ on } \partial \Omega\]

20.2. Variational formulation

We assume that \(f, h, l \in L^2(\Omega)\). The weak formulation of the problem then reads:

Laplacian problem variational formulation

Look for \(u \in H^1_{g,\Gamma_D}(\Omega)\) such that

Variational formulation
\[\displaystyle\int_\Omega \nabla u \cdot \nabla v +\int_{\Gamma_R} u v = \displaystyle \int_\Omega f\ v+ \int_{\Gamma_N} g\ v + \int_{\Gamma_R} l\ v,\quad \forall v \in H^1_{0,\Gamma_D}(\Omega)\]

20.3. Conforming Approximation

We now turn to the finite element approximation using Lagrange finite element. We assume \(\Omega\) to be a segment in 1D, a polygon in 2D or a polyhedron in 3D. We denote \(V_\delta \subset H^1(\Omega)\) an approximation space such that \(V_{g,\delta} \equiv P^k_{c,\delta}\cap H^1_{g,\Gamma_D}(\Omega)\).

The weak formulation reads:

Laplacian problem weak formulation

Look for \(u_\delta \in V_\delta \) such that

\[\displaystyle\int_{\Omega_\delta} \nabla u_{\delta} \cdot \nabla v_\delta +\int_{\Gamma_{R,\delta}} u_\delta\ v_\delta = \displaystyle \int_{\Omega_\delta} f\ v_\delta+ \int_{\Gamma_{N,\delta}} g\ v_\delta + \int_{\Gamma_{R,\delta}} l\ v_\delta,\quad \forall v_\delta \in V_{0,\delta}\]
from now on, we omit \(\delta\) to lighten the notations. Be careful that it appears both the geometrical and approximation level.

20.4. Feel++ Implementation

In Feel++, \(V_{g,\delta}\) is not built but rather \(P^k_{c,\delta}\).

The Dirichlet boundary conditions can be treated using different techniques and we use from now on the elimination technique.

We start with the mesh

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
the keyword auto enables type inference, for more details see Wikipedia C++11 page.

Next the discretization setting by first defining Vh=Pch<k>(mesh) \(\equiv P^k_{c,h}\), then elements of Vh and expressions f, n and g given by command line options or configuration file.

    auto Vh = Pch<2>( mesh );
    auto u = Vh->element("u");
    auto mu = doption(_name="mu");
    auto f = expr( soption(_name="functions.f"), "f" );
    auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
    auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
    auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
    auto g = expr( soption(_name="functions.g"), "g" );
    auto v = Vh->element( g, "g" );

at the following line

    auto v = Vh->element( g, "g" );

v is set to the expression g, which means more precisely that v is the interpolant of g in Vh.

the variational formulation is implemented below, we define the bilinear form a and linear form l and we set strongly the Dirichlet boundary conditions with the keyword on using elimination. If we don’t find Dirichlet, Neumann or Robin in the list of physical markers in the mesh data structure then we impose Dirichlet boundary conditions all over the boundary.

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=f*id(v));
    l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
    l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
    toc("l");

    tic();
    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=mu*gradt(u)*trans(grad(v)) );
    a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
    a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
    //! if no markers Robin Neumann or Dirichlet are present in the mesh then
    //! impose Dirichlet boundary conditions over the entire boundary
    if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
        a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
    toc("a");

    tic();
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);
    toc("a.solve");

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

    tic();
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

    toc("Exporter");
    return 0;

}

We have the following correspondance:

Element sets Domain

elements(mesh)

\(\Omega\)

boundaryfaces(mesh)

\(\partial \Omega\)

markedfaces(mesh,"Dirichlet")

\(\Gamma_D\)

markedfaces(mesh,"Neumann")

\(\Gamma_R\)

markedfaces(mesh,"Robin")

\(\Gamma_R\)

next we solve the algebraic problem

Listing: solve algebraic system
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);

next we compute the \(L^2\) norm of \(u_\delta-g\), it could serve as an \(L^2\) error if \(g\) was manufactured to be the exact solution of the Laplacian problem.

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

and finally we export the results, by default it is in the ensight gold format and the files can be read with Paraview and Ensight. We save both \(u\) and \(g\).

Listing: export Laplacian results
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

20.5. Testcases

The Feel++ Implementation comes with testcases in 2D and 3D.

20.5.1. circle

circle is a 2D testcase where \(\Omega\) is a disk whose boundary has been split such that \(\partial \Omega=\partial \Omega_D \cup \partial \Omega_N \cup \partial \Omega_R\).

Here are some results we can observe after use the following command

cd Testcases/quickstart/circle
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file circle.cfg

This give us some data such as solution of our problem or the mesh used in the application.

ucircle

meshCircle

Solution \(u_\delta\)

Mesh

20.5.2. feelpp2d and feelpp3d

This testcase solves the Laplacian problem in \(\Omega\) an quadrangle or hexadra containing the letters of Feel++

feelpp2d

After running the following command

cd Testcases/quickstart/feelpp2d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file feelpp2d.cfg

we obtain the result \(u_\delta\) and also the mesh

ufeelpp2d

/images/Laplacian/TestCases/Feelpp2d/meshfeelpp2d.png[]

Solution \(u_\delta\)

Mesh

feelpp3d

We can launch this application with the current line

cd Testcases/quickstart/feelpp3d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_3d --config-file feelpp3d.cfg

When it’s finish, we can extract some informations

ufeelpp3d

meshfeelpp3d

Solution \(u_\delta\)

Mesh

21. Levelset

Having the possibility to determine where two regions meeting can be really useful in some scientific domains. That’s why the levelset method is born.

21.1. Levelset introduction

21.1.1. Levelset function

By using a scalar function \phi, define on all regions as the null value is obtained when it’s placed on an interface of two domains.

We denote \Omega_1 and \Omega_2 two domains with \Gamma the interface betwen them. Then \phi can be define as

\phi(\boldsymbol{x}) = \left\{ \begin{array}{cccc} \text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in &\Omega_1 \\ 0, & \boldsymbol{x}& \in &\Gamma\\ -\text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in & \Omega_2 \end{array} \right.

with \text{dist}(\boldsymbol{x}, \Gamma ) = \underset{\boldsymbol{y} \; \in \; \Gamma}{\min}( |\boldsymbol{x} - \boldsymbol{y}| ).

This function \phi had also the following property |\nabla\phi|=1.

Moreover, the unit normal vector \boldsymbol{n} outgoing from the interface and the curvature \mathcal{\kappa} can be obtained from the levelset function.

\boldsymbol{n}=\frac{\nabla\phi}{|\nabla\phi} \\ \mathcal{\kappa}=\nabla \cdot \boldsymbol{n}= \nabla \cdot \frac{\nabla\phi}{|\nabla\phi|}

Now we have exposed the levelset function, we need to define how the levelset will evolve and will spread into all the space. To do this, we use the following advection equation :

\partial_t\phi+\boldsymbol{u}\cdot\nabla\phi=0

where \boldsymbol{u} is an incompressible velocity field.

21.1.2. Heaviside and Dirac functions

We define also the regularized Heaviside function H_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \dfrac{1}{2} \left(1+\dfrac{\phi}{\varepsilon}+\dfrac{\sin\left(\dfrac{\pi \phi}{\varepsilon}\right)}{\pi}\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 1, & \phi \geq \varepsilon \end{array} \right.

and the regularized Dirac function \delta_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \displaystyle\dfrac{1}{2 \varepsilon} \left(1+\cos\left(\dfrac{\pi \phi}{\varepsilon}\right)\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 0, & \phi \geq \varepsilon. \end{array} \right.

The first one gives a different value to each side of the interface ( here 0 in and 1 out ). The second one allow us to define quantities, with value different from 0 at the interface. A typical value of \varepsilon in literature is 1.5h where h is the mesh step size.

It should be noted that these functions allow us to determine respectively the volume and the surface of the interface by V^+_{\varepsilon} = \int_{\Omega} H_\varepsilon \\ S^{\Gamma}_{\varepsilon} = \int_{\Omega} \delta_\varepsilon

22. Solid rotation of a slotted disk

We describe the benchmark proposed by Zalesak.

Computer codes, used for the acquisition of results, are from Vincent Doyeux.

22.1. Problem Description

In order to test our interface propagation method, i.e. the levelset method \phi, we will study the rotation of a slotted disk into a square domain. The geometry can be represented as

Slotted Disk Geometry
Figure 1 : Initial Geometry.

We denote \Omega, the square domain [0,1]\times[0,1]. The center of the slotted disk is placed at (0.5,0.75).
To model the rotation, we will apply an angular velocity, centered in (0.5,0.5), as the disk is back to its initial position after t_f=628.

During this test, we observe three different errors to measure the quality of our method. With these values, two kinds of convergence will be studied : the time convergence, with different time step on an imposed grid and the space one, where the space discretization and the time step are linked by a relation. Several stabilization methods are used such as CIP ( Continuous Interior Penalty ) or SUPG ( Streamline-Upwing/Petrov-Galerkin ).

22.1.1. Boundary conditions

We set a Neumann boundary condition on the boundary of the domain.

22.1.2. Initial conditions

The velocity is imposed as \boldsymbol{u}=\left( \frac{\pi}{314} (50-y),\frac{\pi}{314} (x-50) \right)

Here is the velocity look in the square domain

Velocity Geometry
Figure 2 : Imposed velocity in \Omega.

22.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 35. Fixed and Variable Input Parameters

Name

Description

Nominal Value

Units

r

disk radius

0.15

m

l

slot base

0.05

m

h

slot height

0.25

m

t_f

slotted disk rotation period

628

s

22.3. Outputs

We observe during this benchmark three different errors.
First at all, the mass error, define by e_{\text{m}} = \frac{ \left| m_{\phi_f} - m_{\phi_0}\right| }{m_{\phi_0}} = \frac{ \left| \displaystyle \int_{\Omega} \chi( \phi_f < 0 ) - \displaystyle \int_{\Omega} \chi( \phi_0 < 0 ) \right| }{ \displaystyle \int_{\Omega} \chi( \phi_0 < 0 )}

where \chi is the characteristic function. \chi( f( \phi ) ) = \left\{ \begin{array}{rcl} 1 & \text{ if } & f( \phi ) \neq 0 \\ 0 & \text{ if } & f( \phi ) = 0 \end{array} \right.

However, mass is gain and loose at different emplacements on the mesh, and at the same time, with the level set method.

Secondly, the sign change error e_{\text{sc}} = \sqrt{ \int_\Omega \left( (1-H_0) - (1-H_f) \right)^2 }

with H_0=H_\epsilon(\phi_0) and H_f=H_\epsilon(\phi_f), H_\epsilon the smoothed Heaviside function of thickness 2ε.

This error is better to define the interface displacement. In fact, we can determine where \phi_0\phi_f<0, in other words where the interface has moved.

Finally, we define the classical L^2 error at the interface, as e_{L^2} = \sqrt{ \frac{1}{\displaystyle \int_\Omega \chi( \delta(\phi_0) > 0 ) } \int_\Omega (\phi_0 - \phi_f)^2 \chi( \delta(\phi_0) > 0 ) }.

22.4. Discretization

22.4.1. Time convergence

For this case, we set a fixed grid with mesh step size h=0.04, and so 72314 degree of freedom on a \mathbb{P}^1.

Then, after the disk made one round, we measure the errors obtained from two different discretizations ( BDF2 and Euler ) and compared them.

We repeat this with several time step dt\in \{2.14, 1, 0.5, 0.25, 0.20\}.

Only one stabilization method is used : SUPG

22.4.2. Space convergence

We define the following relation, between time step and mesh step size : dt=C\frac{h}{U_{max}}

where C<1 constant and U_{max} the maximum velocity of \Omega.

From the definition of our velocity, U_{max} is reached at the farthest point from the center of \Omega. In this case, we have U_{max}=0.007, and we set C=0.8.

We use the BDF2 method for time discretization. As in time convergence, we wait one round of the disk to measure the errors and we repeat this test for these values of h: 0.32, 0.16, 0.08, 0.04.

We compare the results from different stabilization methods : CIP, SUPG, GLS ( Galerkin-Least-Squares ) and SGS ( Sub-Grid Scale ).

22.6. Results

22.6.1. Time convergence

Table 36. Time convergence with Euler scheme

dt

e_{L^2}

e_{sc}

e_m

2.14

0.0348851

0.202025

0.202025

1.00

0.0187567

0.147635

0.147635

0.5

0.0098661

0.10847

0.10847

0.25

0.008791

0.0782569

0.0782569

0.20

0.00803373

0.0670677

0.0670677

Table 37. Time convergence with BDF2 scheme

dt

e_{L^2}

e_{sc}

e_m

2.14

0.0118025

0.0906617

0.0492775

1.00

0.00436957

0.0445275

0.0163494

0.5

0.00173637

0.0216359

0.0100621

0.25

0.001003

0.0125971

0.00354644

0.20

0.000949343

0.0117449

0.00317368

Time mass error
Figure 2 : Mass error of Zalesac benchmark
Time sign change
Figure 3 : Sign change error of Zalesac benchmark
Time L2 error
Figure 4 : L^2 error of Zalesac benchmark
Time shape
Figure 5 : Slotted disk shape after a round

22.6.2. Space convergence

stab

h

e_{L^2}

e_{sc}

e_m

CIP

0.32

0.0074

0.072

0.00029

0.16

0.0046

0.055

0.00202

0.08

0.0025

0.033

0.00049

0.04

0.0023

0.020

0.00110

SUPG

0.32

0.012

0.065

0.01632

0.16

0.008

0.049

0.07052

0.08

0.004

0.030

0.00073

0.04

0.001

0.018

0.00831

GLS

0.32

0.013

0.066

0.02499

0.16

0.008

0.051

0.05180

0.08

0.004

0.031

0.00805

0.04

0.001

0.019

0.00672

SGS

0.32

0.012

0.065

0.01103

0.16

0.008

0.050

0.07570

0.08

0.004

0.030

0.00084

0.04

0.001

0.018

0.00850

Space sign change
Figure 6 : Sign change error of Zalesac benchmark
Space L2 error
Figure 7 : L^2 error of Zalesac benchmark
Space Shape
Figure 8 : Slotted disk shape after a round

22.6.3. Conclusion

Let’s begin with time convergence results. Tables shows us that sign change error is better to define the quality of the chosen scheme than the mass error. In fact, the loss of mass somewhere can be nullify by a gain of mass elsewhere. Sign change error shows half an order gain from Euler scheme to BDF2 one, as L^2 errors show us a gain of one order. For the slotted disk shape, BDF2 uses the two previous iterations to obtain the current result, while Euler only need the previous iteration. This explain why we can see an asymmetrical tendency in the first one.

As for space convergence, the different stabilization methods we used give us the same convergence rate equals to 0.6, with close error values, for the sign change error. For the L^2 error case, it’s not as evident as the previous one. Aside the CIP stabilization method with a 0.6 convergence rate, the others show us a convergence rate of 0.9.

22.7. Bibliography

References for this benchmark
  • [Zalesak] Steven T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, Journal of Computational Physics, 1979.

  • [Doyeux] Vincent Doyeux, Modelisation et simulation de systemes multi-fluides. Application aux ecoulements sanguins., Physique Numérique [physics.comp-ph], Université de Grenoble, 2014

:leveloffset:-1

23. The Laplacian

23.1. Problem statement

We are interested in this section in the conforming finite element approximation of the following problem:

Laplacian problem

Look for \(u\) such that

\[\left\{\begin{split} -\Delta u &= f \text{ in } \Omega\\ u &= g \text{ on } \partial \Omega_D\\ \frac{\partial u}{\partial n} &=h \text{ on } \partial \Omega_N\\ \frac{\partial u}{\partial n} + u &=l \text{ on } \partial \Omega_R \end{split}\right.\]
\(\partial \Omega_D\), \(\partial \Omega_N\) and \(\partial \Omega_R\) can be empty sets. In the case \(\partial \Omega_D =\partial \Omega_R = \emptyset\), then the solution is known up to a constant.

In the implementation presented later, \(\partial \Omega_D =\partial \Omega_N = \partial \Omega_R = \emptyset\), then we set Dirichlet boundary conditions all over the boundary. The problem then reads like a standard laplacian with inhomogeneous Dirichlet boundary conditions:

Laplacian Problem with inhomogeneous Dirichlet conditions

Look for \(u\) such that

Inhomogeneous Dirichlet Laplacian problem
\[-\Delta u = f\ \text{ in } \Omega,\quad u = g \text{ on } \partial \Omega\]

23.2. Variational formulation

We assume that \(f, h, l \in L^2(\Omega)\). The weak formulation of the problem then reads:

Laplacian problem variational formulation

Look for \(u \in H^1_{g,\Gamma_D}(\Omega)\) such that

Variational formulation
\[\displaystyle\int_\Omega \nabla u \cdot \nabla v +\int_{\Gamma_R} u v = \displaystyle \int_\Omega f\ v+ \int_{\Gamma_N} g\ v + \int_{\Gamma_R} l\ v,\quad \forall v \in H^1_{0,\Gamma_D}(\Omega)\]

23.3. Conforming Approximation

We now turn to the finite element approximation using Lagrange finite element. We assume \(\Omega\) to be a segment in 1D, a polygon in 2D or a polyhedron in 3D. We denote \(V_\delta \subset H^1(\Omega)\) an approximation space such that \(V_{g,\delta} \equiv P^k_{c,\delta}\cap H^1_{g,\Gamma_D}(\Omega)\).

The weak formulation reads:

Laplacian problem weak formulation

Look for \(u_\delta \in V_\delta \) such that

\[\displaystyle\int_{\Omega_\delta} \nabla u_{\delta} \cdot \nabla v_\delta +\int_{\Gamma_{R,\delta}} u_\delta\ v_\delta = \displaystyle \int_{\Omega_\delta} f\ v_\delta+ \int_{\Gamma_{N,\delta}} g\ v_\delta + \int_{\Gamma_{R,\delta}} l\ v_\delta,\quad \forall v_\delta \in V_{0,\delta}\]
from now on, we omit \(\delta\) to lighten the notations. Be careful that it appears both the geometrical and approximation level.

23.4. Feel++ Implementation

In Feel++, \(V_{g,\delta}\) is not built but rather \(P^k_{c,\delta}\).

The Dirichlet boundary conditions can be treated using different techniques and we use from now on the elimination technique.

We start with the mesh

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
the keyword auto enables type inference, for more details see Wikipedia C++11 page.

Next the discretization setting by first defining Vh=Pch<k>(mesh) \(\equiv P^k_{c,h}\), then elements of Vh and expressions f, n and g given by command line options or configuration file.

    auto Vh = Pch<2>( mesh );
    auto u = Vh->element("u");
    auto mu = doption(_name="mu");
    auto f = expr( soption(_name="functions.f"), "f" );
    auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
    auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
    auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
    auto g = expr( soption(_name="functions.g"), "g" );
    auto v = Vh->element( g, "g" );

at the following line

    auto v = Vh->element( g, "g" );

v is set to the expression g, which means more precisely that v is the interpolant of g in Vh.

the variational formulation is implemented below, we define the bilinear form a and linear form l and we set strongly the Dirichlet boundary conditions with the keyword on using elimination. If we don’t find Dirichlet, Neumann or Robin in the list of physical markers in the mesh data structure then we impose Dirichlet boundary conditions all over the boundary.

    auto l = form1( _test=Vh );
    l = integrate(_range=elements(mesh),
                  _expr=f*id(v));
    l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
    l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
    toc("l");

    tic();
    auto a = form2( _trial=Vh, _test=Vh);
    a = integrate(_range=elements(mesh),
                  _expr=mu*gradt(u)*trans(grad(v)) );
    a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
    a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
    //! if no markers Robin Neumann or Dirichlet are present in the mesh then
    //! impose Dirichlet boundary conditions over the entire boundary
    if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
        a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
    toc("a");

    tic();
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);
    toc("a.solve");

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

    tic();
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

    toc("Exporter");
    return 0;

}

We have the following correspondance:

Element sets Domain

elements(mesh)

\(\Omega\)

boundaryfaces(mesh)

\(\partial \Omega\)

markedfaces(mesh,"Dirichlet")

\(\Gamma_D\)

markedfaces(mesh,"Neumann")

\(\Gamma_R\)

markedfaces(mesh,"Robin")

\(\Gamma_R\)

next we solve the algebraic problem

Listing: solve algebraic system
    //! solve the linear system, find u s.t. a(u,v)=l(v) for all v
    if ( !boption( "no-solve" ) )
        a.solve(_rhs=l,_solution=u);

next we compute the \(L^2\) norm of \(u_\delta-g\), it could serve as an \(L^2\) error if \(g\) was manufactured to be the exact solution of the Laplacian problem.

    cout << "||u_h-g||_L2=" << normL2(_range=elements(mesh), _expr=idv(u)-g) << std::endl;

and finally we export the results, by default it is in the ensight gold format and the files can be read with Paraview and Ensight. We save both \(u\) and \(g\).

Listing: export Laplacian results
    auto e = exporter( _mesh=mesh );
    e->addRegions();
    e->add( "u", u );
    e->add( "g", v );
    e->save();

23.5. Testcases

The Feel++ Implementation comes with testcases in 2D and 3D.

23.5.1. circle

circle is a 2D testcase where \(\Omega\) is a disk whose boundary has been split such that \(\partial \Omega=\partial \Omega_D \cup \partial \Omega_N \cup \partial \Omega_R\).

Here are some results we can observe after use the following command

cd Testcases/quickstart/circle
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file circle.cfg

This give us some data such as solution of our problem or the mesh used in the application.

ucircle

meshCircle

Solution \(u_\delta\)

Mesh

23.5.2. feelpp2d and feelpp3d

This testcase solves the Laplacian problem in \(\Omega\) an quadrangle or hexadra containing the letters of Feel++

feelpp2d

After running the following command

cd Testcases/quickstart/feelpp2d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_2d --config-file feelpp2d.cfg

we obtain the result \(u_\delta\) and also the mesh

ufeelpp2d

/images/Laplacian/TestCases/Feelpp2d/meshfeelpp2d.png[]

Solution \(u_\delta\)

Mesh

feelpp3d

We can launch this application with the current line

cd Testcases/quickstart/feelpp3d
mpirun -np 4 /usr/local/bin/feelpp_qs_laplacian_3d --config-file feelpp3d.cfg

When it’s finish, we can extract some informations

ufeelpp3d

meshfeelpp3d

Solution \(u_\delta\)

Mesh

24. Levelset

Having the possibility to determine where two regions meeting can be really useful in some scientific domains. That’s why the levelset method is born.

24.1. Levelset introduction

24.1.1. Levelset function

By using a scalar function \phi, define on all regions as the null value is obtained when it’s placed on an interface of two domains.

We denote \Omega_1 and \Omega_2 two domains with \Gamma the interface betwen them. Then \phi can be define as

\phi(\boldsymbol{x}) = \left\{ \begin{array}{cccc} \text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in &\Omega_1 \\ 0, & \boldsymbol{x}& \in &\Gamma\\ -\text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in & \Omega_2 \end{array} \right.

with \text{dist}(\boldsymbol{x}, \Gamma ) = \underset{\boldsymbol{y} \; \in \; \Gamma}{\min}( |\boldsymbol{x} - \boldsymbol{y}| ).

This function \phi had also the following property |\nabla\phi|=1.

Moreover, the unit normal vector \boldsymbol{n} outgoing from the interface and the curvature \mathcal{\kappa} can be obtained from the levelset function.

\boldsymbol{n}=\frac{\nabla\phi}{|\nabla\phi} \\ \mathcal{\kappa}=\nabla \cdot \boldsymbol{n}= \nabla \cdot \frac{\nabla\phi}{|\nabla\phi|}

Now we have exposed the levelset function, we need to define how the levelset will evolve and will spread into all the space. To do this, we use the following advection equation :

\partial_t\phi+\boldsymbol{u}\cdot\nabla\phi=0

where \boldsymbol{u} is an incompressible velocity field.

24.1.2. Heaviside and Dirac functions

We define also the regularized Heaviside function H_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \dfrac{1}{2} \left(1+\dfrac{\phi}{\varepsilon}+\dfrac{\sin\left(\dfrac{\pi \phi}{\varepsilon}\right)}{\pi}\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 1, & \phi \geq \varepsilon \end{array} \right.

and the regularized Dirac function \delta_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \displaystyle\dfrac{1}{2 \varepsilon} \left(1+\cos\left(\dfrac{\pi \phi}{\varepsilon}\right)\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 0, & \phi \geq \varepsilon. \end{array} \right.

The first one gives a different value to each side of the interface ( here 0 in and 1 out ). The second one allow us to define quantities, with value different from 0 at the interface. A typical value of \varepsilon in literature is 1.5h where h is the mesh step size.

It should be noted that these functions allow us to determine respectively the volume and the surface of the interface by V^+_{\varepsilon} = \int_{\Omega} H_\varepsilon \\ S^{\Gamma}_{\varepsilon} = \int_{\Omega} \delta_\varepsilon

25. Solid rotation of a slotted disk

We describe the benchmark proposed by Zalesak.

Computer codes, used for the acquisition of results, are from Vincent Doyeux.

25.1. Problem Description

In order to test our interface propagation method, i.e. the levelset method \phi, we will study the rotation of a slotted disk into a square domain. The geometry can be represented as

Slotted Disk Geometry
Figure 1 : Initial Geometry.

We denote \Omega, the square domain [0,1]\times[0,1]. The center of the slotted disk is placed at (0.5,0.75).
To model the rotation, we will apply an angular velocity, centered in (0.5,0.5), as the disk is back to its initial position after t_f=628.

During this test, we observe three different errors to measure the quality of our method. With these values, two kinds of convergence will be studied : the time convergence, with different time step on an imposed grid and the space one, where the space discretization and the time step are linked by a relation. Several stabilization methods are used such as CIP ( Continuous Interior Penalty ) or SUPG ( Streamline-Upwing/Petrov-Galerkin ).

25.1.1. Boundary conditions

We set a Neumann boundary condition on the boundary of the domain.

25.1.2. Initial conditions

The velocity is imposed as \boldsymbol{u}=\left( \frac{\pi}{314} (50-y),\frac{\pi}{314} (x-50) \right)

Here is the velocity look in the square domain

Velocity Geometry
Figure 2 : Imposed velocity in \Omega.

25.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 38. Fixed and Variable Input Parameters

Name

Description

Nominal Value

Units

r

disk radius

0.15

m

l

slot base

0.05

m

h

slot height

0.25

m

t_f

slotted disk rotation period

628

s

25.3. Outputs

We observe during this benchmark three different errors.
First at all, the mass error, define by e_{\text{m}} = \frac{ \left| m_{\phi_f} - m_{\phi_0}\right| }{m_{\phi_0}} = \frac{ \left| \displaystyle \int_{\Omega} \chi( \phi_f < 0 ) - \displaystyle \int_{\Omega} \chi( \phi_0 < 0 ) \right| }{ \displaystyle \int_{\Omega} \chi( \phi_0 < 0 )}

where \chi is the characteristic function. \chi( f( \phi ) ) = \left\{ \begin{array}{rcl} 1 & \text{ if } & f( \phi ) \neq 0 \\ 0 & \text{ if } & f( \phi ) = 0 \end{array} \right.

However, mass is gain and loose at different emplacements on the mesh, and at the same time, with the level set method.

Secondly, the sign change error e_{\text{sc}} = \sqrt{ \int_\Omega \left( (1-H_0) - (1-H_f) \right)^2 }

with H_0=H_\epsilon(\phi_0) and H_f=H_\epsilon(\phi_f), H_\epsilon the smoothed Heaviside function of thickness 2ε.

This error is better to define the interface displacement. In fact, we can determine where \phi_0\phi_f<0, in other words where the interface has moved.

Finally, we define the classical L^2 error at the interface, as e_{L^2} = \sqrt{ \frac{1}{\displaystyle \int_\Omega \chi( \delta(\phi_0) > 0 ) } \int_\Omega (\phi_0 - \phi_f)^2 \chi( \delta(\phi_0) > 0 ) }.

25.4. Discretization

25.4.1. Time convergence

For this case, we set a fixed grid with mesh step size h=0.04, and so 72314 degree of freedom on a \mathbb{P}^1.

Then, after the disk made one round, we measure the errors obtained from two different discretizations ( BDF2 and Euler ) and compared them.

We repeat this with several time step dt\in \{2.14, 1, 0.5, 0.25, 0.20\}.

Only one stabilization method is used : SUPG

25.4.2. Space convergence

We define the following relation, between time step and mesh step size : dt=C\frac{h}{U_{max}}

where C<1 constant and U_{max} the maximum velocity of \Omega.

From the definition of our velocity, U_{max} is reached at the farthest point from the center of \Omega. In this case, we have U_{max}=0.007, and we set C=0.8.

We use the BDF2 method for time discretization. As in time convergence, we wait one round of the disk to measure the errors and we repeat this test for these values of h: 0.32, 0.16, 0.08, 0.04.

We compare the results from different stabilization methods : CIP, SUPG, GLS ( Galerkin-Least-Squares ) and SGS ( Sub-Grid Scale ).

25.6. Results

25.6.1. Time convergence

Table 39. Time convergence with Euler scheme

dt

e_{L^2}

e_{sc}

e_m

2.14

0.0348851

0.202025

0.202025

1.00

0.0187567

0.147635

0.147635

0.5

0.0098661

0.10847

0.10847

0.25

0.008791

0.0782569

0.0782569

0.20

0.00803373

0.0670677

0.0670677

Table 40. Time convergence with BDF2 scheme

dt

e_{L^2}

e_{sc}

e_m

2.14

0.0118025

0.0906617

0.0492775

1.00

0.00436957

0.0445275

0.0163494

0.5

0.00173637

0.0216359

0.0100621

0.25

0.001003

0.0125971

0.00354644

0.20

0.000949343

0.0117449

0.00317368

Time mass error
Figure 2 : Mass error of Zalesac benchmark
Time sign change
Figure 3 : Sign change error of Zalesac benchmark
Time L2 error
Figure 4 : L^2 error of Zalesac benchmark
Time shape
Figure 5 : Slotted disk shape after a round

25.6.2. Space convergence

stab

h

e_{L^2}

e_{sc}

e_m

CIP

0.32

0.0074

0.072

0.00029

0.16

0.0046

0.055

0.00202

0.08

0.0025

0.033

0.00049

0.04

0.0023

0.020

0.00110

SUPG

0.32

0.012

0.065

0.01632

0.16

0.008

0.049

0.07052

0.08

0.004

0.030

0.00073

0.04

0.001

0.018

0.00831

GLS

0.32

0.013

0.066

0.02499

0.16

0.008

0.051

0.05180

0.08

0.004

0.031

0.00805

0.04

0.001

0.019

0.00672

SGS

0.32

0.012

0.065

0.01103

0.16

0.008

0.050

0.07570

0.08

0.004

0.030

0.00084

0.04

0.001

0.018

0.00850

Space sign change
Figure 6 : Sign change error of Zalesac benchmark
Space L2 error
Figure 7 : L^2 error of Zalesac benchmark
Space Shape
Figure 8 : Slotted disk shape after a round

25.6.3. Conclusion

Let’s begin with time convergence results. Tables shows us that sign change error is better to define the quality of the chosen scheme than the mass error. In fact, the loss of mass somewhere can be nullify by a gain of mass elsewhere. Sign change error shows half an order gain from Euler scheme to BDF2 one, as L^2 errors show us a gain of one order. For the slotted disk shape, BDF2 uses the two previous iterations to obtain the current result, while Euler only need the previous iteration. This explain why we can see an asymmetrical tendency in the first one.

As for space convergence, the different stabilization methods we used give us the same convergence rate equals to 0.6, with close error values, for the sign change error. For the L^2 error case, it’s not as evident as the previous one. Aside the CIP stabilization method with a 0.6 convergence rate, the others show us a convergence rate of 0.9.

25.7. Bibliography

References for this benchmark
  • [Zalesak] Steven T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, Journal of Computational Physics, 1979.

  • [Doyeux] Vincent Doyeux, Modelisation et simulation de systemes multi-fluides. Application aux ecoulements sanguins., Physique Numérique [physics.comp-ph], Université de Grenoble, 2014

26. Computational Fluid Dynamics Benchmarks

27. Turek & Hron CFD Benchmark

27.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 Geometry of the reduce model.

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.

27.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.

TurekHron Geometry
Figure 4. Geometry of the Turek & Hron CFD Benchmark

In order to describe the flow, the incompressible Navier-Stokes 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.

27.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} \)

27.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.41-y\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\).

27.3. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 41. Fixed and Variable Input Parameters
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

CFD1 CFD2 CFD3

\(0.2\)

\(1\)

\(2\)

\(m/s\)

27.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\).

27.5. Discretization

To realize these tests, we made the choice to used \(P_N\)-\(P_{N-1}\) Taylor-Hood finite elements, described by Chabannes, to discretize space. With the time discretization, we use BDF, for Backward Differentation Formulation, schemes at different orders \(q\).

27.5.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

Table 42. KSP configuration

type

gmres

relative tolerance

1e-13

max iteration

1000

reuse preconditioner

false

Table 43. SNES configuration

relative tolerance

1e-8

steps tolerance

1e-8

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

Table 44. KSP in SNES configuration

relative tolerance

1e-5

max iteration

1000

max iteration with reuse

CFD1/CFD2 : 100 | CFD3 : 1000

reuse preconditioner

false

reuse preconditioner rebuild at first Newton step

false

Table 45. Preconditioner configuration

type

lu

package

mumps

27.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 feelpp-toolboxes docker reads

Command line to execute CFD1 testcase
$ mpirun -np 4 /usr/local/bin/feelpp_toolbox_fluid_2d --config-file cfd1.cfg

The result files are then stored by default in

Results Directory
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

27.7. Results

Here are results from the different cases studied in this benchmark.

27.7.1. CFD1

Table 46. Results for 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]

27.7.2. CFD2

Table 47. Results for 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]

27.7.3. CFD3

As CFD3 is time-dependent ( 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}(max-min)

  • 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}

Table 48. Results for CFD3
\(\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 (\(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].

TurekHron CFD3 results
Figure 5. Lift and drag forces

27.8. Geometrical Order

Add a section on geometrical order.

27.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.

27.10. Bibliography

References for this benchmark
  • [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure 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.

Unresolved directive in /var/lib/buildkite-agent/builds/sd-87660-1/feelpp/www-dot-feelpp-dot-org/pages/man/04-learning/CFD/README.adoc - include::MultiFluid/README.adoc[]

28. 2D Drops Benchmark

This benchmark has been proposed and realised by Hysing. It allows us to verify our level set code, our Navier-Stokes solver and how they couple together.

Computer codes, used for the acquisition of results, are from Vincent Doyeux, with the use of Chabannes's Navier-Stokes code.

28.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 Navier-Stokes 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.

2D Bubble Geo
Figure 6. Geometry used in 2D Bubble Benchmark

Durig this benchmark, we will study two different cases : the first one with a ellipsoidal bubble and the second one with a squirted bubble.

28.1.1. Boundary conditions

  • On the lateral walls, we imposed slip conditions

\[\begin{eqnarray} \boldsymbol{u}\cdot\boldsymbol{n}&=&0 \\ t\cdot(\nabla\boldsymbol{u}+^t\nabla\boldsymbol{u})\cdot \boldsymbol{n}&=&0 \end{eqnarray}\]
  • On the horizontal walls, no slip conditions are imposed : \(\boldsymbol{u}=0 \)

28.1.2. 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\)

28.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 49. Fixed and Variable Input Parameters
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\)

28.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

\[\boldsymbol{X}_c = \dfrac{ \displaystyle \int_{\Omega_2} \boldsymbol{x}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{ \displaystyle \int_\Omega \boldsymbol{x} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }\]
\[\boldsymbol{U}_c = \dfrac{\displaystyle \int_{\Omega_2} \boldsymbol{u}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{\displaystyle \int_\Omega \boldsymbol{u} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }\]
\[c = \dfrac{\left(4 \pi \displaystyle \int_{\Omega_2} 1 \right)^{\frac{1}{2}}}{ \displaystyle \int_{\Gamma} 1} = \dfrac{ \left(4 \pi \displaystyle \int_{\Omega} (1 - H_\varepsilon(\phi)) \right) ^{\frac{1}{2}}}{ \displaystyle \int_{\Omega} \delta_\varepsilon(\phi)}\]

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.

28.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

28.6. Results

28.6.2. Test 2

We describe the different quantitative results for the two studied cases.

Table 50. Results comparison between benchmark values and our results for the ellipsoidal case

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

Table 51. Results comparison between benchmark values and our results for the squirted case

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

28.7. Bibliography

References for this benchmark
  • [Hysing] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska, Quantitative benchmark computations of two-dimensional 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 multi-fluides, Application aux écoulements sanguins, PhD thesis, Université de Grenoble, 2014.

29. 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}.

29.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:

\[\begin{array}[lll] \rho\rho(\phi(\mathbf{x}) ) \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) + \nabla p - \nabla \cdot \left( \nu(\phi(\mathbf{x})) (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) \right) &=& \rho ( \phi(\mathbf{x}) ) \mathbf{g}, \\ \nabla \cdot \mathbf{u} &=& 0, \\ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi &=& 0, \end{array}\]

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, no-slip 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 two-fluids 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}.

Table 52. Numerical parameters taken for the benchmarks.

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}\).

29.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: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation. summarize for the test 1 the different simulation properties and the table Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom give the carachteristics of each mesh.

Table 53. Numerical parameters used for the test 1 simulations: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation.

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

Table 54. Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom

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 Navier-Stokes 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.

29.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.

29.4. Bibliography

.

References for this benchmark
  • [cottet]

  • [Feelpp] C. Prud’homme et al.

  • [osher] Osher1988, book_Sethian, book_Osher

  • [Franca1992] Franca 1992

30. Turek & Hron CFD Benchmark

30.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 Geometry of the reduce model.

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.

30.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.

TurekHron Geometry
Figure 7. Geometry of the Turek & Hron CFD Benchmark

In order to describe the flow, the incompressible Navier-Stokes 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.

30.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} \)

30.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.41-y\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\).

30.3. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 55. Fixed and Variable Input Parameters
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

CFD1 CFD2 CFD3

\(0.2\)

\(1\)

\(2\)

\(m/s\)

30.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\).

30.5. Discretization

To realize these tests, we made the choice to used \(P_N\)-\(P_{N-1}\) Taylor-Hood finite elements, described by Chabannes, to discretize space. With the time discretization, we use BDF, for Backward Differentation Formulation, schemes at different orders \(q\).

30.5.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

Table 56. KSP configuration

type

gmres

relative tolerance

1e-13

max iteration

1000

reuse preconditioner

false

Table 57. SNES configuration

relative tolerance

1e-8

steps tolerance

1e-8

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

Table 58. KSP in SNES configuration

relative tolerance

1e-5

max iteration

1000

max iteration with reuse

CFD1/CFD2 : 100 | CFD3 : 1000

reuse preconditioner

false

reuse preconditioner rebuild at first Newton step

false

Table 59. Preconditioner configuration

type

lu

package

mumps

30.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 feelpp-toolboxes docker reads

Command line to execute CFD1 testcase
$ mpirun -np 4 /usr/local/bin/feelpp_toolbox_fluid_2d --config-file cfd1.cfg

The result files are then stored by default in

Results Directory
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

30.7. Results

Here are results from the different cases studied in this benchmark.

30.7.1. CFD1

Table 60. Results for 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]

30.7.2. CFD2

Table 61. Results for 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]

30.7.3. CFD3

As CFD3 is time-dependent ( 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}(max-min)

  • 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}

Table 62. Results for CFD3
\(\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 (\(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].

TurekHron CFD3 results
Figure 8. Lift and drag forces

30.8. Geometrical Order

Add a section on geometrical order.

30.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.

30.10. Bibliography

References for this benchmark
  • [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure 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.

Unresolved directive in modeling.adoc - include::../man/04-learning/CFD/MultiFluid/README.adoc[]

31. 2D Drops Benchmark

This benchmark has been proposed and realised by Hysing. It allows us to verify our level set code, our Navier-Stokes solver and how they couple together.

Computer codes, used for the acquisition of results, are from Vincent Doyeux, with the use of Chabannes's Navier-Stokes code.

31.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 Navier-Stokes 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.

2D Bubble Geo
Figure 9. Geometry used in 2D Bubble Benchmark

Durig this benchmark, we will study two different cases : the first one with a ellipsoidal bubble and the second one with a squirted bubble.

31.1.1. Boundary conditions

  • On the lateral walls, we imposed slip conditions

\[\begin{eqnarray} \boldsymbol{u}\cdot\boldsymbol{n}&=&0 \\ t\cdot(\nabla\boldsymbol{u}+^t\nabla\boldsymbol{u})\cdot \boldsymbol{n}&=&0 \end{eqnarray}\]
  • On the horizontal walls, no slip conditions are imposed : \(\boldsymbol{u}=0 \)

31.1.2. 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\)

31.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 63. Fixed and Variable Input Parameters
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\)

31.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

\[\boldsymbol{X}_c = \dfrac{ \displaystyle \int_{\Omega_2} \boldsymbol{x}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{ \displaystyle \int_\Omega \boldsymbol{x} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }\]
\[\boldsymbol{U}_c = \dfrac{\displaystyle \int_{\Omega_2} \boldsymbol{u}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{\displaystyle \int_\Omega \boldsymbol{u} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }\]
\[c = \dfrac{\left(4 \pi \displaystyle \int_{\Omega_2} 1 \right)^{\frac{1}{2}}}{ \displaystyle \int_{\Gamma} 1} = \dfrac{ \left(4 \pi \displaystyle \int_{\Omega} (1 - H_\varepsilon(\phi)) \right) ^{\frac{1}{2}}}{ \displaystyle \int_{\Omega} \delta_\varepsilon(\phi)}\]

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.

31.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

31.6. Results

31.6.2. Test 2

We describe the different quantitative results for the two studied cases.

Table 64. Results comparison between benchmark values and our results for the ellipsoidal case

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

Table 65. Results comparison between benchmark values and our results for the squirted case

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

31.7. Bibliography

References for this benchmark
  • [Hysing] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska, Quantitative benchmark computations of two-dimensional 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 multi-fluides, Application aux écoulements sanguins, PhD thesis, Université de Grenoble, 2014.

32. 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}.

32.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:

\[\begin{array}[lll] \rho\rho(\phi(\mathbf{x}) ) \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) + \nabla p - \nabla \cdot \left( \nu(\phi(\mathbf{x})) (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) \right) &=& \rho ( \phi(\mathbf{x}) ) \mathbf{g}, \\ \nabla \cdot \mathbf{u} &=& 0, \\ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi &=& 0, \end{array}\]

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, no-slip 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 two-fluids 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}.

Table 66. Numerical parameters taken for the benchmarks.

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}\).

32.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: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation. summarize for the test 1 the different simulation properties and the table Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom give the carachteristics of each mesh.

Table 67. Numerical parameters used for the test 1 simulations: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation.

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

Table 68. Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom

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 Navier-Stokes 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.

32.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.

32.4. Bibliography

.

References for this benchmark
  • [cottet]

  • [Feelpp] C. Prud’homme et al.

  • [osher] Osher1988, book_Sethian, book_Osher

  • [Franca1992] Franca 1992

33. Computational Solid Mechanics

Computational solid mechanics ( or CSM ) is a part of mechanics that studies solid deformations or motions under applied forces or other parameters.

33.1. Second Newton’s law

The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows

\[\rho^*_s \frac{\partial^2\boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot (\boldsymbol{F}_s\boldsymbol{\Sigma}_s) = \boldsymbol{f}^t_s\]

It’s define here into a Lagrangian frame.

33.1.1. Variables, symbols and units

Notation

Quantity

Unit

\(\rho_s^*\)

strucure density

\(kg/m^3\)

\(\boldsymbol{\eta}_s\)

displacement

\(m\)

\(\boldsymbol{F}_s\)

deformation gradient

\(\boldsymbol{\Sigma}_s\)

second Piola-Kirchhoff tensor

\(N/m^2\)

\(f_s^t\)

body force

\(N/m^2\)

33.2. Lamé coefficients

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

\[\lambda_s = \frac{E_s\nu_s}{(1+\nu_s)(1-2\nu_s)} \hspace{0.5 cm} , \hspace{0.5 cm} \mu_s = \frac{E_s}{2(1+\nu_s)}\]

33.3. 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.

Reduced Model Geometry
Figure 1 : Geometry of the reduce model

The mathematical problem associated to this reduce model can be describe as

\[ \rho^*_s h \frac{\partial^2 \eta_s}{\partial t^2} - k G_s h \frac{\partial^2 \eta_s}{\partial x^2} + \frac{E_s h}{1-\nu_s^2} \frac{\eta_s}{R_0^2} - \gamma_v \frac{\partial^3 \eta}{\partial x^2 \partial t} = f_s.\]

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.

34. Turek Hron CSM Benchmark

In order to validate our fluid-structure interaction solver, we realize here a benchmark on the deformation of an elastic structure, initially proposed by Turek and Hron.

Computer codes, used for the acquisition of results, are from Vincent Chabannes.

This benchmark is linked to the Turek Hron CFD and Turek Hron FSI benchmarks.

34.1. Problem Description

We consider a solid structure, composed of a hyperelastic bar, bound to one of his extremity \(\Gamma_F^*\) to a rigid stationary circular structure. We denote \(\Gamma_{L}^*=\partial\Omega_s^* \backslash \Gamma_F^*\) the other boundaries. The geometry can be represented as follows

CSM Geometry
Figure 1 : Geometry of the Turek & Hron CSM benchmark.

\(\Omega_s^*\) represent the initial domain, before any deformations. We denote \(\Omega^t_s\) the domain obtained during the deformations application, at the time t.

By the Newton’s second law, we can describe the fundamental equation of the solid mechanic in a Lagrangian frame. Furthermore, we will suppose that the hyperelastic material follows a compressible Saint-Venant-Kirchhoff model. The Lamé coefficients, used in this model, are deducing from the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material.

We will observe, during this simulation, the displacement of \(A\), on the \(x\) and \(y\) axis, when the elastic structure is subjected to its own weight, and compare our results to the reference ones given Turek and Hron.

In the first two cases, CSM1 and CSM2, we want to determine the steady state condition. To find it, a quasi-static algorithm is used, which increase at each step the gravity parameter.

For the third cases, CSM3, we realise a transient simulation, where we will observe the comportment of \(A\) subjected to its weight during a given time span.

34.1.1. Boundary conditions

We set

  • on \(\Gamma_{F}^*\), a condition that imposes this boundary to be fixed : \(\boldsymbol{\eta}_s=0\)

  • on \(\Gamma_{L}^*\), a condition that lets these boundaries be free from constraints : \((\boldsymbol{F}_s\boldsymbol{\Sigma}_s)\boldsymbol{ n }^*_s=\boldsymbol{0}\)

where \(\boldsymbol{n}_s^*\) is the outer unit normal vector from \(\partial \Omega_s^*\).

34.1.2. Initial conditions

The source term \(\boldsymbol{f}_s\), chosen in order to set in motion the elastic structure, is define by

\[\boldsymbol{f}_s=(0,-\rho^*_s g)\]

where \(g\) is a gravitational constant.

The reference document ( Turek and Hron ) don’t specify the time interval used to obtain their results. In our case, it’s chosen as \(\lbrack0,10\rbrack\).

34.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 69. Fixed and Variable Input Parameters
Name Description Nominal Value Units

\(g\)

gravitational constant

2

\(m / s^2\)

\(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\)

\(A\)

control point coordinates

\((0.2,0.2)\)

\(m\)

\(B\)

point coordinates

\((0.15,0.2)\)

\(m\)

\(E_s\)

Young’s modulus

\(1.4\times 10^6\)

\(kg / ms^2\)

\(\nu_s\)

Poisson’s ratio

\(0.4\)

dimensionless

\(\rho^*_s\)

density

\(1000\)

\(kg/ m^3\)

As for solvers we used, Newton’s method is chosen for the non-linear part and a direct method based on LU decomposition is selected for the linear part.

34.3. Outputs

As described before, we have

\[\rho^*_s \frac{\partial^2\boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot (\boldsymbol{F}_s\boldsymbol{\Sigma}_s) = \boldsymbol{f}^t_s\]

We search the displacement \(\boldsymbol{\eta}_s\), on \(\Omega_s^*\), which will satisfy this equation.

In particular, the displacement of the point \(A\) is the one that interests us.

34.4. Discretization

To realize these tests, we made the choice to used Finite Elements Method, with Lagrangian elements of order \(N\) to discretize space.

Newmark-beta method, presented into Chabannes papers, is the one we used for the time discretization. We used this method with \(\gamma=0.5\) and \(\beta=0.25\).

34.4.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

Table 70. KSP configuration

type

gmres

relative tolerance

1e-13

max iteration

1000

reuse preconditioner

true

Table 71. SNES configuration

relative tolerance

1e-8

steps tolerance

1e-8

max iteration

500

max iteration with reuse

10

reuse jacobian

false

reuse jacobian rebuild at first Newton step

true

Table 72. KSP in SNES configuration

relative tolerance

1e-5

max iteration

500

reuse preconditioner

CSM1/CSM2 : false | CSM3 : true

reuse preconditioner rebuild at first Newton step

true

Table 73. Preconditioner configuration

type

lu

package

mumps

34.5. Implementation

To realize the acquisition of the benchmark results, code files contained and using the Feel++ library will be used. Here is a quick look to the different location of them.

First at all, the main code can be found in

    feelpp/applications/models/solid

The configuration file for the CSM3 case, the only one we work on, is located at

    feelpp/applications/models/solid/TurekHron

The result files are then stored by default in

    feel/applications/models/solid/TurekHron/csm3/"OrderDisp""Geometric_order"/"processor_used"

Like that, for the CSM3 case executed on 8 processors, with a \(P_1\) displacement approximation space and a geometric order of 1, the path is

    feel/applications/models/solid/TurekHron/csm3/P1G1/np_8

At least, to retrieve results that interested us for the benchmark and to generate graphs, we use a Python script located at

    feelpp-benchmarking-book/CFD/Turek-Hron/postprocess_cfd.py

34.6. Results

34.6.1. CSM1

\(N_{elt}\)

\(N_{dof}\)

\(x\) displacement \(\lbrack\times 10^{-3}\rbrack\)

\(y\) displacement \(\lbrack\times 10^{-3}\rbrack\)

Reference TurekHron

-7.187

-66.10

1061

4620 (\(P_2\))

-7.039

-65.32

4199

17540 (\(P_2\))

-7.047

-65.37

16495

67464 (\(P_2\))

-7.048

-65.37

1061

10112 (\(P_3\))

-7.046

-65.36

1906

17900 (\(P_3\))

-7.049

-65.37

1061

17726 (\(P_4\))

-7.048

-65.37

All the files used for this case can be found in this rep [ geo file, config file, json file ]

34.6.2. CSM2

\(N_{elt}\)

\(N_{dof}\)

\(x\) displacement \(\lbrack\times 10^{-3}\rbrack\)

\(y\) displacement \(\lbrack\times 10^{-3}\rbrack\)

Reference TurekHron

-0.4690

-16.97

1061

4620 (\(P_2\))

-0.459

-16.77

4201

17548 (\(P_2\))

-0.459

-16.77

16495

67464 (\(P_2\))

-0.459

-16.78

1061

10112 (\(P_3\))

-0.4594

-16.78

16475

150500 (\(P_3\))

-0.460

-16.78

1061

17726 (\(P_4\))

-0.460

-16.78

All the files used for this case can be found in this rep [geo file, config file, json file].

34.6.3. CSM3

The results of the CSM3 benchmark are detailed below.

Table 74. Results for CSM3

\(\Delta t\)

\(N_{elt}\)

\(N_{dof}\)

\(x\) displacement \(\lbrack\times 10^{-3}\rbrack\)

\(y\) displacement \(\lbrack\times 10^{-3}\rbrack\)

/

Reference TurekHron

−14.305 ± 14.305 [1.0995]

−63.607 ± 65.160 [1.0995]

0.02

4199

17536(\(P_2\))

-14.585 ± 14.590 [1.0953]

-63.981 ± 65.521 [1.0930]

4199

38900(\(P_3\))

-14.589 ± 14.594 [1.0953]

-63.998 ± 65.522 [1.0930]

1043

17536(\(P_4\))

-14.591 ± 14.596 [1.0953]

-64.009 ± 65.521 [1.0930]

4199

68662(\(P_4\))

-14.590 ± 14.595 [1.0953]

-64.003 ± 65.522 [1.0930]

0.01

4199

17536(\(P_2\))

-14.636 ± 14.640 [1.0969]

-63.937 ± 65.761 [1.0945]

4199

38900(\(P_3\))

-14.642 ± 14.646 [1.0969]

-63.949 ± 65.771 [1.0945]

1043

17536(\(P_4\))

-14.645 ± 14.649 [1.0961]

-63.955 ± 65.778 [1.0945]

4199

68662(\(P_4\))

-14.627 ± 14.629 [1.0947]

-63.916 ± 65.739 [1.0947]

0.005

4199

17536(\(P_2\))

-14.645 ± 14.645 [1.0966]

-64.083 ± 65.521 [1.0951]

4199

38900(\(P_3\))

-14.649 ± 14.650 [1.0966]

-64.092 ± 65.637 [1.0951]

1043

17536(\(P_4\))

-14.652 ± 14.653 [1.0966]

-64.099 ± 65.645 [1.0951]

4199

68662(\(P_4\))

-14.650 ± 14.651 [1.0966]

-64.095 ± 65.640 [1.0951]

fullviewCSM

\text{Figure 2: x and y displacements}

All the files used for this case can be found in this rep [ geo file, config file, json file ]

34.6.4. Conclusion

To obtain these data, we used several different mesh refinements and different polynomial approximations for the displacement on the time interval \(\lbrack 0,10 \rbrack\).

Our results are pretty similar to those from Turek and Hron, despite a small gap. This gap can be caused by the difference between our time interval and the one used for the reference acquisitions.

34.7. Bibliography

References for this benchmark
  • [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure 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é de Grenoble, 2013.

35. Solenoïd Benchmark

An analytical solution of the elasticity equation exists in an infinitely long solenoid. Such a geometry has the advantage to be axisymetrical, we can first expresse the solution with cylindrical coordinates.

35.1. Definition

35.1.1. Equilibrium equation in cylindrical coordinates

The analytical solution comes from the electromagnetism equations, indeed this case of a soleinoid crossed by an electric current can model the behaviour of a resistive magnet. The volumic forces \(\mathbf{f}\) of the equilibrium equation are consequently dependent on the current density \(\mathbf{J}\) and the induced magnetic field \(\mathbf{B}\), which gives :

\[\nabla\cdot\bar{\bar{\sigma}} + \mathbf{J}\times\mathbf{B} = \mathbf{0}\]

We will use the following notations:

\[\mathbf{u} = \begin{pmatrix} u \\ v \\ w \\ \end{pmatrix} \quad \mathbf{J} = \begin{pmatrix} J_{r} \\ J_{\theta} \\ J_{z} \\ \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} B_{r} \\ B_{\theta} \\ B_{z} \\ \end{pmatrix}, \quad \text{and} \quad \bar{\bar{\sigma}} = \begin{pmatrix} \sigma_{r r} & \sigma_{r \theta} & \sigma_{r z} \\ \sigma_{\theta r} & \sigma_{\theta \theta} & \sigma_{\theta z} \\ \sigma_{z r} & \sigma_{z \theta} & \sigma_{z z} \\ \end{pmatrix}\]

We consider here that the current distribution in the soleinoid is uniform. That means that \(J_{\theta}\) is constant, and \(J_r = J_z = 0\). The volumic forces \(\mathbf{J} \times \mathbf{B}\) becomes :

\[\mathbf{J}\times\mathbf{B} = \begin{pmatrix} J_{\theta} B_z \\ 0 \\ - J_{\theta} B_r \end{pmatrix}\]

The revwriting of the equilibrium equation in cylindrical coordinates gives :

\[\left\{ \begin{array}{l} \displaystyle{ \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\partial \sigma_{rz}}{\partial z} + \frac{1}{r}( \sigma_{rr} - \sigma_{\theta \theta} ) + J_{\theta} B_z = 0 }\\[0.4cm] \displaystyle{ \frac{\partial \sigma_{\theta r}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{\theta \theta}}{\partial \theta} + \frac{\partial \sigma_{\theta z}}{\partial z} + \frac{2}{r} \sigma_{\theta r} = 0 }\\[0.4cm] \displaystyle{ \frac{\partial \sigma_{zr}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{\theta z}}{\partial \theta} + \frac{\partial \sigma_{zz}}{\partial z} + \frac{1}{r} \sigma_{rz} - J_{\theta} B_r = 0 }\\ \end{array} \right.\]

The axisymetrical properties of this geometries means that the displacement is invariant with respect to \(\theta\).
Furthermore, the soleinoid we consider has the particularity to be infinitely long, so there is no displacement along \(z\) axis.
We can consequently get rid of all derivatives \(\frac{\partial \cdot}{\partial \theta}\) and \(\frac{\partial \cdot}{\partial z}\).
Finally, we shall note that components \(\sigma_{\theta r}\) and \(\sigma_{zr}\) of the stress tensor \(\bar{\bar{\sigma}}\) are expressed from Hooke’s law only from the components \(v\) and \(w\) of the displacement vector \(\mathbf{u}\), which nullify the two last equations.

Thus, the equilibrium equation is reduced to:

\[-\sigma_{\theta \theta} + \frac{\partial}{\partial r}(r \sigma_{rr}) = -r J_{\theta} B_z\]

We need to express this equation in terms of the displacement \(\mathbf{u}\), This can be done using Hooke’s law which links the stress tensor to the tensor of small deformation by:

\[\bar{\bar{\sigma}}=2\mu\bar{\bar{\varepsilon}} + \lambda Tr(\bar{\bar{\varepsilon}})Id\]

where \(\mu\) and \(\lambda\) are the Lamé coefficients, and the tensor of small deformation is given in cylindrical coordinates by:

\[\bar{\bar{\varepsilon}} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T) = \begin{pmatrix} \frac{\partial u}{\partial r} & \frac{1}{2}\left( \frac{1}{r}\frac{\partial u}{\partial \theta} - \frac{v}{r} + \frac{\partial v}{\partial r} \right) & \frac{1}{2} \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial r} \right) \\ \frac{1}{2}\left( \frac{1}{r}\frac{\partial u}{\partial \theta} - \frac{v}{r} + \frac{\partial v}{\partial r} \right) & \frac{1}{r}\frac{\partial v}{\partial \theta} + \frac{u}{r} & \frac{1}{2} \left( \frac{1}{r}\frac{\partial w}{\partial \theta} + \frac{\partial v}{\partial z} \right) \\ \frac{1}{2} \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial r} \right) & \frac{1}{2} \left( \frac{1}{r}\frac{\partial w}{\partial \theta} + \frac{\partial v}{\partial z} \right) & \frac{\partial w}{\partial z} \end{pmatrix}\]

Then, using the last two definitions and the properties of the solenoid (axisymmetric and infinitely long), we can rewrite the equilibrium equation as:

\[\frac{\partial}{\partial r}\left( r \frac{\partial u}{\partial r} \right) - \frac{u}{r} = - \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}r J_{\theta} B_z\]

35.1.2. Analytical solution

We want to find an analytical solution of the form :

\[u_{cyl}(r) = C_1 r + \frac{C_2}{r} + u_p(r)\]

where \(C_1\) and \(C_2\) are constants, and \(u_p(r)\) a particular solution of the equilibrium equation in cylindrical coordinates.

From Ampére’s theorem and considering that the soleinoid is axisymetrical and infinitely long, we deduce that \(B_r = 0\), and \(B_z\) depends only on \(r\), such that :

\[B_z(r_1) - B_z(r) = \frac{1}{\mu} \int_{r_1}^r J_{\theta} dr\]

with \(r_1\) the internal radius of the soleinoid.

For a uniform distribution of current in the solenoid (\(j_{\theta}\) constant), we deduce that \(B_z\) can be expressed as :

\[B_z(r) = B_z(r_1) - \frac{\Delta B_z}{\alpha - 1} \left( \frac{r}{r_1} - 1 \right)\]

Replacing \(b_z\) with his expression in the equilibrium equation, this gives :

\[ \frac{\partial}{\partial r}\left( r \frac{\partial u}{\partial r} \right) - \frac{u}{r} = - \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}r J_{\theta} \left( B_z(r_1) - \frac{\Delta B_z}{\alpha - 1} \left( \frac{r}{r_1} - 1 \right) \right)\]

where \(r_2\) is the external radius, \(\alpha = \frac{r_2}{r_1}\) and \(\Delta b_z = B_z(r_1) - B_z(r_2)\).

A particular solution \(u_p(r)\) for this equation is given by:

\[u_p(r) = \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}r_1 J_{\theta} \left[ -\frac{r_1}{3}\left( B_z(r_1) + \frac{\Delta B_z}{\alpha - 1} \right) \left( \frac{r}{r_1}\right)^2 + \frac{r_1}{8}\frac{\Delta B_z}{\alpha - 1} \left( \frac{r}{r_1}\right)^3 \right]\]

The constants \(C_1\) and \(C_2\) are set by the boundary conditions, we consider here that there is no surface forces. That gives \(\bar{\bar{\sigma}}\cdot\mathbf{n} = 0\) on internal and external radius, that is \(\sigma_{rr}(r_1)=\sigma_{rr}(r_2)=0\).

Using the definition of \(u_{cyl}\) in

\[\sigma_{rr} = \frac{E}{(1+\nu)(1-2\nu)} \left[ (1-\nu)\frac{\partial u}{\partial r} + \nu \frac{u}{r} \right]\]

we can solve the system to find the constants:

\[\begin{align} C_1 = \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}J_{\theta}r_1 &\left[ \left( B_z(r_1) + \frac{\Delta B_z}{\alpha - 1}\right) \left( \frac{2 - \nu}{3} \right) \left( 1 - \frac{r_2^2}{(r_2^2 - r_1^2)}\left( 1 - \frac{r_2}{r_1} \right) \right) \right. \\ &+ \left. \frac{\Delta B_z}{\alpha - 1}\left( \frac{2\nu - 3}{8} \right) \left( 1 - \frac{r_2^2}{(r_2^2 - r_1^2)} \left( 1 - \left( \frac{r_2}{r_1} \right)^2 \right) \right) \right] \end{align}\]
\[C_2 = \frac{-r_1^3 r_2^2 J_{\theta}}{(r_2^2 - r_1^2)}\frac{(1+\nu)}{E(1-\nu)} \left[ \left( B_z(r_1) + \frac{\Delta B_z}{\alpha - 1}\right) \left( \frac{2 - \nu}{3} \right)\left( 1 - \frac{r_2}{r_1} \right) +\frac{\Delta B_z}{\alpha - 1}\left( \frac{2\nu - 3}{8} \right)\left( 1 - \left( \frac{r_2}{r_1} \right)^2 \right) \right]\]

The final step is to translate this analytical solution \(u_{cyl}(r)\) into cartesian coordinates to obtain the analytical cartesian displacement \(\mathbf{u}_{cart}\):

\[\mathbf{u}_{cart}(x,y)= \begin{pmatrix} cos( \theta )u_{cyl}(\sqrt{x^2 + y^2})\\ sin( \theta )u_{cyl}(\sqrt{x^2 + y^2})\\ 0 \end{pmatrix} = \begin{pmatrix} \frac{x}{\sqrt{x^2 + y^2}}u_{cyl}(\sqrt{x^2 + y^2}) \\ \frac{y}{\sqrt{x^2 + y^2}}u_{cyl}(\sqrt{x^2 + y^2}) \\ 0 \end{pmatrix}\]

35.1.3. Geometry

We use a solenoïd of thickness one with \(r_1=1\) and \(r_2=2\) and with a length sufficiently important (\(l=10\,r_2\)) so that the influence of the top and of the bottom of the geometry, which are supposed not to exist, is close to zero.

35.1.4. Boundary conditions

The boundary conditions taken into account for the analytical solution have to be reproduced for the simulation. That means null pressure forces on internal and external radius, and displacement set to zero (Dirichlet) on the top and on the bottom to keep only the radial component.

We set:

  • \(\mathbf{u} = 0\) on \(\Gamma_{top}\cup\Gamma_{bottom}\)

  • \(\bar{\bar{\sigma}}\cdot \mathbf{n} = 0\) on \(\Gamma_{int}\cup\Gamma_{ext}\)

  • \(\mathbf{f} = \begin{pmatrix} \frac{x}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ \frac{y}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ 0 \end{pmatrix}\) in \(\Omega\)

35.2. Inputs

We use the following parameters:

Table 75. Inputs
Name Value

\(E\)

\(2.1e^6\)

\(\nu\)

\(0.33\)

\(\mathbf{f}\)

\(\begin{pmatrix} \frac{x}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ \frac{y}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ 0 \end{pmatrix}\)

35.3. Output

We compare the radial component of the displacement on the segment \(z=l/2\), \(y=0\) and \(x\in \lbrack 1,2\rbrack \).

35.4. Results

Here are the analytical and the computed \(x\) component of the displacement. This has been obtain with a characteristic size of \(0.1\) and \(646 233\) dofs.

We can see that the errors grows as we approach the external radius. But the max of the error is \(5e^{-4}\) and it converges as the characteristic size decreases.

36. NAFEMS LE1 Benchmarck

This benchmark is extract from the Abaqus Benchmarks Manual.

36.1. Definition

We focus on the LE1 benchmarks in particular.

36.1.1. Geometry

The geometry is given here by :
geoLE10

36.1.2. Boundary conditions

We set:

  • \(u_y = 0\) on DC

  • \(u_x = 0\) on AB

  • \(\bar{\bar{\varepsilon}}\cdot\mathbf{n}=1e^7\) on BC.

36.2. Inputs

We have the following parameters:

Table 76. Inputs
Name Value

\(E\)

\(210\, GPa\)

\(\nu\)

\(0.3\)

\(\rho\)

\(7800\, kg/m^2\)

36.3. Outputs

We want to compare the value of \(\sigma_{yy}\) at the point D. The reference value is \(92.7\, MPa\).

36.4. Results

The value of \(\sigma_{yy}\) at the point D is \(94.09\, MPa\) for \(32 000\) dofs, which is \(1.49%\) higher than the target.

One possibility to get a more accurate output is to use a mixed formulation, where the stress tensor would also be an unknown.

37. NAFEMS LE10 Benchmarck

This benchmark is extract from the Abaqus Benchmarks Manual.

37.1. Definition

We focus on the LE10 benchmarks in particular.

37.1.1. Geometry

The geometry is given here by :

geoLE10

where the thickness is \(0.6\).

In addition, we define the point E which is the midpoint of CC' and E' the midpoint of BB'.

37.1.2. Boundary conditions

We set:

  • \(u_y = 0\) on DCD’C'

  • \(u_x = 0\) on ABA’B'

  • \(u_x = u_y = 0\) on BCB’C'

  • \(\bar{\bar{\varepsilon}}\cdot\mathbf{n}=-1e^6\) on the top surface.

37.2. Inputs

We have the following parameters:

Table 77. Inputs
Name Value

\(E\)

\(210\, GPa\)

\(\nu\)

\(0.3\)

\(\rho\)

\(7800\, kg/m^2\)

37.3. Outputs

We want to compare the value of \(\sigma_{yy}\) at the point D. The reference value is \(5.38\, MPa\).

37.4. Results

The value of \(\sigma_{yy}\) at the point D is \(5.53\, MPa\) for \(300 000\) dofs, which is \(2.78%\) higher than the target.

One possibility to get a more accurate output is to use a mixed formulation, where the stress tensor would also be an unknown.

38. Thermal Building Environment

39. ISO 10211:2007 Thermal bridges in building construction

39.1. Introduction

ISO 10211:2007 sets out the specifications for a three-dimensional and a two-dimensional geometrical model of a thermal bridge for the numerical calculation of:

  1. heat flows, in order to assess the overall heat loss from a building or part of it;

  2. minimum surface temperatures, in order to assess the risk of surface condensation.

These specifications include the geometrical boundaries and subdivisions of the model, the thermal boundary conditions, and the thermal values and relationships to be used.

ISO 10211:2007 is based upon the following assumptions:

  1. all physical properties are independent of temperature;

  2. there are no heat sources within the building element.

ISO 10211:2007 can also be used for the derivation of linear and point thermal transmittances and of surface temperature factors. More information here.

39.2. Implementation

Only the 2D specifications have been implemented.

40. Running the testcase

$ mpirun -np 4 /usr/local/bin/feelpp_thermodyn_2d --config-file thermo2dCase2.cfg

41. Thermal Building Environment

42. ISO 10211:2007 Thermal bridges in building construction

42.1. Introduction

ISO 10211:2007 sets out the specifications for a three-dimensional and a two-dimensional geometrical model of a thermal bridge for the numerical calculation of:

  1. heat flows, in order to assess the overall heat loss from a building or part of it;

  2. minimum surface temperatures, in order to assess the risk of surface condensation.

These specifications include the geometrical boundaries and subdivisions of the model, the thermal boundary conditions, and the thermal values and relationships to be used.

ISO 10211:2007 is based upon the following assumptions:

  1. all physical properties are independent of temperature;

  2. there are no heat sources within the building element.

ISO 10211:2007 can also be used for the derivation of linear and point thermal transmittances and of surface temperature factors. More information here.

42.2. Implementation

Only the 2D specifications have been implemented.

43. Running the testcase

$ mpirun -np 4 /usr/local/bin/feelpp_thermodyn_2d --config-file thermo2dCase2.cfg

44. ISO 10211:2007 Thermal bridges in building construction

44.1. Introduction

ISO 10211:2007 sets out the specifications for a three-dimensional and a two-dimensional geometrical model of a thermal bridge for the numerical calculation of:

  1. heat flows, in order to assess the overall heat loss from a building or part of it;

  2. minimum surface temperatures, in order to assess the risk of surface condensation.

These specifications include the geometrical boundaries and subdivisions of the model, the thermal boundary conditions, and the thermal values and relationships to be used.

ISO 10211:2007 is based upon the following assumptions:

  1. all physical properties are independent of temperature;

  2. there are no heat sources within the building element.

ISO 10211:2007 can also be used for the derivation of linear and point thermal transmittances and of surface temperature factors. More information here.

44.2. Implementation

Only the 2D specifications have been implemented.

45. Running the testcase

$ mpirun -np 4 /usr/local/bin/feelpp_thermodyn_2d --config-file thermo2dCase2.cfg

46. Heat and Fluid Toolbox Examples

In the toolbox we couple fluid and heat transfer. We model both force and natural (or free) convection.

  1. link:NaturalConvection/README.adoc

47. Natural Convection

The examples in this directory deal with processes due solely to natural convection.

48. Natural Convection in a Cavity

This is a standard benchmark with many data available.

48.1. Description

The goal of this project is to simulate the fluid flow under natural convection: the heated fluid circulates towards the low temperature under the action of density and gravity differences. The phenomenon is important, in the sense it models evacuation of heat, generated by friction forces for example, with a cooling fluid.

We shall put in place a simple convection problem in order to study the phenomenon without having to handle the difficulties of more complex domaines. We describe then some necessary transformations to the equations, then we define quantities of interest to be able to compare the simulations with different parameter values.

48.2. Geometry

Natural Convection Cavity
Figure 1 : Geometry of natural convection benchmark.

To study the convection, we use a model problem: it consists in a rectangular tank of height 1 and width W, in which the fluid is enclosed, see figure Geometry of natural convection benchmark.. We wish to know the fluid velocity \mathbf{u}, the fluid pressure p and fluid temperature \theta.

We introduce the adimensionalized Navier-Stokes and heat equations parametrized by the Grashof and Prandtl numbers. These parameters allow to describe the various regimes of the fluid flow and heat transfer in the tank when varying them.

The adimensionalized steady incompressible Navier-Stokes equations reads:

\begin{split} \mathbf{u}\cdot\nabla \mathbf{u} + \nabla p - \frac{1}{\sqrt{\text{Gr}}} \Delta \mathbf{u} &= \theta \mathbf{e}_2 \\ \nabla \cdot \mathbf{u} &= 0\ \text{sur}\ \Omega\\ \mathbf{u} &= \mathbf{0}\ \text{sur}\ \partial \Omega \end{split}

where \mathrm{Gr} is the Grashof number, \mathbf{u} the adimensionalized velocity and p adimensionalized pressure and \theta the adimensionalized temperature. The temperature is in fact the difference between the temperature in the tank and the temperature T_0 on boundary \Gamma_1.

The heat equation reads:

\begin{split} \mathbf{u} \cdot \nabla \theta -\frac{1}{\sqrt{\text{Gr}}{\mathrm{Pr}}} \Delta \theta &= 0\\ \theta &= 0\ \text{sur}\ \Gamma_1\\ \frac{\partial \theta}{\partial n} &= 0\ \text{sur}\ \Gamma_{2,4}\\ \frac{\partial \theta}{\partial n} &= 1\ \text{sur}\ \Gamma_3 \end{split}

where \mathrm{Pr} is the Prandtl number.

48.3. Influence of parameters

what are the effects of the Grashof and Prandtl numbers ? We remark that both terms with these parameters appear in front of the \Delta parameter, they thus act on the diffusive terms. If we increase the Grashof number or the Prandtl number the coefficients multiplying the diffusive terms decrease, and this the convection, that is to say the transport of the heat via the fluid, becomes dominant. This leads also to a more difficult and complex flows to simulate, see figure [fig:heatns:2]. The influence of the Grashof and Prandtl numbers are different but they generate similar difficulties and flow configurations. Thus we look only here at the influence of the Grashof number which shall vary in [1, 1e7].

flow grashof
Figure 2Velocity norm with respect to Grashof

48.4. Quantities of interest

We would like to compare the results of many simulations with respect to the Grashof defined in the previous section. We introduce two quantities which will allow us to observe the behavior of the flow and heat transfer.

48.4.1. Mean temperature

We consider first the mean temperature on boundary \Gamma_3

T_3 = \int_{\Gamma_3} \theta

This quantity should decrease with increasing Grashof because the fluid flows faster and will transport more heat which will cool down the heated boundary \Gamma_3. We observe this behavior on the figure [fig:heatns:3].

Mean temperature with respect to the Grashof number

48.4.2. Flow rate

Another quantity of interest is the flow rate through the middle of the tank. We define a segment \Gamma_f as being the vertical top semi-segment located at W/2 with height 1/2, see figure [fig:heatns:1]. The flow rate, denoted \mathrm{D}_f, reads \mathrm{D}_f = \int_{\Gamma_f} \mathbf{u} \cdot \mathbf{e}_1

where \mathbf{e}_1=(1,0). Note that the flow rate can be negative or positive depending on the direction in which the fluid flows.

As a function of the Grashof, we shall see a increase in the flow rate. This is true for small Grashof, but starting at 1e3 the flow rate decreases. The fluid is contained in a boundary layer which is becoming smaller as the Grashof increases.

image::debit_grashof.png[Behavior of the flow rate with respect to the Grashof number; h = 0.02, \mathbb{P}_3 for the velocity, \mathbb{P}_2 for the pressure and \mathbb{P}_1 for the temperature.]

48.5. Running the model

$ mpirun -np 4 /usr/local/bin/feelpp_toolbox_fluid_2d --config-file cfd2d.cfg

49. Fluid Structure Interaction

We will interest now to the different interactions a fluid and a structure can have together with specific conditions.

49.1. Fluid structure model

To describe and solve our fluid-structure interaction problem, we need to define a model, which regroup structure model and fluid model parts.

We have then in one hand the fluid equations, and in the other hand the structure equations.

FSI
Figure 1 : FSI case example

The solution of this model are \((\mathcal{A}^t, \boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s)\).

49.2. ALE

Generally, the solid mechanic equations are expressed in a Lagrangian frame, and the fluid part in Eulerian frame. To define and take in account the fluid domain displacement, we use a technique name ALE ( Arbitrary Lagrangian Eulerian ). This allow the flow to follow the fluid-structure interface movements and also permit us to have a different deformation velocity than the fluid one.

Let denote \(\Omega^{t_0}\) the calculation domain, and \(\Omega^t\) the deformed domain at time \(t\). As explain before, we want to conserve the Lagrangian and Eulerian characteristics of each part, and to do this, we introduce \(\mathcal{A}^t\) the ALE map.

This map give us the position of \(x\), a point in the deformed domain at time \(t\) from the position of \(x^*\) in the initial configuration \(\Omega^*\).

ALE
Figure 2 : ALE map

\(\mathcal{A}^t\) is a homeomorphism, i.e. a continuous and bijective application we can define as

\[\begin{eqnarray*} \mathcal{A}^t : \Omega^* &\longrightarrow & \Omega^{t} \\ \mathbf{x}^* &\longmapsto & \mathbf{x} \left(\mathbf{x}^*,t \right) = \mathcal{A}^t \left(\mathbf{x}^*\right) \end{eqnarray*}\]

We denote also \(\forall \mathbf{x}^* \in \Omega^*\), the application :

\[\begin{eqnarray*} \mathbf{x} : \left[t_0,t_f \right] &\longrightarrow & \Omega^t \\ t &\longmapsto & \mathbf{x} \left(\mathbf{x}^*,t \right) \end{eqnarray*}\]

This ALE map can then be retrieve into the fluid-structure model.

50. 2D FSI elastic tube

This test case has originally been realised by [Pena], [Nobile] and [GerbeauVidrascu] only with the free outlet condition.

Computer codes, used for the acquisition of results, are from Vincent [Chabannes]

50.1. Problem Description

We interest here to the case of bidimensional blood flow modelisation. We want to reproduce and observe pressure wave spread into a canal with a fluid-structure interaction model.

Elastic Tube Geometry
Figure 1 : Geometry of two-dimension elastic tube.

The figure above shows us the initial geometry we will work on. The canal is represent by a rectangle with width and height, respectively equal to 6 and 1 cm. The upper and lower walls are mobile and so, can be moved by flow action.

By using the 1D reduced model, named generalized string and explained by [Chabannes], we didn’t need to define the elastic domain ( for the vascular wall ) here. So the structure domain is \(\Omega_s^*=\Gamma_{fsi}^*\)

During this benchmark, we will compare two different cases : the free outlet condition and the Windkessel model. The first one, as its name said, impose a free condition on the fluid at the end of the domain. The second one is used to model more realistically an flow outlet into our case. The chosen time step is \(\Delta t=0.0001\)

50.1.1. Boundary conditions

We set :

  • on \(\Gamma_f^{i,*}\) the pressure wave pulse \[ \boldsymbol{\sigma}_{f} \boldsymbol_f = \left\{ \begin{aligned} & \left(-\frac{2 \cdot 10^4}{2} \left( 1 - \cos \left( \frac{ \pi t} {2.5 \cdot 10^{-3}} \right) \right), 0\right)^T \quad & \text{ if } t < 0.005 \\ & \boldsymbol{0} \quad & \text{ else } \end{aligned} \right. \]

  • on \(\Gamma_f^{o,*}\)

    • Case 1 : free outlet : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f =0\)

    • Case 2 : Windkessel model ( \(P_0\) proximal pressure, see [Chabannes] ) : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f = -P_0\boldsymbol{n}_f\)

  • on \(\Gamma_f^{i,*} \cup \Gamma_f^{o,*}\) a null displacement : \(\boldsymbol{\eta}_f=0\)

  • on \(\Gamma^*_{fsi}\) : \(\eta_s=0\)

  • We add also the specific coupling conditions, obtained from the axisymmetric reduced model, on \(\Gamma^*_{fsi}\)

\[\dot{\eta}_s \boldsymbol{e}_r - \boldsymbol{u}_f = \boldsymbol{0}\]
\[f_s + \left(J_{\mathcal{A}_f^t} \boldsymbol{F}_{\mathcal{A}_f^t}^{-T} \hat{\boldsymbol{\sigma}}_f \boldsymbol{n}^*_f\right) \cdot \boldsymbol{e}_r = 0\]
\[\boldsymbol{\varphi}_s^t - \mathcal{A}_f^t = \boldsymbol{0}\]

50.2. Inputs

Table 78. Fixed and Variable Input Parameters

Name

Description

Nominal Value

Units

\(E_s\)

Young’s modulus

\(0.75\)

\(dynes.cm^{-2}\)

\(\nu_s\)

Poisson’s ratio

\(0.5\)

dimensionless

\(h\)

walls thickness

0.1

\(cm\)

\(\rho_s\)

structure density

\(1.1\)

\(g.cm^{-3}\)

\(R_0\)

tube radius

\(0.5\)

\(cm\)

\(G_s\)

shear modulus

\(10^5\)

\(Pa\)

\(k\)

Timoshenko’s correction factor

\(2.5\)

dimensionless

\(\gamma_v\)

viscoelasticity parameter

\(0.01\)

dimensionless

\(\mu_f\)

viscosity

\(0.003\)

\(poise\)

\(\rho_f\)

density

\(1\)

\(g.cm^{-3}\)

\(R_p\)

proximal resistance

\(400\)

\(R_d\)

distal resistance

\(6.2 \times 10^3\)

\(C_d\)

capacitance

\(2.72 \times 10^{-4}\)

50.3. Outputs

After solving the fluid struture model, we obtain \((\mathcal{A}^t, \boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s)\)

with \(\mathcal{A}^t\) the ALE map, \(\boldsymbol{u}_f\) the fluid velocity, \(p_f\) the fluid pressure and \(\boldsymbol\eta_s\) the structure displacement

50.4. Discretization

\(\mathcal{F}\) is the set of all mesh faces, we denote \(\mathcal{F}_{stab}\) the face we stabilize

\[\mathcal{F}_{stab} = \mathcal{F} \bigcap \left( \left\{ (x,y) \in \mathrm{R}^2: (x - 0.3) \leqslant 0 \right\} \cup \left\{ (x,y) \in \mathrm{^2: (x - 5.7) } \geqslant 0 \right\} \right)\]

In fact, after a first attempt, numerical instabilities can be observed at the fluid inlet. These instabilities, caused by pressure wave, and especially by the Neumann condition, make our fluid-structure solver diverge.

To correct them, we choose to add a stabilization term, obtain from the stabilized CIP formulation ( see [Chabannes], Chapter 6 ).

As this stabilization bring an important cost with it, by increasing the number of non-null term into the problem matrix, we only apply it at the fluid entrance, where the instabilities are located.

Now we present the different situations we worked on.

Config

Fluid

Structure

\(N_{elt}\)

\(N_{geo}\)

\(N_{dof}\)

\(N_{elt}\)

\(N_{geo}\)

\(N_{dof}\)

\((1)\)

\(342\)

\(3~(P4P3)\)

\(7377\)

\(58\)

\(1\)

\(176~(P3)\)

\((2)\)

\(342\)

\(4~(P5P4)\)

\(11751\)

\(58\)

\(1\)

\(234~(P4)\)

For the fluid time discretization, BDF, at order \(2\), is the method we use.

And Newmark-beta method is the one we choose for the structure time discretization, with parameters \(\gamma=0.5\) and \(\beta=0.25\).

These methods can be retrieved in [Chabannes] papers.

50.4.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

KSP

case

fluid

solid

type

gmres

relative tolerance

\(1e-13\)

max iteration

\(30\)

\(10\)

reuse preconditioner

true

false

SNES

case

fluid

solid

relative tolerance

\(1e-8\)

steps tolerance

\(1e-8\)

max iteration

\(50\)

max iteration with reuse

\(50\)

reuse jacobian

false

reuse jacobian rebuild at first Newton step

false

true

KSP in SNES

case

fluid

solid

relative tolerance

\(1e-5\)

max iteration

\(1000\)

max iteration with reuse

\(1000\)

reuse preconditioner

true

false

reuse preconditioner rebuild at first Newton step

false

PC

case

fluid

solid

type

LU

package

mumps

FSI

solver method

fix point

tolerance

\(1e-6\)

max iterations

\(1\)

50.5. Implementation

To realize the acquisition of the benchmark results, code files contained and using the Feel++ library will be used. Here is a quick look to the different location of them.

Let’s start with the main code, that can be retrieve in

    feelpp/applications/models/fsi

The configuration file associated to this test is named wavepressure2d.cfg and is located at

    feelpp/applications/models/fsi/wavepressure2d

The result files are then stored by default in

    applications/models/fsi/wavepressure2d/P2P1G1-P1G1/np_1

50.6. Results

The two following pictures have their pressure and velocity magnitude amplify by 5.

Elastic Tube Free outlet
Figure 2 : Results with free outlet conditon
Elastic Tube Windkessel
Figure 3 : Results with the Windkessel model
Inflow Outflow
Figure 4 : Evolution of the inflow and the outflow
Displacement Magnitude
Figure 5 : Maximum displacement magnitude

To draw the next two figures, we define 60 sections \(\{x_i\}_{i=0}^{60}\) with \(x_i=0.1i\).

Average Pressure
Figure 5 : Average pressure with the free outlet and the Windkessel model
Flow Rate
Figure 7 : Flow rate with the free outlet and the Windkessel model
FSI 1
Figure 8 : Implicit and semi-implicit FSI coupling comparison
FSI 2
Figure 9 : Implicit and semi-implicit FSI coupling comparison

All the files used for this case can be found in this rep [geo file, config file, fluid json file, solid json file].

50.6.1. Conclusion

Let’s begin with results with the free outlet condition ( see figure 2 ). These pictures show us how the pressure wave progresses into the tube. We can denote an increase of the fluid velocity at the end of the tube. Also, the wave eases at the same place. For the simulation with the Windkessel model, we observe a similar comportment at the beginning ( see figure 3 ). However, the outlet is more realistic than before. In fact, the pressure seems to propagate more naturally with this model. + In the two cases, the velocity field is disturbed at the fluid-structure interface. A mesh refinement around this region increases the quality. However, this is not crucial for the blood flow simulation.

Now we can interest us to the quantitative results.

The inflow and outflow evolution figure ( see figure 4 ) shows us similarities for the two tests at the inlet. At the outlet, in contrast, the flow increases for the free outlet condition. In fact, when the pressure wave arrived at the outlet of the tube, it is reflected to the other way. In the same way, when the reflected wave arrived at the inlet, it is reflected again. The Windkessel model reduce significantly this phenomenon. Some residues stay due to 0D coupling model and structure fixation.

We also have calculate the maximum displacement magnitude for the two model ( see figure 5 ). The same phenomenons explained ahead are retrieve here. We denote that, for the free outlet, the structure undergoes movements during the test time, caused by the wave reflection. The Windkessel model reduces these perturbations thanks to the 0D model.

The average pressure and the fluid flow ( see figure 6 and 7 ) show us the same non-physiological phenomenons as before. The results we obtain are in accordance with the ones proposed by [Nobile].

To end this benchmark, we will compare the two resolution algorithms used with the fluid-structure model : the implicit and the semi-implicit ones. The second configuration with Windkessel model is used for the measures.

We have then the fluid flow and the displacement magnitude ( figure 8 ) curves, which superimposed on each other. So the accuracy obtained by the semi-implicit method seems good here.

The performances of the two algorithms ( figure 9 ) are expressed from number of iterations and CPU time at each step time. The semi-implicit method is a bit ahead of the implicit one on number of iterations. However, the CPU time is smaller for 2 or 3 time, due to optimization in this method. First an unique ALE map estimation is need. Furthermore, linear terms of the Jacobian matrix, residuals terms and dependent part of the ALE map can be stored and reused at each iteration.

50.7. Bibliography

References for this benchmark
  • [Pena] G. Pena, Spectral element approximation of the incompressible Navier-Stokes equations evolving in a moving domain and applications, École Polytechnique Fédérale de Lausanne, November 2009.

  • [Nobile] F. Nobile, Numerical approximation of fluid-structure interaction problems with application to haemodynamics, École Polytechnique Fédérale de Lausanne, Switzerland, 2001.

  • [GerbeauVidrascu] J.F. Gerbeau, M. Vidrascu, A quasi-newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows, 2003.

  • [Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Université de Grenoble, 2013.

51. 3D FSI elastic tube

Computer codes, used for the acquisition of results, are from Vincent [Chabannes].

51.1. Problem Description

As in the 2D case, the blood flow modelisation, by observing a pressure wave progression into a vessel, is the subjet of this benchmark. But this time, instead of a two-dimensional model, we use a three-dimensional model, with a cylinder

Elastic Tube Geometry
Figure 1 : Geometry of three-dimensional elastic tube.

This represents the domains into the initial condition, with \(\Omega_f\) and \(\Omega_s\) respectively the fluid and the solid domain. The cylinder radius equals to \(r+\epsilon\), where \(r\) is the radius of the fluid domain and \(\epsilon\) the part of the solid domain.

Elastic Tube Geometry
Figure 2 : Sections of the three-dimensional elastic tube.

\(\Gamma^*_{fsi}\) is the interface between the fluid and solid domains, whereas \(\Gamma^{e,*}_s\) is the interface between the solid domain and the exterior. \(\Gamma_f^{i,*}\) and \(\Gamma_f^{o,*}\) are respectively the inflow and the outflow of the fluid domain. Likewise, \(\Gamma_s^{i,*}\) and \(\Gamma_s^{o,*}\) are the extremities of the solid domain.

During this benchmark, we will study two different cases, named BC-1 and BC-2, that differ from boundary conditions. BC-2 are conditions imposed to be more physiological than the ones from BC-1. So we waiting for more realistics based results from BC-2.

51.1.1. Boundary conditions

  • on \(\Gamma_f^{i,*}\) the pressure wave pulse \[ \boldsymbol{\sigma}_{f} \boldsymbol_f = \left\{ \begin{aligned} & \left(-\frac{1.3332 \cdot 10^4}{2} \left( 1 - \cos \left( \frac{ \pi t} {1.5 \cdot 10^{-3}} \right) \right), 0 \right)^T \quad & \text{ if } t < 0.003 \\ & \boldsymbol{0} \quad & \text{ else } \end{aligned} \right. \]

  • We add the coupling conditions on \(\Gamma^*_{fsi}\)

\[ \frac{\partial \boldsymbol{\eta_{s}} }{\partial t} - \boldsymbol{u}_f \circ \mathcal{A}_{f}^t = \boldsymbol{0} , \quad \text{ on } \Gamma_{fsi}^* \times \left[t_i,t_f \right] \quad \boldsymbol{(1)}\]
\[ \boldsymbol{F}_{s} \boldsymbol{\Sigma}_{s} \boldsymbol{n}^*_s + J_{\mathcal{A}_{f}^t} \hat{\boldsymbol{\sigma}}_f \boldsymbol{F}_{\mathcal{A}_{f}^t}^{-T} \boldsymbol{n}^*_f = \boldsymbol{0} , \quad \text{ on } \Gamma_{fsi}^* \times \left[t_i,t_f \right] \quad \boldsymbol{(2)}\]
\[ \boldsymbol{\varphi}_s^t - \mathcal{A}_{f}^t = \boldsymbol{0} , \quad \text{ on } \Gamma_{fsi}^* \times \left[t_i,t_f \right] \quad \boldsymbol{(3)}\]

Then we have two different cases :

  • Case BC-1

    • on \(\Gamma_f^{o,*}\) : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f =0\)

    • on \(\Gamma_s^{i,*} \cup \Gamma_s^{o,*}\) a null displacement : \(\boldsymbol{\eta}_s=0\)

    • on \(\Gamma^{e,*}_{s}\) : \(\boldsymbol{F}_s\boldsymbol{\Sigma}_s\boldsymbol{n}_s^*=0\)

    • on \(\Gamma_f^{i,*}U \Gamma_f^{o,*}\) : \(\mathcal{A}^t_f=\boldsymbol{\mathrm{x}}^*\)

  • Case BC-2

    • on \(\Gamma_f^{o,*}\) : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f = -P_0\boldsymbol{n}_f\)

    • on \(\Gamma_s^{i,*}\) a null displacement \(\boldsymbol{\eta}_s=0\)

    • on \(\Gamma^{e,*}_{s}\) : \(\boldsymbol{F}_s\boldsymbol{\Sigma}_s\boldsymbol{n}_s^* + \alpha \boldsymbol{\eta}_s=0\)

    • on \(\Gamma^{o,*}_{s}\) : \(\boldsymbol{F}_s\boldsymbol{\Sigma}_s\boldsymbol{n}_s^* =0\)

    • on \(\Gamma_f^{i,*}\) : \(\mathcal{A}^t_f=\boldsymbol{\mathrm{x}}^*\)

    • on \(\Gamma_f^{o,*}\) : \(\nabla \mathcal{A}^t_f \boldsymbol{n}_f^*=\boldsymbol{n}_f^*\)

51.1.2. Initial conditions

The chosen time step is \(\Delta t=0.0001\)

51.2. Inputs

Table 79. Fixed and Variable Input Parameters
Name Description Nominal Value Units

\(E_s\)

Young’s modulus

\(3 \times 10^6 \)

\(dynes.cm^{-2}\)

\(\nu_s\)

Poisson’s ratio

\(0.3\)

dimensionless

\(r\)

fluid tube radius

0.5

\(cm\)

\(\epsilon\)

solid tube radius

0.1

\(cm\)

\(L\)

tube length

5

\(cm\)

\(A\)

A coordinates

(0,0,0)

\(cm\)

\(B\)

B coordinates

(5,0,0)

\(cm\)

\(\mu_f\)

viscosity

\(0.03\)

\(poise\)

\(\rho_f\)

density

\(1\)

\(g.cm^{-3}\)

\(R_p\)

proximal resistance

\(400\)

\(R_d\)

distal resistance

\(6.2 \times 10^3\)

\(C_d\)

capacitance

\(2.72 \times 10^{-4}\)

51.3. Outputs

After solving the fluid struture model, we obtain \((\mathcal{A}^t, \boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s)\)

with \(\mathcal{A}^t\) the ALE map, \(\boldsymbol{u}_f\) the fluid velocity, \(p_f\) the fluid pressure and \(\boldsymbol\eta_s\) the structure displacement.

51.4. Discretization

Here are the different configurations we worked on. The parameter Incomp define if we use the incompressible constraint or not.

Config

Fluid

Structure

\(N_{elt}\)

\(N_{geo}\)

\(N_{dof}\)

\(N_{elt}\)

\(N_{geo}\)

\(N_{dof}\)

Incomp

\((1)\)

\(13625\)

\(1~(P2P1)\)

\(69836\)

\(12961\)

\(1\)

\(12876~(P1)\)

No

\((2)\)

\(13625\)

\(1~(P2P1)\)

\(69836\)

\(12961\)

\(1\)

\(81536~(P1)\)

Yes

\((3)\)

\(1609\)

\(2~(P3P2)\)

\(30744\)

\(3361\)

\(2\)

\(19878~(P2)\)

No

For the structure time discretization, we will use Newmark-beta method, with parameters \(\gamma=0.5\) and \(\beta=0.25\).

And for the fluid time discretization, BDF, at order \(2\), is the method we choose.

These two methods can be found in [Chabannes] papers.

51.5. Implementation

To realize the acquisition of the benchmark results, code files contained and using the Feel++ library will be used. Here is a quick look to the different location of them.

Let’s start with the main code, that can be retrieve in

    feelpp/applications/models/fsi

The configuration file associated to this test is named wavepressure3d.cfg and is located at

    feelpp/applications/models/fsi/wavepressure3d

The result files are then stored by default in

    applications/models/fsi/wavepressure3d/P2P1G1-P1G1/np_1

All the files used for this case can be found in this rep [geo file, config file, fluid json file, solid json file].

51.7. Bibliography

References for this benchmark
  • [Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Université de Grenoble, 2013.

52. Turek Hron FSI Benchmark

In order to validate our fluid-structure interaction solver, a benchmark, initially proposed by [TurekHron], is realized.

Computer codes, used for the acquisition of results, are from Vincent [Chabannes].

Note This benchmark is linked to the Turek Hron CFD and Turek Hron CSM benchmarks.

52.1. Problem Description

In this case, we want to study the flow of a laminar incrompressible flow around an elastic obstacle, fixed to a rigid cylinder, we combine here the study realized on Turek Hron CFD benchmark and Turek Hron CSM benchmark. So a 2D model, representative of this state, is considered.

TurekHron Geometry
Figure 1 : Geometry of the Turek & Hron FSI Benchmark.

We denote \(\Omega_f^*\) the fluid domain, represented by a rectangle of dimension \( [0,2.5] \times [0,0.41] \). This domain is characterized by its density \(\rho_f\) and its dynamic viscosity \(\mu_f\). In this case, the fluid material we used is glycerin.

Furthermore, the model chosen to describe the flow is the incompressible Navier-Stockes model. It is defined by the conservation of momentum equation and the conservation of mass equation. We also have the material constitutive equation to take in account in this domain, as well as the harmonic extension operator for the fluid movement.

On the other side, the structure domain, named \(\Omega_s^*\), represent the hyperelastic bar, bound to the stationary cylinder. As we want to work in a Lagrangian frame, and by Newton’s second law, the fundamental equation of the solid mechanic will be used. For the model, that our hyperelastic material follows, we use the Saint-Venant-Kirchhoff one, define with Lamé coefficients. These coefficients are obtained from the the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material.

All of these are then gather into the fluid-structure coupling system.

52.1.1. Boundary conditions

We set

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

  • on \(\Gamma_{wall}^* \cup \Gamma_{circ}^*\), a homogeneous Dirichlet condition : \(\boldsymbol{u}_f=\boldsymbol{0};\)

  • on \(\Gamma_{out}^*\), a Neumann condition : \(\boldsymbol{\sigma}_f\boldsymbol{n}_f=\boldsymbol{0};\)

  • on \(\Gamma_{fixe}^*\), a condition that imposes this boundary to be fixed : \(\boldsymbol{\eta}_s=0;\)

  • on \(\Gamma_{fsi}^*\) :

\[\frac{\partial \boldsymbol{\eta_{s}} }{\partial t} - \boldsymbol{u}_f \circ \mathcal{A}^t_f = \boldsymbol{0} \quad \left( \text{Velocities continuity}\right)\]
\[ \boldsymbol{F}_{s} \boldsymbol{\Sigma}_{s} \boldsymbol{n}^*_s + J_{\mathcal{A}^t_f} \boldsymbol{F}_{\mathcal{A}^t_f}^{-T} \hat{\boldsymbol{\sigma}}_f \boldsymbol{n}^*_f = \boldsymbol{0} \quad \left( \text{ Constraint continuity}\right)\]
\[ \boldsymbol{\varphi}_s^t - \mathcal{A}^t_f = \boldsymbol{0} \quad \left( \text{Geometric continuity}\right)\]

where \(\boldsymbol{n}_f\) is the outer unit normal vector from \(\partial \Omega_f^*\).

52.1.2. Initial conditions

Fluid

In order to describe the flow inlet by \(\Gamma_{in}\), a parabolic velocity profile is used. It can be expressed by

\[ v_{cst} = 1.5 \bar{U} \frac{4}{0.1681} y \left(0.41-y\right)\]

where \(\bar{U}\) is the mean inflow velocity.

However, to impose a progressive increase of this velocity profile, 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.

Finally, we don’t want a source term, so \(f_f\equiv 0\).

52.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 80. Fixed and Variable Input Parameters
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\)

\(A\)

control point coordinates

\((0.2,0.2)\)

\(m\)

\(B\)

point coordinates

\((0.15,0.2)\)

\(m\)

\(E_s\)

Young’s modulus

\(5.6 \times 10^6\)

\(kg.m^{-1}s^{-2}\)

\(\nu_s\)

Poisson’s ratio

\(0.4\)

dimensionless

\(\rho_s\)

structure density

\(1000\)

\(kg.m^{-3}\)

\(\nu_f\)

kinematic viscosity

\(1\times 10^{-3}\)

\(m^2.s^{-1}\)

\(\mu_f\)

dynamic viscosity

\(1\)

\(kg.m^{-1}.s^{-1}\)

\(\rho_f\)

density

\(1000\)

\(kg.m^{-3}\)

\(f_f\)

source term

0

\(\bar{U}\)

mean inflow velocity

2

\(m.s^{-1}\)

As for solvers we used, Newton’s method is chosen for the non-linear part and a direct method based on LU decomposition is selected for the linear part.

52.3. Outputs

The quantities we observe during this benchmark are on one hand the lift and drag forces ( respectively \(F_L\) and \(F_D\) ), as well as the displacement, on \(x\) and \(y\) axis, of the point A is the second value that interest us here.

They are the solution of the link::../README.adoc[fluid-structure system].

This system also give us the ALE map \(\mathcal{A}_f^t\).

52.4. Discretization

To realize these tests, we made the choice to used \(P_N~-~P_{N-1}\) Taylor-Hood finite elements to discretize the space.

For the fluid time discretization, BDF, at order \(q\), is the method we have chosen.

And finally Newmark-beta method is the one we used for the structure time discretization, with parameters \(\gamma=0.5\) and \(\beta=0.25\).

These methods can be retrieved in [Chabannes] papers.

52.4.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

KSP

case

fluid

solid

type

gmres

relative tolerance

\(1e-13\)

max iteration

\(1000\)

reuse preconditioner

true

SNES

case

fluid

solid

relative tolerance

\(1e-8\)

steps tolerance

\(1e-8\)

max iteration

\(50\)

max iteration with reuse

\(50\)

\(10\)

reuse jacobian

true

false

reuse jacobian rebuild at first Newton step

false

true

SNES

case

fluid

solid

relative tolerance

\(1e-5\)

max iteration

\(1000\)

max iteration with reuse

\(1000\)

\(10\)

reuse preconditioner

true

reuse preconditioner rebuild at first Newton step

true

PC

case

fluid

solid

type

LU

package

mumps

FSI

solver method

fix point with Aitken relaxation

tolerance

\(1e-6\)

max iterations

\(1000\)

initial \(\theta\)

\(0.98\)

minimum \(\theta\)

\(1e-12\)

52.5. Implementation

To realize the acquisition of the benchmark results, code files contained and using the Feel++ library will be used. Here is a quick look to the different location of them.

Let’s start with the main code, that can be retrieve in

    feelpp/applications/models/fsi

The FSI3 configuration file is located at

    feelpp/applications/models/fsi/TurekHron

The result files are then stored by default in

    feel/applications/models/fsi/TurekHron/fsi3/
    <velocity_space><pression_space><Geometric_order>-<OrderDisp><Geometric_order>/np_<processor_used>

For example, for the FSI3 case executed on \(4\) processors, with a \(P_2\) velocity approximation space, a \(P_1\) pressure approximation space, a geometric order of \(1\) for fluid part and a \(P_1\) displacement approximation space and geometric order equals to \(1\) for solid part, the path is

    feel/applications/models/fsi/TurekHron/fsi3/P2P1G1-P1G1/np_4

52.6. Results

First at all, we will discretize the simulation parameters for the different cases studied.

Table 81. Simulation parameters

\(N_{elt}\)

\(N_{dof}\)

\( [P^N_c(\Omega_{f,\delta}]^2 \times P^{N-1}_c(\Omega_{f,\delta}) \times V^{N-1}_{s,\delta}\)

\(\Delta t\)

Turek and Hron

15872

304128

0.00025

(1)

1284

27400

\( [P^4_c(\Omega_{f,(h,3)}]^2 \times P^3_c(\Omega_{f,(h,3)}) \times V^3_{s,(h,3)}\)

0.005

(2)

2117

44834

\( [P^4_c(\Omega_{f,(h,3)}]^2 \times P^3_c(\Omega_{f,(h,3)}) \times V^3_{s,(h,3)}\)

0.005

(3)

4549

95427

\( [P^4_c(\Omega_{f,(h,3)}]^2 \times P^3_c(\Omega_{f,(h,3)}) \times V^3_{s,(h,3)}\)

0.005

(4)

17702

81654

\( [P^2_c(\Omega_{f,(h,1)}]^2 \times P^1_c(\Omega_{f,(h,1)}) \times V^1_{s,(h,1)}\)

0.0005

Then the FSI3 benchmark results are detailed below.

Table 82. Results for FSI3

\(x\) displacement \( [\times 10^{-3}] \)

\(y\) displacement \( [\times 10^{-3}] \)

Drag

Lift

[TurekHron]

-2.69 ± 2.53 [10.9]

1.48 ± 34.38 [5.3]

457.3 ± 22.66 [10.9]

2.22 ± 149.78 [5.3]

[Breuer]

464.5 ± 40.50

6.00 ± 166.00 [5.5]

[TurekHron2]

-2.88 ± 2.72 [10.9]

1.47 ± 34.99 [5.5]

460.5 ± 27.74 [10.9]

2.50 ± 153.91 [5.5]

[MunschBreuer]

-4.54 ± 4.34 [10.1]

1.50 ± 42.50 [5.1]

467.5 ± 39.50 [10.1]

16.2 ± 188.70 [5.1]

[Gallinger]

474.9 ± 28.10

3.90 ± 165.90 [5.5]

[Sandboge]

-2.83 ± 2.78 [10.8]

1.35 ± 34.75 [5.4]

458.5 ± 24.00 [10.8]

2.50 ± 147.50 [5.4]

(1)

-2.86 ± 2.74 [10.9]

1.31 ± 34.71 [5.4]

459.7 ± 29.97 [10.9]

4.46 ± 172.53 [5.4]

(2)

-2.85 ± 2.72 [10.9]

1.35 ± 34.62 [5.4]

459.2 ± 29.62 [10.9]

3.53 ± 172.73 [5.4]

(3)

-2.88 ± 2.75 [10.9]

1.35 ± 34.72 [5.4]

459.3 ± 29.84 [10.9]

3.19 ± 171.20 [5.4]

(4)

-2.90 ± 2.77 [11.0]

1.33 ± 34.90 [5.5]

457.9 ± 31.79 [11.0]

8.93 ± 216.21 [5.5]

All the files used for this case can be found in this rep [geo file, config file, fluid json file, solid json file].

52.6.1. Conclusion

Our first three results are quite similar to given references values. That show us that high order approximation order for space and time give us accurate values, while allow us to use less degree of freedom.

However, the lift force seems to undergo some disturbances, compared to reference results, and it’s more noticeable in our fourth case. This phenomenon is describe by [Beuer], where they’re explaining these disturbances are caused by Aitken dynamic relaxation, used in fluid structure relation for the fixed point algorithm.

In order to correct them, they propose to lower the fixed point tolerance, but this method also lowers calculation performances. An other method to solve this deviation is to use a fixed relaxation parameter \(\theta\). In this case, the optimal \theta seems to be equal to \(0.5\).

52.7. Bibliography

References for this benchmark
  • [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure 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é de Grenoble, 2013.

  • [Breuer] M. Breuer, G. De Nayer, M. Münsch, T. Gallinger, and R. Wüchner, Fluid–structure interaction using a partitioned semi-implicit predictor–corrector coupling scheme for the application of clarge-eddy simulation, Journal of Fluids and Structures, 2012.

  • [TurekHron2] S. Turek, J. Hron, M. Madlik, M. Razzaq, H. Wobker, and JF Acker, Numerical simulation and benchmarking of a monolithic multigrid solver for fluid-structure interaction problems with application to hemodynamics, Fluid Structure Interaction II, pages 193–220, 2010.

  • [MunschBreuer] M. Münsch and M. Breuer, Numerical simulation of fluid–structure interaction using eddy–resolving schemes, Fluid Structure Interaction II, pages 221–253, 2010.

  • [Gallinger] T.G. Gallinger, Effiziente Algorithmen zur partitionierten Lösung stark gekoppelter Probleme der Fluid-Struktur-Wechselwirkung, Shaker, 2010.

  • [Sandboge] R. Sandboge, Fluid-structure interaction with openfsitm and md nastrantm structural solver, Ann Arbor, 1001 :48105, 2010.

53. Lid-driven Cavity Flow benchmark

This test case has originally been proposed by [MokWall].

Computer codes, used for the acquisition of results, are from Vincent [Chabannes].

This benchmark has aso been realized by [Gerbeau], [Vàzquez], [Kuttler] and [Kassiotis].

53.1. Problem Description

We study here an incompressible fluid flowing into a cavity, where its walls are elastic. We use the following geometry to represent it.

LidDriven Geometry
Figure 1 : Geometry of lid-driven cavity flow.

The domain \(\Omega_f^*\) is define by a square \( [0,1]^2 \), \(\Gamma^{i,*}_f\) and \(\Gamma^{o,*}_f\) are respectively the flow entrance and the flow outlet. A constant flow velocity, following the \(x\) axis, will be imposed on \(\Gamma_f^{h,*}\) border, while a null flow velocity will be imposed on \(\Gamma_f^{f,*}\). This last point represent also a non-slip condition for the fluid.

Furthermore, we add a structure domain, at the bottom of the fluid one, named \(\Omega_s^*\). It is fixed by his two vertical sides \(\Gamma_s^{f,*}\), and we denote by \(\Gamma_f^{w,*}\) the border which will interact with \(\Omega_f^*\).

During this test, we will observe the displacement of a point \(A\), positioned at \((0;0.5)\), into the \(y\) direction, and compare our results to ones found in other references.

53.1.1. Boundary conditions

Before enunciate the boundary conditions, we need to describe a oscillatory velocity, following the x axis and dependent of time.

\[ v_{in} = 1-cos\left( \frac{2\pi t}{5} \right)\]

Then we can set

  • on \(\Gamma^{h,*}_f\), an inflow Dirichlet condition : \(\boldsymbol{u}_{f} = ( v_{in}, 0 )\)

  • on \(\Gamma^{i,*}_f\), an inflow Dirichlet condition : \(\boldsymbol{u}_{f} = ( v_{in}(8 y -7) ,0)\)

  • on \(\Gamma^{f,*}_f\), a homogeneous Dirichlet condition : \(\boldsymbol{u}_{f} = \boldsymbol{0}\)

  • on \(\Gamma^{o,*}_f\), a Neumann condition : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f = \boldsymbol{0}\)

  • on \(\Gamma^{f,*}_s\), a condition that imposes this boundary to be fixed : \(\boldsymbol{\eta}_{s} = \boldsymbol{0}\)

  • on \(\Gamma^{e,*}_s\), a condition that lets these boundaries be free from constraints : \(\boldsymbol{F}_{s} \boldsymbol{\Sigma}_s \boldsymbol{n}_s = \boldsymbol{0}\)

To them we also add the fluid-structure coupling conditions on \(\Gamma_{fsi}^*\) :

\[\frac{\partial \boldsymbol{\eta_{s}} }{\partial t} - \boldsymbol{u}_f \circ \mathcal{A}^t_f = \boldsymbol{0} \quad \left( \text{Velocities continuity}\right)\]
\[ \boldsymbol{F}_{s} \boldsymbol{\Sigma}_{s} \boldsymbol{n}^*_s + J_{\mathcal{A}^t_f} \boldsymbol{F}_{\mathcal{A}^t_f}^{-T} \hat{\boldsymbol{\sigma}}_f \boldsymbol{n}^*_f = \boldsymbol{0} \quad \left( \text{ Constraint continuity}\right)\]
\[ \boldsymbol{\varphi}_s^t - \mathcal{A}^t_f = \boldsymbol{0} \quad \left( \text{Geometric continuity}\right)\]

53.1.2. Initial conditions

To realize the simulations, we used a time step \(\Delta t\) equals to \(0.01\) s.

53.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 83. Fixed and Variable Input Parameters

Name

Description

Nominal Value

Units

\(E_s\)

Young’s modulus

\(250\)

\(N.m^{-2}\)

\(\nu_s\)

Poisson’s ratio

\(0\)

dimensionless

\(\rho_s\)

structure density

\(500\)

\(kg.m^{-3}\)

\(\mu_f\)

viscosity

\(1\times 10^{-3}\)

\(m^2.s^{-1}\)

\(\rho_f\)

density

\(1\)

\(kg.m^{-3}\)

53.3. Outputs

\((\boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s, \mathcal{A}_f^t)\), checking the fluid-structure system, are the output of this problem.

53.4. Discretization

To discretize space, we used \(P_N~-~P_{N-1}\) Taylor-Hood finite elements.

For the structure time discretization, Newmark-beta method is the one we used. And for the fluid time discretization, we used BDF, at order \(q\).

53.5. Implementation

All the codes files are into FSI

53.6. Results

We begin with a \(P_2~-~P_1\) approximation for the fluid with a geometry order equals at \(1\), and a fluid-structure stable interface.

Then we retry with a \(P_3~-~P_2\) approximation for the fluid with a geometry order equals at \(2\), and a fluid-structure stable interface.

Finally we launch it with the same conditions as before, but with a deformed interface.

53.6.1. Conclusion

First at all, we can see that the first two tests offer us similar results, despite different orders uses. At contrary, the third result set are better than the others.

The elastic wall thinness, in the stable case, should give an important refinement on the fluid domain, and so a better fluid-structure coupling control. However, the deformed case result are closer to the stable case made measure.

53.7. Bibliography

References for this benchmark
  • [MokWall] DP Mok and WA Wall, Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures, Trends in computational structural mechanics, Barcelona, 2001.

  • [Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Université de Grenoble, 2013.

  • [Gerbeau] J.F. Gerbeau, M. Vidrascu, et al, A quasi-newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows, 2003.

  • [Vazquez] J.G. Valdés Vazquèz et al, Nonlinear analysis of orthotropic membrane and shell structures including fluid-structure interaction, 2007.

  • [KuttlerWall] U. Kuttler and W.A. Wall, Fixed-point fluid–structure interaction solvers with dynamic relaxation, Computational Mechanics, 2008.

  • [Kassiotis] C. Kassiotis, A. Ibrahimbegovic, R. Niekamp, and H.G. Matthies, Nonlinear fluid–structure interaction problem ,part i : implicit partitioned algorithm, nonlinear stability proof and validation examples, Computational Mechanics, 2011.

54. Optimization problems

This section presents some mathematical optimization problems. The following examples source codes are located in "doc/manual/opt/".

55. Laplacian nlopt

Here we propose to solve the following problem, finding the parameters b,c,d solution of the following homogeneous heat equation

\begin{cases} \begin{aligned} -\nabla\cdot(\kappa\nabla u) &=0 \quad \text{on}\;\Omega \\ u &= g \quad \text{on}\;\partial\Omega \end{aligned} \end{cases}

with \kappa\rightarrow\kappa(\kappa_0,\kappa_1,\kappa_2) and where \Omega denotes the whole domain, \partial\Omega denotes the boundary of the domain such that these parameters verify the observation u_\text{obs} .

In other words, the goal is to minimize an objective function J(\kappa_0,\kappa_1,\kappa_2)

\min_{\kappa_0,\kappa_1,\kappa_2\in\mathbb{R}} J = \frac{1}{2}\int_{U} (u-u_\text{obs})^2

with U\in\Omega an open subset of all observable measures.

In practice, observations are often boundary measurement taken on a subset of the domain boundary U\subseteq\partial\Omega . Thus arise concerns about existence and uniqueness of solution for the inverse problem, especially if observations does not cover the whole boundary.

For the following example, we take U=\Omega and we define

\kappa=1+{\kappa_0}_{\vert_{I_0}} + {\kappa_1}_{\vert_{I_1}} + {\kappa_2}_{\vert_{I_2}}

where I_0, I_1, I_2 are different subset (inclusions) of \Omega .

55.1. Adjoint method

We add a small perturbation to the given parameter b, c, d denoted

\begin{aligned} \kappa_0^\delta = \kappa_0 + \alpha\delta \kappa_0 \\ \kappa_1^\delta = \kappa_1 + \alpha\delta \kappa_1 \\ \kappa_2^\delta = \kappa_2 + \alpha\delta \kappa_2 \end{aligned}

for \alpha\in\mathbb{R} . Then u^\delta=u+\alpha\delta u is the perturbed solution of the equation

\begin{cases} \begin{aligned} \displaystyle -\nabla\cdot(\kappa^\delta \nabla u^\delta) &= 0 \quad \text{on}\;\Omega \\ u^\delta &= 0 \quad \text{on}\;\partial\Omega \end{aligned} \end{cases}

where the diffusion coefficient for the perturbation \kappa^\delta \rightarrow \kappa( \kappa_0^\delta, \kappa_1^\delta, \kappa_2^\delta )

Lets note \delta J the shift for the cost function J such that

\begin{aligned} \delta J &=J(\kappa_0^\delta, \kappa_1^\delta, \kappa_2^\delta) - J(\kappa_0,\kappa_1,\kappa_2) \\ &=\frac{1}{2}\int_U (u^\delta - u_\text{obs})^2 - (u-u_\text{obs})^2 \\ &=\frac{1}{2}\int_U (u^\delta - u)(u^\delta+u-2 u_\text{obs}) \\ &=\frac{1}{2}\int_U (u^\delta - u)( (u^\delta-u_\text{objs}) + (u - u_\text{obs}) ) \\ \end{aligned} If we divide by \alpha and look when \alpha\rightarrow 0 , then we have

\delta J \approx \int_U \delta u(u - u_\text{obs})

Now we desire a method to compute the gradient. We consider the following linear tangent model deduced from two previous system of equations. We have

-\nabla\cdot(\kappa^\delta \nabla u^\delta) = -\nabla\cdot(\kappa \nabla u)

We introduce a state variable p (the adjoint) chosen appropriately as we will see further. We multiply the previous model by this variable and integrate on the domain.

-\int_\Omega \nabla\cdot(\kappa^\delta \nabla u^\delta) p = -\int_\Omega \nabla\cdot(\kappa \nabla u) p

Then writing the weak form we have

\int_\Omega \kappa^\delta \nabla u^\delta \nabla p -\int_\Omega\kappa^\delta \frac{\partial u^\delta}{\partial\mathbf n}p = \int_\Omega \kappa\nabla u\nabla p -\int_\Omega\kappa\frac{\partial u}{\partial\mathbf n}p

we expand the perturbed diffusion coefficient as \kappa^\delta=\kappa +\alpha\delta\kappa and rearrange terms by u^\delta -u ( =\alpha\delta u ). We divide by \alpha as previously and makes \alpha\rightarrow 0 , we have

\int_\Omega \kappa \nabla \delta u \nabla p -\int_{\partial\Omega} \kappa \frac{\partial \delta u }{\partial\mathbf n}p \approx -\int_\Omega \delta\kappa\nabla u \nabla p +\int_{\partial\Omega} \delta\kappa \frac{\partial u}{\partial\mathbf n}p

We integrate by part a second time

-\int_\Omega \delta u \nabla\cdot(\kappa \nabla p) +\int_{\partial\Omega} \left( \kappa \frac{\partial p }{\partial\mathbf n}\delta u - \kappa \frac{\partial \delta u }{\partial\mathbf n}p \right) \approx \int_\Omega u \nabla\cdot(\delta\kappa \nabla p) -\int_{\partial\Omega} \left( \delta\kappa \frac{\partial p}{\partial\mathbf n}u +\delta\kappa \frac{\partial u}{\partial\mathbf n}p \right)

we obtain finally for p=0 on the boundary

-\int_\Omega \delta u \nabla\cdot(\kappa \nabla p) \approx \int_\Omega u\nabla\cdot(\delta\kappa \nabla p)

Now to compute the gradient, we recall the previous equation for J and we deduce from build from the previous equation the following model (adjoint equation) for the state variable p

\begin{cases} \begin{aligned} \displaystyle -\nabla\cdot(\kappa\nabla p) &= (u - u_\text{objs}) \\ u &= 0 \quad \text{on}\;\partial\Omega \end{aligned} \end{cases}

If we multiply this equation by u^\delta-u and integrate over the domain, we obtain

\begin{aligned} -\int_\Omega \nabla\cdot(\kappa\nabla p)(u^\delta-u) &= \int_\Omega (u - u_\text{objs})(u^\delta - u) \\ &\approx \int_\Omega (u - u_\text{objs})\delta u \\ &\approx \delta J \end{aligned}

Now if we write the weak form of this equation, we have

\delta J \approx\int_\Omega \kappa\nabla p\nabla \delta u -\int_{\partial\Omega} \kappa\frac{\partial p}{\partial\mathbf n} \delta u

We apply the boundary condition

\begin{aligned} \delta J &=\int_\Omega \kappa\nabla p\nabla\delta u \\ &=-\int_\Omega \delta\kappa\nabla p\nabla u \end{aligned}

We recall \kappa chosen as

\kappa=1+{k_0}_{\vert_{I_0}} + {k_1}_{\vert_{I_1}} + {k_2}_{\vert_{I_2}}

Then the we deduce the gradient from the derivatives

\nabla J(\kappa_0,\kappa_1,\kappa_2) = \left( -\int_{I_0} \nabla p\nabla u, -\int_{I_1} \nabla p\nabla u, -\int_{I_2} \nabla p\nabla u \right)

So now we have a method to compute the gradient.

56. Applications

57. Mesh

As the title suggest, this section will be dedicate to mesh applications. It is composed by

57.1. feelpp_mesh_partitioner

feelpp_mesh_partitioner is a simple application which can generate a partitioned mesh and save it in a Feel++ specific json+hdf5 file format. The generated mesh can then be loaded very efficiently in parallel.

57.1.1. Options

feelpp_mesh_partitioner requires some options.

Table 84. Table of command-line feelpp_mesh_partitioner options

Name

Description

Default value

dim

dimension of the mesh

3

shape

shape of the mesh elements

Simplex

part

number of desired partitions

ifile

name or path to the mesh

ofile

output filename prefix (without extension)

odir

output directory

57.1.2. Examples

We are now using the feelpp/toolboxes:latest docker images as described in Getting Started to demonstrate feelpp_mesh_partitioner usage. We use the meshes in src/feelpp/data/gmsh/primitives in the docker image.

Starting the docker feelpp/toolboxes:latest
docker run -ti -v $HOME/feel:/feel feelpp/toolboxes:latest

57.1.3. Generating a single mesh partitioning

We generate a mesh partitioned on 4 cores with the following command

feelpp_mesh_partitioner --part 4  --ifile src/feelpp/data/gmsh/primitives/torus.geo --ofile torus

You should have in the current directory 4 files

ls torus*
torus.geo  torus.h5  torus.json  torus.msh

Now the file torus.json can be loaded in a Feel++ application distributed on 4 cores.

Generating multiple mesh partitioning

Often we are interested in a set of partitioned meshes in order to do a speed-up study. feelpp_mesh_partitioner is the right tool for that. Let’s consider a set of partitions \(\mathcal{P}=\{2, 4, 8, 16, 32 \}\).

feelpp_mesh_partitioner --part 2 4 8 16 32  --ifile src/feelpp/data/gmsh/primitives/torus.geo --odir torus

You should have in the directory torus (thanks to the odir option) 5 partitioned meshes

ls torus*
torus_p16.h5    torus_p2.h5    torus_p32.h5    torus_p4.h5    torus_p8.h5
torus_p16.json  torus_p2.json  torus_p32.json  torus_p4.json  torus_p8.json
The mesh filenames contain the partition information.

58. Reading

When we have partitoned a mesh with feelpp_mesh_partitioner, we would naturally like to read the output file in order to extract data and work on. That’s the purpose of mesh.cpp.

58.1. Code

From the file given by the option mesh.filename, the programm can extract the data and load the associated mesh into our execution. With this, we can give a json file, obtained from feelpp_mesh_partitioner and it will retrieve the parts of the original mesh.

Appendix A: Bibliography