1. Fluid Structure Interaction Models
The Fluid Structure models are formed from the combination of a Solid model and a Fluid model.
1.2. Fluid structure coupling conditions
In order to have a correct fluidstructure 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 fluidstruture coupling conditions, respectively velocities continuity, constraint continuity and geometric continuity.
1.2.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
1.2.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 PiolaKirchhoff tensor 
\(N.m^{2}\) 
\(\mathcal{A}_f^t\) 
Arbitrary Lagrangian Eulerian ( ALE ) map 
dimensionless 
and
Examples
2. Fluid Structure Interaction
We will interest now to the different interactions a fluid and a structure can have together with specific conditions.
2.1. Fluid structure model
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)\).
2.2. ALE
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.
3. 2D FSI elastic tube
This test case has originally been realised by [Pena], [Nobile] and [GerbeauVidrascu] only with the free outlet condition.
Computer codes, used for the acquisition of results, are from Vincent [Chabannes]
3.1. Problem Description
We interest here to the case of bidimensional blood flow modelisation. We want to reproduce and observe pressure wave spread into a canal with a fluidstructure interaction model.
The figure above shows us the initial geometry we will work on. The canal is represent by a rectangle with width and height, respectively equal to 6 and 1 cm. The upper and lower walls are mobile and so, can be moved by flow action.
By using the 1D reduced model, named generalized string and explained by [Chabannes], we didn’t need to define the elastic domain ( for the vascular wall ) here. So the structure domain is \(\Omega_s^*=\Gamma_{fsi}^*\)
During this benchmark, we will compare two different cases : the free outlet condition and the Windkessel model. The first one, as its name said, impose a free condition on the fluid at the end of the domain. The second one is used to model more realistically an flow outlet into our case. The chosen time step is \(\Delta t=0.0001\)
3.1.1. Boundary conditions
We set :

on \(\Gamma_f^{i,*}\) the pressure wave pulse \[ \boldsymbol{\sigma}_{f} \boldsymbol{n}_f = \left\{ \begin{aligned} & \left(\frac{2 \cdot 10^4}{2} \left( 1  \cos \left( \frac{ \pi t} {2.5 \cdot 10^{3}} \right) \right), 0\right)^T \quad & \text{ if } t < 0.005 \\ & \boldsymbol{0} \quad & \text{ else } \end{aligned} \right. \]

on \(\Gamma_f^{o,*}\)

Case 1 : free outlet : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f =0\)

Case 2 : Windkessel model ( \(P_0\) proximal pressure, see [Chabannes] ) : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f = P_0\boldsymbol{n}_f\)


on \(\Gamma_f^{i,*} \cup \Gamma_f^{o,*}\) a null displacement : \(\boldsymbol{\eta}_f=0\)

on \(\Gamma^*_{fsi}\) : \(\eta_s=0\)

We add also the specific coupling conditions, obtained from the axisymmetric reduced model, on \(\Gamma^*_{fsi}\)
3.2. Inputs
Name 
Description 
Nominal Value 
Units 
\(E_s\) 
Young’s modulus 
\(0.75\) 
\(dynes.cm^{2}\) 
\(\nu_s\) 
Poisson’s ratio 
\(0.5\) 
dimensionless 
\(h\) 
walls thickness 
0.1 
\(cm\) 
\(\rho_s\) 
structure density 
\(1.1\) 
\(g.cm^{3}\) 
\(R_0\) 
tube radius 
\(0.5\) 
\(cm\) 
\(G_s\) 
shear modulus 
\(10^5\) 
\(Pa\) 
\(k\) 
Timoshenko’s correction factor 
\(2.5\) 
dimensionless 
\(\gamma_v\) 
viscoelasticity parameter 
\(0.01\) 
dimensionless 
\(\mu_f\) 
viscosity 
\(0.003\) 
\(poise\) 
\(\rho_f\) 
density 
\(1\) 
\(g.cm^{3}\) 
\(R_p\) 
proximal resistance 
\(400\) 

\(R_d\) 
distal resistance 
\(6.2 \times 10^3\) 

