Modeling
1. Notations and units
Notation 
Quantity 
Unit 
\(\boldsymbol{\eta}_s\) 
displacement 
\(m\) 
\(\rho_s\) 
density 
\(kg.m^{3}\) 
\(\lambda_s\) 
first Lamé coefficients 
\(N.m^{2}\) 
\(\mu_s\) 
second Lamé coefficients 
\(N.m^{2}\) 
\(E_s\) 
Young modulus 
\(kg.m^{1}.s^{2}\) 
\(\nu_s\) 
Poisson’s ratio 
dimensionless 
\(\boldsymbol{F}_s\) 
deformation gradient 

\(\boldsymbol{\Sigma}_s\) 
second PiolaKirchhoff tensor 

\(f_s^t\) 
body force 

strain tensor \(\boldsymbol{F}_s = \boldsymbol{I} + \nabla \boldsymbol{\eta}_s\)

CauchyGreen tensor \(\boldsymbol{C}_s = \boldsymbol{F}_s^{T} \boldsymbol{F}_s\)

GreenLagrange tensor
2. Equations
The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows
2.1. Linear elasticity
2.2. Hyperelasticity
2.2.1. SaintVenantKirchhoff
2.2.2. NeoHookean
Isochoric part : \(\boldsymbol{\Sigma}_s^\text{iso}\)
Name  \(\mathcal{W}_S(J_s)\)  \(\boldsymbol{\Sigma}_s^{\text{iso}}\) 

NeoHookean 
\(\mu_s J^{2/3}(\boldsymbol{I}  \frac{1}{3} \text{tr}(\boldsymbol{C}) \ \boldsymbol{C}^{1}) \) 
Volumetric part : \(\boldsymbol{\Sigma}_s^\text{vol}\)
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)\) 
2.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.
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.
3. CSM Toolbox
3.1. Models
The solid mechanics model can be selected in json file :
"Model":"HyperElasticity"
Model  Name in json 

Linear Elasticity 

Hyper Elasticity 

When materials are closed to incompressibility formulation in displacement/pressure are available.
Model  Name  Volumic law 

SaintVenantKirchhoff 

classic, simo1985 
NeoHookean 

classic, simo1985 
option: mechanicalproperties.compressible.volumic_law
3.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
"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}\).
3.3. Boundary Conditions
Name  Options  Type 

Dirichlet 
faces, edges and componentwise 
"Dirichlet" 
Neumann 
scalar, vectorial 
"Neumann_scalar" or "Neumann_vectorial" 
Pressure follower , 
Nonlinear boundary condition set in deformed domain 
TODO 
Robin 
TODO 
TODO 
3.4. Body forces
Name  Options  Type 

Expression 
Vectorial 
"VolumicForces" 
3.5. Post Process
3.5.1. Exports for visualisation
The fields allowed to be exported in the Fields
section are:

displacement

velocity

acceleration

stress or normalstress

pressure

materialproperties

pid

fsi

VonMises

Tresca

principalstresses

all
3.5.2. Measures

Points

Maximum

Minimum

VolumeVariation
Points
Same syntax as FluidMechanics with available Fields :

displacement

velocity

acceleration

pressure

principalstress0

principalstress1

principalstress2

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>
3.6. Run simulations
programme avalaible :

feelpp_toolbox_solid_2d

feelpp_toolbox_solid_3d
mpirun np 4 feelpp_toolbox_solid_2d configfile=<myfile.cfg>
Benchmarks
4. Turek Hron CSM Benchmark
In order to validate our fluidstructure 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. 
4.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
\(\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 SaintVenantKirchhoff 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 quasistatic 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.
4.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^*\).
4.1.2. Initial conditions
The source term \(\boldsymbol{f}_s\), chosen in order to set in motion the elastic structure, is define by
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\).
4.2. Inputs
The following table displays the various fixed and variables parameters of this testcase.
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 nonlinear part and a direct method based on LU decomposition is selected for the linear part.
4.3. Outputs
As described before, we have
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.
4.4. Discretization
To realize these tests, we made the choice to used Finite Elements Method, with Lagrangian elements of order \(N\) to discretize space.
Newmarkbeta 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\).
4.4.1. Solvers
Here are the different solvers ( linear and nonlinear ) used during results acquisition.
type 
gmres 
relative tolerance 
1e13 
max iteration 
1000 
reuse preconditioner 
true 
relative tolerance 
1e8 
steps tolerance 
1e8 
max iteration 
500 
max iteration with reuse 
10 
reuse jacobian 
false 
reuse jacobian rebuild at first Newton step 
true 
relative tolerance 
1e5 
max iteration 
500 
reuse preconditioner 
CSM1/CSM2 : false  CSM3 : true 
reuse preconditioner rebuild at first Newton step 
true 
type 
lu 
package 
mumps 
4.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
feelppbenchmarkingbook/CFD/TurekHron/postprocess_cfd.py
4.6. Results
4.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 ]
4.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].
4.6.3. CSM3
The results of the CSM3 benchmark are detailed below.
\(\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] 
\(\) \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 ]
4.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.
4.7. Bibliography

[TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluidstructure 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.
5. 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.
5.1. Definition
5.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 :
We will use the following notations:
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 :
The revwriting of the equilibrium equation in cylindrical coordinates gives :
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:
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:
where \(\mu\) and \(\lambda\) are the Lamé coefficients, and the tensor of small deformation is given in cylindrical coordinates by:
Then, using the last two definitions and the properties of the solenoid (axisymmetric and infinitely long), we can rewrite the equilibrium equation as:
5.1.2. Analytical solution
We want to find an analytical solution of the form :
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 :
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 :
Replacing \(b_z\) with his expression in the equilibrium equation, this gives :
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:
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
we can solve the system to find the constants:
The final step is to translate this analytical solution \(u_{cyl}(r)\) into cartesian coordinates to obtain the analytical cartesian displacement \(\mathbf{u}_{cart}\):
5.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.
5.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\)
5.2. Inputs
We use the following parameters:
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}\) 
5.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 \).
5.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.
6. NAFEMS LE1 Benchmarck
This benchmark is extract from the Abaqus Benchmarks Manual.
6.1. Definition
We focus on the LE1 benchmarks in particular.
6.1.1. Geometry
The geometry is given here by :
6.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.
6.2. Inputs
We have the following parameters:
Name  Value 

\(E\) 
\(210\, GPa\) 
\(\nu\) 
\(0.3\) 
\(\rho\) 
\(7800\, kg/m^2\) 
6.3. Outputs
We want to compare the value of \(\sigma_{yy}\) at the point D. The reference value is \(92.7\, MPa\).
6.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.
7. NAFEMS LE10 Benchmarck
This benchmark is extract from the Abaqus Benchmarks Manual.
7.1. Definition
We focus on the LE10 benchmarks in particular.
7.1.1. Geometry
The geometry is given here by :
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'.
7.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.
7.2. Inputs
We have the following parameters:
Name  Value 

\(E\) 
\(210\, GPa\) 
\(\nu\) 
\(0.3\) 
\(\rho\) 
\(7800\, kg/m^2\) 
7.3. Outputs
We want to compare the value of \(\sigma_{yy}\) at the point D. The reference value is \(5.38\, MPa\).
7.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.