Feel++

Modeling with Feel++

Advection Models

JSON file

Let’s describe how the json file is build.

Properties

Table 1. Table of Models for model option
Model Name

Advection

Advection

Stokes

Stokes

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 2. Boundary conditions
Name Options Type

Dirichlet

faces, edges and component-wise

"Dirichlet"

Neumann

scalar, vectorial

"Neumann_scalar" or "Neumann_vectorial"

Action

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

Feel++ code

Here is the code

First at all, we define our model type with

typedef FeelModels::Advection<
        Simplex<FEELPP_DIM,1>,
        Lagrange<OrderAdvection, Scalar,Continuous,PointSetFekete> > model_type;

This definition allows us to create our advection model object FM like this

auto Adv = model_type::New("advection");

The method New retrieve all data from the configuration and json files, as well build a mesh if need.

With this object, we can initialize our model parameters, such as boundaries conditions. Data on our model and on the numeric solver are then save and print on the terminal. This is made by

Adv->init();
Adv->printAndSaveInfo();

Now that our model is completed, we can solve the associated problem. To begin the resolution

Adv->isStationary()

determine if our model is stationary or not.

If it is, then we need to solve our system only one time and export the obtained results.

Adv->solve();
Adv->exportResults();

If it’s not, our model is time reliant, and a loop on time is necessary. Our model is then solved and the results are export at each time step.

 for ( ; !Adv->timeStepBase()->isFinished(); Adv->updateTimeStep() )
 {
   Adv->solve();
   Adv->exportResults();
 }

Code

{% include "../Examples/advection_model.cpp" %}

Configuration file

The config file is used to define options, linked to our case, we would have the possibility to change at will. It can be, for example, files paths as follows

#      advection
[advection]
geofile=$cfgdir/ring2d.geo
filename=$cfgdir/ring2d.json

[exporter]
directory=simul/advect2d/ring/data
format=ensightgold

It can also be resolution dependent parameters such as mesh elements size, methods used to define our problem and solvers.

#       time
[advection.bdf]
order=2
[ts]
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" %}

Json file

First at all, we define some general information like the name ( and short name ) and the model we would like to use

"Name": "Ring2d",
"ShortName": "Ring2d",
"Model": "Advection",

In this case, we have only boundary conditions to define. Here, we impose homogeneous Dirichlet conditions.

"BoundaryConditions":
    {
        "advection":
        {
            "Dirichlet":
            {
                "Bottom":
                {
                    "expr":"0"
                },
                "Left":
                {
                    "expr":"0"
                },
                "InnerCircle":
                {
                    "expr":"0"
                },
                "OuterCircle":
                {
                    "expr":"0"
                }
            }
        }
    }

Code

{% include "../Examples/ring2d.json" %}

Compilation/Execution

Once you’ve a build dir, you just have to realise the command make at

{buildir}/applications/models/advection

This will generate an executable named feelpp_application_advection_2d. To execute it, you need to give the path of the cfg file associated to your case, with --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

CFD Toolbox

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"

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

2.1. Generalised Newtonian fluid

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

The viscosity law is activated by:

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

3. Boundary Conditions

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

3.1. Dirichlet on velocity

A Dirichlet condition on velocity field reads:

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

3.2. Dirichlet on pressure

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

3.3. Neumann

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

3.4. Slip

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

3.5. Inlet

The boundary condition at inlets allow to define a velocity profile on a set of marked faces \(\Gamma_{\mathrm{inlet}}\) in fluid mesh:

\boldsymbol{u}_f = - g \ \boldsymbol{ n } \quad \text{ on } \Gamma_{\mathrm{inlet}}

The function \(g\) is computed from flow velocity profiles:

Constant profile

\text{Find } g \in C^0(\Gamma) \text{ such that } \\ \begin{eqnarray} g &=& \beta \quad &\text{ in } \Gamma \setminus \partial\Gamma \\ g&=&0 \quad &\text{ on } \partial\Gamma \end{eqnarray}

Parabolic profile

\text{Find } g \in H^2(\Gamma) \text{ such that : } \\ \begin{eqnarray} \Delta g &=& \beta \quad &\text{ in } \Gamma \\ g&=&0 \quad &\text{ on } \partial\Gamma \end{eqnarray}

where \(\beta\) is a constant determined by adding a constraint to the inflow:

velocity_max

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

flow_rate

\(\int_\Gamma ( g \ \boldsymbol{n} ) \cdot \boldsymbol{n} = \alpha\)

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

3.6. Outlet flow