\(C_d\) 
capacitance 
\(2.72 \times 10^{4}\) 
3.3. Outputs
After solving the fluid struture model, we obtain \((\mathcal{A}^t, \boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s)\)
with \(\mathcal{A}^t\) the ALE map, \(\boldsymbol{u}_f\) the fluid velocity, \(p_f\) the fluid pressure and \(\boldsymbol\eta_s\) the structure displacement
3.4. Discretization
\(\mathcal{F}\) is the set of all mesh faces, we denote \(\mathcal{F}_{stab}\) the face we stabilize
In fact, after a first attempt, numerical instabilities can be observed at the fluid inlet. These instabilities, caused by pressure wave, and especially by the Neumann condition, make our fluidstructure solver diverge.
To correct them, we choose to add a stabilization term, obtain from the stabilized CIP formulation ( see [Chabannes], Chapter 6 ).
As this stabilization bring an important cost with it, by increasing the number of nonnull term into the problem matrix, we only apply it at the fluid entrance, where the instabilities are located.
Now we present the different situations we worked on.
Config 
Fluid 
Structure 

\(N_{elt}\) 
\(N_{geo}\) 
\(N_{dof}\) 
\(N_{elt}\) 
\(N_{geo}\) 
\(N_{dof}\) 

\((1)\) 
\(342\) 
\(3~(P4P3)\) 
\(7377\) 
\(58\) 
\(1\) 
\(176~(P3)\) 

\((2)\) 
\(342\) 
\(4~(P5P4)\) 
\(11751\) 
\(58\) 
\(1\) 
\(234~(P4)\) 
For the fluid time discretization, BDF, at order \(2\), is the method we use.
And Newmarkbeta method is the one we choose for the structure time discretization, with parameters \(\gamma=0.5\) and \(\beta=0.25\).
These methods can be retrieved in [Chabannes] papers.
3.4.1. Solvers
Here are the different solvers ( linear and nonlinear ) used during results acquisition.
KSP 

case 
fluid 
solid 
type 
gmres 

relative tolerance 
\(1e13\) 

max iteration 
\(30\) 
\(10\) 
reuse preconditioner 
true 
false 
SNES 

case 
fluid 
solid 
relative tolerance 
\(1e8\) 

steps tolerance 
\(1e8\) 

max iteration 
\(50\) 

max iteration with reuse 
\(50\) 

reuse jacobian 
false 

reuse jacobian rebuild at first Newton step 
false 
true 
KSP in SNES 

case 
fluid 
solid 
relative tolerance 
\(1e5\) 

max iteration 
\(1000\) 

max iteration with reuse 
\(1000\) 

reuse preconditioner 
true 
false 
reuse preconditioner rebuild at first Newton step 
false 
PC 

case 
fluid 
solid 
type 
LU 

package 
mumps 
FSI 

solver method 
fix point 
tolerance 
\(1e6\) 
max iterations 
\(1\) 
3.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.
Let’s start with the main code, that can be retrieve in
feelpp/applications/models/fsi
The configuration file associated to this test is named wavepressure2d.cfg and is located at
feelpp/applications/models/fsi/wavepressure2d
The result files are then stored by default in
applications/models/fsi/wavepressure2d/P2P1G1P1G1/np_1
3.6. Results
The two following pictures have their pressure and velocity magnitude amplify by 5.
To draw the next two figures, we define 60 sections \(\{x_i\}_{i=0}^{60}\) with \(x_i=0.1i\).
All the files used for this case can be found in this rep [geo file, config file, fluid json file, solid json file].
3.6.1. Conclusion
Let’s begin with results with the free outlet condition ( see figure 2 ). These pictures show us how the pressure wave progresses into the tube. We can denote an increase of the fluid velocity at the end of the tube. Also, the wave eases at the same place. For the simulation with the Windkessel model, we observe a similar comportment at the beginning ( see figure 3 ). However, the outlet is more realistic than before. In fact, the pressure seems to propagate more naturally with this model. + In the two cases, the velocity field is disturbed at the fluidstructure interface. A mesh refinement around this region increases the quality. However, this is not crucial for the blood flow simulation.
Now we can interest us to the quantitative results.
The inflow and outflow evolution figure ( see figure 4 ) shows us similarities for the two tests at the inlet. At the outlet, in contrast, the flow increases for the free outlet condition. In fact, when the pressure wave arrived at the outlet of the tube, it is reflected to the other way. In the same way, when the reflected wave arrived at the inlet, it is reflected again. The Windkessel model reduce significantly this phenomenon. Some residues stay due to 0D coupling model and structure fixation.
We also have calculate the maximum displacement magnitude for the two model ( see figure 5 ). The same phenomenons explained ahead are retrieve here. We denote that, for the free outlet, the structure undergoes movements during the test time, caused by the wave reflection. The Windkessel model reduces these perturbations thanks to the 0D model.
The average pressure and the fluid flow ( see figure 6 and 7 ) show us the same nonphysiological phenomenons as before. The results we obtain are in accordance with the ones proposed by [Nobile].
To end this benchmark, we will compare the two resolution algorithms used with the fluidstructure model : the implicit and the semiimplicit ones. The second configuration with Windkessel model is used for the measures.
We have then the fluid flow and the displacement magnitude ( figure 8 ) curves, which superimposed on each other. So the accuracy obtained by the semiimplicit method seems good here.
The performances of the two algorithms ( figure 9 ) are expressed from number of iterations and CPU time at each step time. The semiimplicit method is a bit ahead of the implicit one on number of iterations. However, the CPU time is smaller for 2 or 3 time, due to optimization in this method. First an unique ALE map estimation is need. Furthermore, linear terms of the Jacobian matrix, residuals terms and dependent part of the ALE map can be stored and reused at each iteration.
3.7. Bibliography

