Computational Fluid Dynamics

Welcome to the Computational Fluid Dynamics(CFD) toolbox. This module explains how to use the Laminar Flow interface to model and simulate fluid mechanics for laminar, incompressible fluids.

The engineering community often uses the term CFD, computational fluid dynamics, to refer to the numerical simulation of fluids.

1. Features

2D and 3D
Stokes flows
Laminar Incompressible flows
Steady and unsteady simulations
BDF time schemes up to order 4
Moving domain support using Arbitrary Lagrangian Eulerian (ALE) formulation
Support high order geometry including in context of moving domain
Stabilization pressure and advection dominated using Galerking Least Square
Boundary conbditions : no-slip, slip(symmetry), inflow(not necessarily aligned with an axis), pressure, outflow

2. CFD Toolbox

2.1. Models

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

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

Listing : select fluid model

"Models": { "equations":"Navier-Stokes" }

2.2. Materials

Table 1. Materials properties defined in the heat toolbox (symbols are given without components suffix)
Name	Symbol	Shape	Description
density	rho	scalar	density
dynamic-viscosity	mu	scalar	dynamic viscosity
turbulent-dynamic-viscosity	mu_t	scalar	turbulent dynamic viscosity
consistency-index	mu_k	scalar	consistency index
power-law-index	mu_power_law_n	scalar	power law index
viscosity-min	mu_min	scalar
viscosity-max	mu_max	scalar
viscosity-zero-shear	mu_0	scalar
viscosity-infinite-shear	mu_inf	scalar
carreau-law-lambda	mu_carreau_law_lambda	scalar
carreau-law-n	mu_carreau_law_n	scalar
carreau-yasuda-law-lambda	mu_carreau_yasuda_law_lambda	scalar
carreau-yasuda-law-n	mu_carreau_yasuda_law_n	scalar
carreau-yasuda-law-a	mu_carreau_yasuda_law_a	scalar

2.3. Fields and symbols expressions

Table 2. Fields available
Name	Description	Shape
velocity	the velocity of fluid	vectorial
pressure	the pressure	scalar

Table 3. Symbols expressions by considering `fluid` as the keyword used with the toolbox
Symbol	Expression	Description
fluid_U_0	\(u_0\)	evaluate the component X of fluid velocity
fluid_U_1	\(u_1\)	evaluate the component Y of fluid velocity
fluid_U_2	\(u_2\)	evaluate the component Z of fluid velocity
fluid_U_magnitude	\(\lVert \boldsymbol{u} \rVert = \sqrt{\boldsymbol{u} \cdot \boldsymbol{u}}\)	evaluate the magnitude of fluid velocity
fluid_curl_U	\(\nabla \wedge \boldsymbol{u}\)	evaluate the curl of fluid velocity
fluid_curl_U_magnitude	\(\lVert \nabla \wedge \boldsymbol{u} \rVert\)	evaluate the magnitude of curl of fluid velocity
fluid_div_U	\(\nabla \cdot \boldsymbol{u}\)	evaluate the div of fluid velocity

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 4. 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":
{
    "<name>"
    {
        "markers":"[marker1,marker2]",
        "rho":"1.0e3",
        "mu":"1.0"
    }
}

2.3.1. Generalised Newtonian fluid

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

The viscosity law is activated by:

Table 5. 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

Viscosity law

options

unit

power_law

power_law.k

power_law.n

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

l/g

carreau_law

viscosity.zero_shear

viscosity.infinite_shear

carreau_law.lambda

carreau_law.n

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

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

dimensionless

2.4. Boundary Conditions

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

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

2.4.2. Dirichlet on pressure

\[\begin{align} p & = g \\ \boldsymbol{u}_f \times {\boldsymbol{ n }} & = \boldsymbol{0} \end{align}\]

2.4.3. Neumann

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

2.4.4. Slip

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

2.4.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 7. Inlet flow options
Option	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

2.4.6. Outlet flow

Table 8. 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 9. 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

2.4.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>",
                // ...
            }
        }
    }
}

2.4.8. Options summary

Table 10. 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

2.5. Body forces

Body forces are also defined in the fluid section of BoundaryConditions in the json file.

"BoundaryConditions":
{
    "fluid":{
        "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.

2.6. Post Processing

"PostProcess":
{
    "Exports":
    {
        "fields":["field1","field2",...]
    },
    "Measures":
    {
        "<measure type>":
        {
            "label":
            {
                "<range type>":"value",
                "fields":["field1","field3"]
            }
        }
    }
}

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

2.6.2. Measures

Points
Force
FlowRate
Pressure
VelocityDivergence

The following fluid variables are available

Name

Description

fluid_U_magnitude

magnitude of the velocity

fluid_Ux

x component of the velocity

fluid_Uy

y component of the velocity

fluid_Uz

z component of the velocity

fluid_P

pressure

You can use these expressions in Statistics

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"]

2.6.3. Export user functions

A function defined by a symbolic expression can be represented as a finite element field thanks to nodal projection. This function can be exported.

"Functions":
{
    "toto":{ "expr":"x*y:x:y"},
    "toto2":{ "expr":"0.5*ubar*x*y:x:y:ubar"},
    "totoV":{ "expr":"{2*x,y}:x:y"}
},
"PostProcess":
{
   "Exports":
   {
       "fields":["velocity","pressure","pid","totoV","toto","toto2"]
   }
}

2.7. Stabilization methods

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

2.7.2. CIP family

Documentation pending

2.8. Run simulation

The computational fluid dynamics applications available are

2D CFD toolbox : feelpp_toolbox_fluid_2d
3D CFD toolbox : feelpp_toolbox_fluid_3d

Here is an example of execution of the 2D CFD toolbox on 4 processors using the configuration script <myfile.cfg>

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