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:

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.

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.

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

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

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.

programme avalaible :

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

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

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

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

where \(\eta_s\), the radial displacement that satisfy this equation, \(k\) is the Timoshenko’s correction factor, and \(\gamma_v\) is a viscoelasticity parameter. The material is defined by its density \(\rho_s^*\), its Young’s modulus \(E_s\), its Poisson’s ratio \(\nu_s\) and its shear modulus \(G_s\)

At the end, we take \( \eta_s=0\text{ on }\partial\Omega_s^*\) as a boundary condition, which will fixed the wall to its extremities.

The solid mechanics model can be selected in json file :

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

option: mechanicalproperties.compressible.volumic_law

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

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

programme avalaible :

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 :

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

Let’s describe how the json file is build.

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

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

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

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

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

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.

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

where \(k\) is the thermal conductivity.

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

We consider that the thermal exchanges are limited to the cooled region, which amounts to apply an homogeneous Neumann condition on all other boundaries.

In the thermo-electric model, we have 4 parameters

We gather in the following table parameter ranges, nominal values as well as units for \(\sigma, k, T_c, h\).

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

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

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

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

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

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:

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

We start with the mesh

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.

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.

next we solve the algebraic problem

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.

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

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

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

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

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

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

:leveloffset:-1

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.

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.

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.

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

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

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

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

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

The result files are then stored by default in

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

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

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.

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

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 :

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.

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

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

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

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.

This is the parameters associate to the two cases, which interest us here.

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:

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

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

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.

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.

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.

.

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

It’s define here into a Lagrangian frame.

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

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

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

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

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

where \(\eta_s\), the radial displacement that satisfy this equation, \(k\) is the Timoshenko’s correction factor, and \(\gamma_v\) is a viscoelasticity parameter. The material is defined by its density \(\rho_s^*\), its Young’s modulus \(E_s\), its Poisson’s ratio \(\nu_s\) and its shear modulus \(G_s\)

At the end, we take \( \eta_s=0\text{ on }\partial\Omega_s^*\) as a boundary condition, which will fixed the wall to its extremities.

inlude::FSI/README.adoc[]

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:

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:

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

Only the 2D specifications have been implemented.

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

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

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:

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

We start with the mesh

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.

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.

next we solve the algebraic problem

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.

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

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

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

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

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

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

:leveloffset:-1

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

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

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:

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

We start with the mesh

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.