[Pena] G. Pena, Spectral element approximation of the incompressible NavierStokes equations evolving in a moving domain and applications, École Polytechnique Fédérale de Lausanne, November 2009.

[Nobile] F. Nobile, Numerical approximation of fluidstructure interaction problems with application to haemodynamics, École Polytechnique Fédérale de Lausanne, Switzerland, 2001.

[GerbeauVidrascu] J.F. Gerbeau, M. Vidrascu, A quasinewton algorithm based on a reduced model for fluidstructure interaction problems in blood flows, 2003.

[Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Université de Grenoble, 2013.
4. 3D FSI elastic tube
Computer codes, used for the acquisition of results, are from Vincent [Chabannes].
4.1. Problem Description
As in the 2D case, the blood flow modelisation, by observing a pressure wave progression into a vessel, is the subjet of this benchmark. But this time, instead of a twodimensional model, we use a threedimensional model, with a cylinder
This represents the domains into the initial condition, with \(\Omega_f\) and \(\Omega_s\) respectively the fluid and the solid domain. The cylinder radius equals to \(r+\epsilon\), where \(r\) is the radius of the fluid domain and \(\epsilon\) the part of the solid domain.
\(\Gamma^*_{fsi}\) is the interface between the fluid and solid domains, whereas \(\Gamma^{e,*}_s\) is the interface between the solid domain and the exterior. \(\Gamma_f^{i,*}\) and \(\Gamma_f^{o,*}\) are respectively the inflow and the outflow of the fluid domain. Likewise, \(\Gamma_s^{i,*}\) and \(\Gamma_s^{o,*}\) are the extremities of the solid domain.
During this benchmark, we will study two different cases, named BC1 and BC2, that differ from boundary conditions. BC2 are conditions imposed to be more physiological than the ones from BC1. So we waiting for more realistics based results from BC2.
4.1.1. Boundary conditions

on \(\Gamma_f^{i,*}\) the pressure wave pulse \[ \boldsymbol{\sigma}_{f} \boldsymbol{n}_f = \left\{ \begin{aligned} & \left(\frac{1.3332 \cdot 10^4}{2} \left( 1  \cos \left( \frac{ \pi t} {1.5 \cdot 10^{3}} \right) \right), 0 \right)^T \quad & \text{ if } t < 0.003 \\ & \boldsymbol{0} \quad & \text{ else } \end{aligned} \right. \]

We add the coupling conditions on \(\Gamma^*_{fsi}\)
Then we have two different cases :

Case BC1

on \(\Gamma_f^{o,*}\) : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f =0\)

on \(\Gamma_s^{i,*} \cup \Gamma_s^{o,*}\) a null displacement : \(\boldsymbol{\eta}_s=0\)

on \(\Gamma^{e,*}_{s}\) : \(\boldsymbol{F}_s\boldsymbol{\Sigma}_s\boldsymbol{n}_s^*=0\)

on \(\Gamma_f^{i,*}\(\)U \(\)\Gamma_f^{o,*}\) : \(\mathcal{A}^t_f=\boldsymbol{\mathrm{x}}^*\)


Case BC2

on \(\Gamma_f^{o,*}\) : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f = P_0\boldsymbol{n}_f\)

on \(\Gamma_s^{i,*}\) a null displacement \(\boldsymbol{\eta}_s=0\)

on \(\Gamma^{e,*}_{s}\) : \(\boldsymbol{F}_s\boldsymbol{\Sigma}_s\boldsymbol{n}_s^* + \alpha \boldsymbol{\eta}_s=0\)

