CFD Toolbox

Computational Fluid Dynamics

Benchmarks are available at CFD Benchmarks.

1. Notations and units

Notation Quantity Unit

\(\rho_f\)

fluid density

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

\(\boldsymbol{u}_f\)

fluid velocity

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

\(\boldsymbol{\sigma}_f\)

fluid stress tensor

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

\(\boldsymbol{f}^t_f\)

source term

\(kg \cdot m^{-3} \cdot s^{-1}\)

\(p_f\)

pressure fields

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

\(\mu_f\)

dynamic viscosity

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

\(\bar{U}\)

characteristic inflow velocity

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

\(\nu\)

kinematic viscosity

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

\(L\)

characteristic length

\(m\)

2. Equations

Navier-Stokes model is used to model incompressible Newtonian fluid. It can be described by these conservative laws :

Momentum conservation equation
\[\rho_{f} \left. \frac{\partial\mathbf{u}_f}{\partial t} \right|_\mathrm{x} + \rho_{f} \left( \boldsymbol{u}_{f} \cdot \nabla_{\mathrm{x}} \right) \boldsymbol{u}_{f} - \nabla_{\mathrm{x}} \cdot \boldsymbol{\sigma}_{f} = \boldsymbol{f}^t_f , \quad \text{ in } \Omega^t_f \times \left[t_i,t_f \right]\]
Mass conservation equation
\[\nabla_{\mathrm{x}} \cdot \boldsymbol{u}_{f} = 0, \quad \text{ in } \Omega^t_f \times \left[t_i,t_f \right]\]

we complete this set of equations with the fluid constitutive law

\[ \boldsymbol{\sigma}_{f} = -p_f \boldsymbol{I} + 2\mu_f D(\boldsymbol{u}_{f}) \]

with strain tensor \(D(\boldsymbol{u}_{f})\) defined by :

\[ D(\boldsymbol{u}_{f}) = \frac{1}{2} (\nabla_\mathrm{x} \boldsymbol{u}_f + (\nabla_\mathrm{x} \boldsymbol{u}_f)^T) \]

An alternative model is the Stokes model. It is valid in the case of small Reynolds number. It corresponds to the same formulation than Navier-Stokes equations but without the convective term \(\left( \boldsymbol{u}_{f} \cdot \nabla_{\mathrm{x}} \right) \boldsymbol{u}_{f}\) .

Multifluid flows

We introduce here the methodology associated to multifluid flows.

3. Introduction

Let’s define a bounded domain \(\Omega \subset \mathbb{R}^p\) (\(p=2,3\)) decomposed into two subdomains \(\Omega_1\) and \(\Omega_2\). We denote \(\Gamma\) the interface between the two partitions. The goal of the level set method is to track implicitly the interface \(\Gamma(t)\) moving at a velocity \(\mathbf{u}\). The level set method has been described in [osher] and its main ingredient is a continuous scalar function \(\phi\) (the /level set/ function) defined on the whole domain. This function is chosen to be positive in \(\Omega_1\), negative in \(\Omega_2\) and zero on \(\Gamma\). The motion of the interface is then described by the advection of the level set function with a divergence free velocity field \(\mathbf{u}\):

\[ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = 0,\quad \nabla \cdot \mathbf{u} = 0.\]

A convenient choice for \(\phi\) is a signed distance function to the interface. Indeed, the property \(|\nabla \phi| = 1\) of distance functions eases the numerical solution and gives a convenient support for delta and Heaviside functions).

Nevertheless, it is known that the advection equation \eqref{eq:advection} does not conserve the property \(|\nabla \phi|=1\). Thus, when \(|\nabla \phi|\) is far from \(1\) we have to reset \(\phi\) as a distance function without moving the interface. To do so we can either use an Hamilton-Jacobi method or the fast marching method (see \cite{Winkelmann2007} for details about the fast marching method).

In two-fluid flow simulations, we need to define some quantities related to the interface such as the density, the viscosity, or some interface forces. To this end, we introduce the smoothed Heaviside and delta functions:

\[ H_{\varepsilon}(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \displaystyle\frac{1}{2} \left(1+\frac{\phi}{\varepsilon}+\frac{\sin(\frac{\pi \phi}{\varepsilon})}{\pi}\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 1, & \phi \geq \varepsilon. \end{array} \right.\]
\[ \delta_{\varepsilon}(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq - \varepsilon,\\ \displaystyle\frac{1}{2 \varepsilon} \left(1+\cos(\frac{\pi \phi}{\varepsilon})\right), & -\varepsilon \leq \phi \leq \varepsilon, \\ 0, & \phi \geq \varepsilon. \end{array} \right.\]

where \(\varepsilon\) is a parameter defining a ``numerical thickness'' of the interface. A typical value of \(\varepsilon\) is \(1.5 h\) where \(h\) is the mesh size of elements crossed by the iso-value \(0\) of the level set function.

The Heaviside function is used to define parameters having different values on each subdomains. For example, we define the density of two-fluid flow as \(\rho = \rho_2 + (\rho_1-\rho_2) H_{\varepsilon}(\phi)\) (we use a similar expression for the viscosity \(\nu\)). Regarding the delta function, it is used to define quantities on the interface. In particular, in the variational formulations, we replace integrals over the interface \(\Gamma\) by integrals over the entire domain \(\Omega\) using the smoothed delta function. If \(\phi\) is a signed distance function, we have : \(\int_{\Gamma} 1 \simeq \int_{\Omega} \delta_{\varepsilon}(\phi)\). If \(\phi\) is not close enough to a distance function, then \(\int_{\Gamma} 1 \simeq \int_{\Omega} |\nabla \phi| \delta_{\varepsilon}(\phi)\) which still tends to the measure of \(\Gamma\) as \(\varepsilon\) vanishes. However, if \(\phi\) is not a distance function, the support of \(\delta_{\varepsilon}\) can have a different size on each side of the interface. More precisely, the support of \(\delta_{\varepsilon}\) is narrowed on the side where \(|\nabla \phi|>1\) and enlarged on regions where \(|\nabla \phi|<1\). It has been shown in [cottet] that replacing \(\phi\) by \(\frac{\phi}{|\nabla \phi|}\) has the property that \(|\nabla \frac{\phi}{|\nabla \phi|}| \simeq 1\) near the interface and has the same iso-value \(0\) as \(\phi\). Thus, replacing \(\phi\) by \(\frac{\phi}{|\nabla \phi|}\) as support of the delta function does not move the interface. Moreover, the spread interface has the same size on each part of the level-set \(\phi=0\). It reads \(\int_{\Gamma} 1 \simeq \int_{\Omega} \delta_{\varepsilon}(\frac{\phi}{| \nabla \phi|})\). The same technique is used for the Heaviside function.

5. Numerical implementation and coupling with the fluid solver

We use Feel++ to discretize and solve the problem. Equation \eqref{eq:advection} is solved using a stabilized finite element method. We have implemented several stabilization methods such as Streamline Upwind Diffusion (SUPG), Galerkin Least Square (GLS) and Subgrid Scale (SGS). A general review of these methods is available in [Franca92]. Other available methods include the Continuous Interior Penalty method (CIP) are described in \cite{Burman2006}. The variational formulation at the semi-discrete level for the stabilized equation \eqref{eq:advection} reads, find \(\phi_h \in {\mathbb R}_h^k\) such that \(\forall \psi_h \in {\mathbb R}_h^k\) :

\[ \left(\int_{\Omega} \frac{\partial \phi_h}{\partial t} \psi_h + \int_{\Omega} (\mathbf{u}_h \cdot \nabla \phi_h) \psi_h\right) + S(\phi_h, \psi_h) = 0,\]

where \(S(\cdot, \cdot)\) stands for the stabilization bilinear form (see section \ref{sec:membr-inext} for description of \({\mathbb R}_h^k\) and \(\mathbf{u}_h\)). In our case, we use a BDF2 scheme which needs the solution at the two previous time step to compute the one at present time. For the first time step computation or after a reinitialization we use an Euler scheme.

6. Two-fluid flow

The two-fluid flow problem can be expressed by

\[\begin{align} \frac{D( \rho_\phi \boldsymbol{u} )}{Dt} - \boldsymbol{\nabla} \cdot ( 2 \mu_\phi \boldsymbol{D}({\boldsymbol{u}}) ) + \boldsymbol{\nabla} p = \boldsymbol{f}_\phi \\ \nabla \cdot \boldsymbol{u} = 0 \\ \partial_t \phi + \boldsymbol{u} \cdot \boldsymbol{\nabla} \phi =0 \end{align}\]

Where \(\boldsymbol{f}_\phi\) is the force obtain by projection of the density of interfacial forces on the domain \(\Omega\)

\[\boldsymbol{f}_{\phi} = \boldsymbol{g}(\phi, \boldsymbol{ n }, \kappa) \delta_\varepsilon(\phi)\]

CFD Toolbox

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

8. 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"
    }
}

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

9. Boundary Conditions

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

9.1. Dirichlet on velocity

A Dirichlet condition on velocity field reads:

\[ \boldsymbol{u}_f = \boldsymbol{g} \quad \text{ on } \Gamma \]

or only a component of vector \(\boldsymbol{u}_f =(u_f^1,u_f^2 ,u_f^3 )\)

\[ u_f^i = g \quad \text{ on } \Gamma \]

Several methods are available to enforce the boundary condition:

  • elimination

  • Nitsche

  • Lagrange multiplier

9.2. Dirichlet on pressure

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

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

9.4. Slip

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

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

9.6. Outlet flow

Table 5. Outlet flow options
Option Value Default value Description

model

free,windkessel

free

select an outlet modeling

9.6.1. Free outflow

\[ \boldsymbol{\sigma}_{f} \boldsymbol{ n } = \boldsymbol{0} \]

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

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

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

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

11. Post Processing

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

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

11.2. Measures

  • Points

  • Force

  • FlowRate

  • Pressure

  • VelocityDivergence

Points

In order to evaluate velocity or pressure at specific points and save the results in .csv file, the user must define:

  • "<tag>" representing this data in the .csv file

  • the coordinate of point

  • the fields evaluated ("pressure" or "velocity")

"Points":
{
  "<tag>":
  {
    "coord":"{0.6,0.2,0}",
    "fields":"pressure"
  },
 "<tag>":
  {
    "coord":"{0.15,0.2,0}",
    "fields":"velocity"
  }
}
Flow rate

The flow rate can be evaluated and save on .csv file. The user must define:

  • "<tag>" representing this data in the .csv file

  • "<face_marker>" representing the name of marked face

  • the fluid direction ("interior_normal" or "exterior_normal") of the flow rate.

"FlowRate":
{
    "<tag>":
    {
        "markers":"<face_marker>",
        "direction":"interior_normal"
    },
    "<tag>":
    {
        "markers":"<face_marker>",
        "direction":"exterior_normal"
    }
}
Forces

compute lift and drag

"Forces":["fsi-wall","fluid-cylinder"]

11.3. Export user functions

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

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

12. Stabilization methods

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

12.2. CIP family

Documentation pending

13. Run simulation

programme avalaible :

  • feelpp_toolbox_fluid_2d

  • feelpp_toolbox_fluid_3d

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