Table 7. Outlet flow options
Option Value Default value Description

model

free,windkessel

free

select an outlet modeling

3.6.1. Free outflow

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

3.6.2. Windkessel model

We use a 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 8. 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

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

3.8. Options summary

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

4. Body forces

Body forces are also defined in BoundaryConditions category in json file.


"VolumicForces":
{
    "<marker>":
    {
        "expr":"{0,0,-gravityCst*7850}:gravityCst"
    }
}

The marker corresponds to mesh elements marked with this tag. If the marker is an empty string, it corresponds to all elements of the mesh.

5. Post Processing

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

5.1. Exports for vizualisation

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

  • velocity

  • pressure

  • displacement

  • vorticity

  • stress or normal-stress

  • wall-shear-stress

  • density

  • viscosity

  • pid

  • alemesh

5.2. Measures

  • Points

  • Force

  • FlowRate

  • Pressure

  • VelocityDivergence

Points

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

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

  • the coordinate of point

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

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

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

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

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

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

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

compute lift and drag

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

5.3. Export user functions

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

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

6. Stabilization methods

6.1. GLS family

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

6.2. CIP family

Documentation pending

7. Run simulation

programme avalaible :

  • feelpp_toolbox_fluid_2d

  • feelpp_toolbox_fluid_3d

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

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

8. 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}\]

9. Hyper-elasticity

9.1. Saint-Venant-Kirchhoff

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

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

9.2.1. Isochoric part : \(\boldsymbol{\Sigma}_s^\text{iso}\)

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

9.2.2. Volumetric part : \(\boldsymbol{\Sigma}_s^\text{vol}\)

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

10. Axisymmetric reduced model

We interest us here to a 1D reduced model, named generalized string.

The axisymmetric form, which will interest us here, is a tube of length \(L\) and radius \(R_0\). It is oriented following the axis \(z\) and \(r\) represent the radial axis. The reduced domain, named \(\Omega_s^*\) is represented by the dotted line. So the domain, where radial displacement \(\eta_s\) is calculated, is \(\Omega_s^*=\lbrack0,L\rbrack\).

We introduce then \(\Omega_s^{'*}\), where we also need to estimate a radial displacement as before. The unique variance is this displacement direction.

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.

CSM Toolbox

11. Models

The solid mechanics model can be selected in json file :

Listing : select solid model
"Model":"Hyper-Elasticity"
Table 13. 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 14. 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

12. Materials

The Lamé coefficients are deducing from the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material we work on and can be express

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

13. Boundary Conditions

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

14. Body forces

Table 16. Volumic forces
Name Options Type

Expression

Vectorial

"VolumicForces"

15. Post Process

15.1. Exports for visualisation

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

  • displacement

  • velocity

  • acceleration

  • stress or normal-stress

  • pressure

  • material-properties

  • pid

  • fsi

  • Von-Mises

  • Tresca

  • principal-stresses

  • all

15.2. Measures

  • Points

  • Maximum

  • Minimum

  • VolumeVariation

15.2.1. 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, …​

15.2.2. Maximum/Minimum

The Maximum and minimum can be evaluated and save on .csv file. User need to define (i) <Type> ("Maximum" or "Minimum"), (ii) "<tag>" representing this data in the .csv file, (iii) "<marker>" representing the name of marked entities and (iv) the field where extremum is computed.

"<Type>":
{
    "<tag>":
    {
        "markers":"marker>",
        "fields":["displacement","velocity"]
    }
}

15.2.3. VolumeVariation

"VolumeVariation":<marker>

16. Run simulations

programme avalaible :

  • feelpp_toolbox_solid_2d

  • feelpp_toolbox_solid_3d

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

Fluid Structure Interaction Models

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

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

18.1. Fluid structure coupling conditions with 1D reduced model

For the coupling conditions, between the 2D fluid and 1D structure, we need to modify the original ones \(\boldsymbol{(1)},\boldsymbol{(2)}, \boldsymbol{(3)}\) by

\[\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)}\]

18.2. Variables, symbols and units

Notation

Quantity

Unit

\(\boldsymbol{u}_f\)

fluid velocity

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

\(\boldsymbol{\sigma}_f\)

fluid stress tensor

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

\(\boldsymbol{\eta}_s\)

displacement

\(m\)

\(\boldsymbol{F}_s\)

deformation gradient

dimensionless

\(\boldsymbol{\Sigma}_s\)

second 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})\]

Thermo-Electric Toolbox

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

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

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

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

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

Unresolved directive in README.adoc - include::Electromagnetism/README.adoc[]