on \(\Gamma^{o,*}_{s}\) : \(\boldsymbol{F}_s\boldsymbol{\Sigma}_s\boldsymbol{n}_s^* =0\)

on \(\Gamma_f^{i,*}\) : \(\mathcal{A}^t_f=\boldsymbol{\mathrm{x}}^*\)

on \(\Gamma_f^{o,*}\) : \(\nabla \mathcal{A}^t_f \boldsymbol{n}_f^*=\boldsymbol{n}_f^*\)

4.1.2. Initial conditions
The chosen time step is \(\Delta t=0.0001\)
4.2. Inputs
Name  Description  Nominal Value  Units 

\(E_s\) 
Young’s modulus 
\(3 \times 10^6 \) 
\(dynes.cm^{2}\) 
\(\nu_s\) 
Poisson’s ratio 
\(0.3\) 
dimensionless 
\(r\) 
fluid tube radius 
0.5 
\(cm\) 
\(\epsilon\) 
solid tube radius 
0.1 
\(cm\) 
\(L\) 
tube length 
5 
\(cm\) 
\(A\) 
A coordinates 
(0,0,0) 
\(cm\) 
\(B\) 
B coordinates 
(5,0,0) 
\(cm\) 
\(\mu_f\) 
viscosity 
\(0.03\) 
\(poise\) 
\(\rho_f\) 
density 
\(1\) 
\(g.cm^{3}\) 
\(R_p\) 
proximal resistance 
\(400\) 

\(R_d\) 
distal resistance 
\(6.2 \times 10^3\) 

\(C_d\) 
capacitance 
\(2.72 \times 10^{4}\) 
4.3. Outputs
After solving the fluid struture model, we obtain \((\mathcal{A}^t, \boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s)\)
with \(\mathcal{A}^t\) the ALE map, \(\boldsymbol{u}_f\) the fluid velocity, \(p_f\) the fluid pressure and \(\boldsymbol\eta_s\) the structure displacement.
4.4. Discretization
Here are the different configurations we worked on. The parameter Incomp define if we use the incompressible constraint or not.
Config 
Fluid 
Structure 

\(N_{elt}\) 
\(N_{geo}\) 
\(N_{dof}\) 
\(N_{elt}\) 
\(N_{geo}\) 
\(N_{dof}\) 
Incomp 

\((1)\) 
\(13625\) 
\(1~(P2P1)\) 
\(69836\) 
\(12961\) 
\(1\) 
\(12876~(P1)\) 
No 

\((2)\) 
\(13625\) 
\(1~(P2P1)\) 
\(69836\) 
\(12961\) 
\(1\) 
\(81536~(P1)\) 
Yes 

\((3)\) 
\(1609\) 
\(2~(P3P2)\) 
\(30744\) 
\(3361\) 
\(2\) 
\(19878~(P2)\) 
No 
For the structure time discretization, we will use Newmarkbeta method, with parameters \(\gamma=0.5\) and \(\beta=0.25\).
And for the fluid time discretization, BDF, at order \(2\), is the method we choose.
These two methods can be found in [Chabannes] papers.
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.
Let’s start with the main code, that can be retrieve in
feelpp/applications/models/fsi
The configuration file associated to this test is named wavepressure3d.cfg and is located at
feelpp/applications/models/fsi/wavepressure3d
The result files are then stored by default in
applications/models/fsi/wavepressure3d/P2P1G1P1G1/np_1
All the files used for this case can be found in this rep [geo file, config file, fluid json file, solid json file].
5. Liddriven Cavity Flow benchmark
This test case has originally been proposed by [MokWall].
Computer codes, used for the acquisition of results, are from Vincent [Chabannes].
This benchmark has aso been realized by [Gerbeau], [Vàzquez], [Kuttler] and [Kassiotis].
5.1. Problem Description
We study here an incompressible fluid flowing into a cavity, where its walls are elastic. We use the following geometry to represent it.
The domain \(\Omega_f^*\) is define by a square \( [0,1]^2 \), \(\Gamma^{i,*}_f\) and \(\Gamma^{o,*}_f\) are respectively the flow entrance and the flow outlet. A constant flow velocity, following the \(x\) axis, will be imposed on \(\Gamma_f^{h,*}\) border, while a null flow velocity will be imposed on \(\Gamma_f^{f,*}\). This last point represent also a nonslip condition for the fluid.
Furthermore, we add a structure domain, at the bottom of the fluid one, named \(\Omega_s^*\). It is fixed by his two vertical sides \(\Gamma_s^{f,*}\), and we denote by \(\Gamma_f^{w,*}\) the border which will interact with \(\Omega_f^*\).
During this test, we will observe the displacement of a point \(A\), positioned at \((0;0.5)\), into the \(y\) direction, and compare our results to ones found in other references.
5.1.1. Boundary conditions
Before enunciate the boundary conditions, we need to describe a oscillatory velocity, following the \(\)x\(\) axis and dependent of time.
Then we can set

