Toolbox is available at FSI Toolbox. 
We will interest now to the different interactions a fluid and a structure can have together with specific conditions.
To describe and solve our fluidstructure interaction problem, we need to define a model, which regroup structure model and fluid model parts.
We have then in one hand the fluid equations, and in the other hand the structure equations.
The solution of this model are \((\mathcal{A}^t, \boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s)\).
Generally, the solid mechanic equations are expressed in a Lagrangian frame, and the fluid part in Eulerian frame. To define and take in account the fluid domain displacement, we use a technique name ALE ( Arbitrary Lagrangian Eulerian ). This allow the flow to follow the fluidstructure interface movements and also permit us to have a different deformation velocity than the fluid one.
Let denote \(\Omega^{t_0}\) the calculation domain, and \(\Omega^t\) the deformed domain at time \(t\). As explain before, we want to conserve the Lagrangian and Eulerian characteristics of each part, and to do this, we introduce \(\mathcal{A}^t\) the ALE map.
This map give us the position of \(x\), a point in the deformed domain at time \(t\) from the position of \(x^*\) in the initial configuration \(\Omega^*\).
\(\mathcal{A}^t\) is a homeomorphism, i.e. a continuous and bijective application we can define as
We denote also \(\forall \mathbf{x}^* \in \Omega^*\), the application :
This ALE map can then be retrieve into the fluidstructure model.
In order to validate our fluidstructure interaction solver, a benchmark, initially proposed by [TurekHron], is realized.
Computer codes, used for the acquisition of results, are from Vincent [Chabannes].
Note This benchmark is linked to the Turek Hron CFD and Turek Hron CSM benchmarks.
In this case, we want to study the flow of a laminar incrompressible flow around an elastic obstacle, fixed to a rigid cylinder, we combine here the study realized on Turek Hron CFD benchmark and Turek Hron CSM benchmark. So a 2D model, representative of this state, is considered.
We denote \(\Omega_f^*\) the fluid domain, represented by a rectangle of dimension \( [0,2.5] \times [0,0.41] \). This domain is characterized by its density \(\rho_f\) and its dynamic viscosity \(\mu_f\). In this case, the fluid material we used is glycerin.
Furthermore, the model chosen to describe the flow is the incompressible NavierStockes model. It is defined by the conservation of momentum equation and the conservation of mass equation. We also have the material constitutive equation to take in account in this domain, as well as the harmonic extension operator for the fluid movement.
On the other side, the structure domain, named \(\Omega_s^*\), represent the hyperelastic bar, bound to the stationary cylinder. As we want to work in a Lagrangian frame, and by Newton’s second law, the fundamental equation of the solid mechanic will be used. For the model, that our hyperelastic material follows, we use the SaintVenantKirchhoff one, define with Lamé coefficients. These coefficients are obtained from the the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material.
All of these are then gather into the fluidstructure coupling system.
We set
on \(\Gamma_{in}^*\), an inflow Dirichlet condition : \( \boldsymbol{u}_f=(v_{in},0);\)
on \(\Gamma_{wall}^* \cup \Gamma_{circ}^*\), a homogeneous Dirichlet condition : \(\boldsymbol{u}_f=\boldsymbol{0};\)
on \(\Gamma_{out}^*\), a Neumann condition : \(\boldsymbol{\sigma}_f\boldsymbol{n}_f=\boldsymbol{0};\)
on \(\Gamma_{fixe}^*\), a condition that imposes this boundary to be fixed : \(\boldsymbol{\eta}_s=0;\)
on \(\Gamma_{fsi}^*\) :
where \(\boldsymbol{n}_f\) is the outer unit normal vector from \(\partial \Omega_f^*\).
In order to describe the flow inlet by \(\Gamma_{in}\), a parabolic velocity profile is used. It can be expressed by
where \(\bar{U}\) is the mean inflow velocity.
However, to impose a progressive increase of this velocity profile, we define
\[ v_{in} = \left\{ \begin{aligned} & v_{cst} \frac{1\cos\left( \frac{\pi}{2} t \right) }{2} \quad & \text{ if } t < 2 \\ & v_{cst} \quad & \text{ otherwise } \end{aligned} \right. \]
With t the time.
Finally, we don’t want a source term, so \(f_f\equiv 0\).
The following table displays the various fixed and variables parameters of this testcase.
Name  Description  Nominal Value  Units 

\(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 
\(5.6 \times 10^6\) 
\(kg.m^{1}s^{2}\) 
\(\nu_s\) 
Poisson’s ratio 
\(0.4\) 
dimensionless 
\(\rho_s\) 
structure density 
\(1000\) 
\(kg.m^{3}\) 
\(\nu_f\) 
kinematic viscosity 
\(1\times 10^{3}\) 
\(m^2.s^{1}\) 
\(\mu_f\) 
dynamic viscosity 
\(1\) 
\(kg.m^{1}.s^{1}\) 
\(\rho_f\) 
density 
\(1000\) 
\(kg.m^{3}\) 
\(f_f\) 
source term 
0 

\(\bar{U}\) 
mean inflow velocity 
2 
\(m.s^{1}\) 
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.
The quantities we observe during this benchmark are on one hand the lift and drag forces ( respectively \(F_L\) and \(F_D\) ), as well as the displacement, on \(x\) and \(y\) axis, of the point A is the second value that interest us here.
They are the solution of the fluidstructure system.
This system also give us the ALE map \(\mathcal{A}_f^t\).
To realize these tests, we made the choice to used \(P_N~~P_{N1}\) TaylorHood finite elements to discretize the space.
For the fluid time discretization, BDF, at order \(q\), is the method we have chosen.
And finally Newmarkbeta method is the one we used for the structure time discretization, with parameters \(\gamma=0.5\) and \(\beta=0.25\).
These methods can be retrieved in [Chabannes] papers.
Here are the different solvers ( linear and nonlinear ) used during results acquisition.
KSP 

case 
fluid 
solid 
type 
gmres 

relative tolerance 
\(1e13\) 

max iteration 
\(1000\) 

reuse preconditioner 
true 
SNES 

case 
fluid 
solid 
relative tolerance 
\(1e8\) 

steps tolerance 
\(1e8\) 

max iteration 
\(50\) 

max iteration with reuse 
\(50\) 
\(10\) 
reuse jacobian 
true 
false 
reuse jacobian rebuild at first Newton step 
false 
true 
SNES 

case 
fluid 
solid 
relative tolerance 
\(1e5\) 

max iteration 
\(1000\) 

max iteration with reuse 
\(1000\) 
\(10\) 
reuse preconditioner 
true 

reuse preconditioner rebuild at first Newton step 
true 
PC 

case 
fluid 
solid 
type 
LU 

package 
mumps 
FSI 

solver method 
fix point with Aitken relaxation 
tolerance 
\(1e6\) 
max iterations 
\(1000\) 
initial \(\theta\) 
\(0.98\) 
minimum \(\theta\) 
\(1e12\) 
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.
Let’s start with the main code, that can be retrieve in
feelpp/applications/models/fsi
The FSI3 configuration file is located at
feelpp/applications/models/fsi/TurekHron
The result files are then stored by default in
feel/applications/models/fsi/TurekHron/fsi3/ <velocity_space><pression_space><Geometric_order><OrderDisp><Geometric_order>/np_<processor_used>
For example, for the FSI3 case executed on \(4\) processors, with a \(P_2\) velocity approximation space, a \(P_1\) pressure approximation space, a geometric order of \(1\) for fluid part and a \(P_1\) displacement approximation space and geometric order equals to \(1\) for solid part, the path is
feel/applications/models/fsi/TurekHron/fsi3/P2P1G1P1G1/np_4
First at all, we will discretize the simulation parameters for the different cases studied.
\(N_{elt}\) 
\(N_{dof}\) 
\( [P^N_c(\Omega_{f,\delta}]^2 \times P^{N1}_c(\Omega_{f,\delta}) \times V^{N1}_{s,\delta}\) 
\(\Delta t\) 

15872 
304128 
0.00025 

(1) 
1284 
27400 
\( [P^4_c(\Omega_{f,(h,3)}]^2 \times P^3_c(\Omega_{f,(h,3)}) \times V^3_{s,(h,3)}\) 
0.005 
(2) 
2117 
44834 
\( [P^4_c(\Omega_{f,(h,3)}]^2 \times P^3_c(\Omega_{f,(h,3)}) \times V^3_{s,(h,3)}\) 
0.005 
(3) 
4549 
95427 
\( [P^4_c(\Omega_{f,(h,3)}]^2 \times P^3_c(\Omega_{f,(h,3)}) \times V^3_{s,(h,3)}\) 
0.005 
(4) 
17702 
81654 
\( [P^2_c(\Omega_{f,(h,1)}]^2 \times P^1_c(\Omega_{f,(h,1)}) \times V^1_{s,(h,1)}\) 
0.0005 
Then the FSI3 benchmark results are detailed below.
\(x\) displacement \( [\times 10^{3}] \) 
\(y\) displacement \( [\times 10^{3}] \) 
Drag 
Lift 

2.69 ± 2.53 [10.9] 
1.48 ± 34.38 [5.3] 
457.3 ± 22.66 [10.9] 
2.22 ± 149.78 [5.3] 

464.5 ± 40.50 
6.00 ± 166.00 [5.5] 

2.88 ± 2.72 [10.9] 
1.47 ± 34.99 [5.5] 
460.5 ± 27.74 [10.9] 
2.50 ± 153.91 [5.5] 

4.54 ± 4.34 [10.1] 
1.50 ± 42.50 [5.1] 
467.5 ± 39.50 [10.1] 
16.2 ± 188.70 [5.1] 

474.9 ± 28.10 
3.90 ± 165.90 [5.5] 

2.83 ± 2.78 [10.8] 
1.35 ± 34.75 [5.4] 
458.5 ± 24.00 [10.8] 
2.50 ± 147.50 [5.4] 

(1) 
2.86 ± 2.74 [10.9] 
1.31 ± 34.71 [5.4] 
459.7 ± 29.97 [10.9] 
4.46 ± 172.53 [5.4] 
(2) 
2.85 ± 2.72 [10.9] 
1.35 ± 34.62 [5.4] 
459.2 ± 29.62 [10.9] 
3.53 ± 172.73 [5.4] 
(3) 
2.88 ± 2.75 [10.9] 
1.35 ± 34.72 [5.4] 
459.3 ± 29.84 [10.9] 
3.19 ± 171.20 [5.4] 
(4) 
2.90 ± 2.77 [11.0] 
1.33 ± 34.90 [5.5] 
457.9 ± 31.79 [11.0] 
8.93 ± 216.21 [5.5] 
All the files used for this case can be found in this rep [geo file, config file, fluid json file, solid json file].
Our first three results are quite similar to given references values. That show us that high order approximation order for space and time give us accurate values, while allow us to use less degree of freedom.
However, the lift force seems to undergo some disturbances, compared to reference results, and it’s more noticeable in our fourth case. This phenomenon is describe by [Beuer], where they’re explaining these disturbances are caused by Aitken dynamic relaxation, used in fluid structure relation for the fixed point algorithm.
In order to correct them, they propose to lower the fixed point tolerance, but this method also lowers calculation performances. An other method to solve this deviation is to use a fixed relaxation parameter \(\theta\). In this case, the optimal \(\)\theta\(\) seems to be equal to \(0.5\).
[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.
[Breuer] M. Breuer, G. De Nayer, M. Münsch, T. Gallinger, and R. Wüchner, Fluid–structure interaction using a partitioned semiimplicit predictor–corrector coupling scheme for the application of clargeeddy simulation, Journal of Fluids and Structures, 2012.
[TurekHron2] S. Turek, J. Hron, M. Madlik, M. Razzaq, H. Wobker, and JF Acker, Numerical simulation and benchmarking of a monolithic multigrid solver for fluidstructure interaction problems with application to hemodynamics, Fluid Structure Interaction II, pages 193–220, 2010.
[MunschBreuer] M. Münsch and M. Breuer, Numerical simulation of fluid–structure interaction using eddy–resolving schemes, Fluid Structure Interaction II, pages 221–253, 2010.
[Gallinger] T.G. Gallinger, Effiziente Algorithmen zur partitionierten Lösung stark gekoppelter Probleme der FluidStrukturWechselwirkung, Shaker, 2010.
[Sandboge] R. Sandboge, Fluidstructure interaction with openfsitm and md nastrantm structural solver, Ann Arbor, 1001 :48105, 2010.