on \(\Gamma^{h,*}_f\), an inflow Dirichlet condition : \(\boldsymbol{u}_{f} = ( v_{in}, 0 )\)

on \(\Gamma^{i,*}_f\), an inflow Dirichlet condition : \(\boldsymbol{u}_{f} = ( v_{in}(8 y 7) ,0)\)

on \(\Gamma^{f,*}_f\), a homogeneous Dirichlet condition : \(\boldsymbol{u}_{f} = \boldsymbol{0}\)

on \(\Gamma^{o,*}_f\), a Neumann condition : \(\boldsymbol{\sigma}_{f} \boldsymbol{n}_f = \boldsymbol{0}\)

on \(\Gamma^{f,*}_s\), a condition that imposes this boundary to be fixed : \(\boldsymbol{\eta}_{s} = \boldsymbol{0}\)

on \(\Gamma^{e,*}_s\), a condition that lets these boundaries be free from constraints : \(\boldsymbol{F}_{s} \boldsymbol{\Sigma}_s \boldsymbol{n}_s = \boldsymbol{0}\)
To them we also add the fluidstructure coupling conditions on \(\Gamma_{fsi}^*\) :
5.1.2. Initial conditions
To realize the simulations, we used a time step \(\Delta t\) equals to \(0.01\) s.
5.2. Inputs
The following table displays the various fixed and variables parameters of this testcase.
Name 
Description 
Nominal Value 
Units 
\(E_s\) 
Young’s modulus 
\(250\) 
\(N.m^{2}\) 
\(\nu_s\) 
Poisson’s ratio 
\(0\) 
dimensionless 
\(\rho_s\) 
structure density 
\(500\) 
\(kg.m^{3}\) 
\(\mu_f\) 
viscosity 
\(1\times 10^{3}\) 
\(m^2.s^{1}\) 
\(\rho_f\) 
density 
\(1\) 
\(kg.m^{3}\) 
5.3. Outputs
\((\boldsymbol{u}_f, p_f, \boldsymbol{\eta}_s, \mathcal{A}_f^t)\), checking the fluidstructure system, are the output of this problem.
5.4. Discretization
To discretize space, we used \(P_N~~P_{N1}\) TaylorHood finite elements.
For the structure time discretization, Newmarkbeta method is the one we used. And for the fluid time discretization, we used BDF, at order \(q\).
5.5. Implementation
All the codes files are into FSI
5.6. Results
We begin with a \(P_2~~P_1\) approximation for the fluid with a geometry order equals at \(1\), and a fluidstructure stable interface.
Then we retry with a \(P_3~~P_2\) approximation for the fluid with a geometry order equals at \(2\), and a fluidstructure stable interface.
Finally we launch it with the same conditions as before, but with a deformed interface.
5.6.1. Conclusion
First at all, we can see that the first two tests offer us similar results, despite different orders uses. At contrary, the third result set are better than the others.
The elastic wall thinness, in the stable case, should give an important refinement on the fluid domain, and so a better fluidstructure coupling control. However, the deformed case result are closer to the stable case made measure.
5.7. Bibliography

[MokWall] DP Mok and WA Wall, Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures, Trends in computational structural mechanics, Barcelona, 2001.

[Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Université de Grenoble, 2013.

[Gerbeau] J.F. Gerbeau, M. Vidrascu, et al, A quasinewton algorithm based on a reduced model for fluidstructure interaction problems in blood flows, 2003.

[Vazquez] J.G. Valdés Vazquèz et al, Nonlinear analysis of orthotropic membrane and shell structures including fluidstructure interaction, 2007.

[KuttlerWall] U. Kuttler and W.A. Wall, Fixedpoint fluid–structure interaction solvers with dynamic relaxation, Computational Mechanics, 2008.

[Kassiotis] C. Kassiotis, A. Ibrahimbegovic, R. Niekamp, and H.G. Matthies, Nonlinear fluid–structure interaction problem ,part i : implicit partitioned algorithm, nonlinear stability proof and validation examples, Computational Mechanics, 2011.