Turek Hron FSI Benchmark

In order to validate our fluid-structure 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.

Problem Description

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.

Figure 1 : Geometry of the Turek & Hron FSI Benchmark.

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 Navier-Stockes 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 Saint-Venant-Kirchhoff 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 fluid-structure coupling system.

Boundary conditions

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}^*$ :

$\frac{\partial \boldsymbol{\eta_{s}} }{\partial t} - \boldsymbol{u}_f \circ \mathcal{A}^t_f = \boldsymbol{0} \quad \left( \text{Velocities continuity}\right)$
$\boldsymbol{F}_{s} \boldsymbol{\Sigma}_{s} \boldsymbol{n}^*_s + J_{\mathcal{A}^t_f} \boldsymbol{F}_{\mathcal{A}^t_f}^{-T} \hat{\boldsymbol{\sigma}}_f \boldsymbol{n}^*_f = \boldsymbol{0} \quad \left( \text{ Constraint continuity}\right)$
$\boldsymbol{\varphi}_s^t - \mathcal{A}^t_f = \boldsymbol{0} \quad \left( \text{Geometric continuity}\right)$

where $\boldsymbol{n}_f$ is the outer unit normal vector from $\partial \Omega_f^*$.

Initial conditions

Fluid

In order to describe the flow inlet by $\Gamma_{in}$, a parabolic velocity profile is used. It can be expressed by

$v_{cst} = 1.5 \bar{U} \frac{4}{0.1681} y \left(0.41-y\right)$

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

Inputs

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

Table 1. Fixed and Variable Input Parameters
Name Description Nominal Value Units

$l$

elastic structure length

$0.35$

$m$

$h$

elastic structure height

$0.02$

$m$

$r$

$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 non-linear part and a direct method based on LU decomposition is selected for the linear part.

Outputs

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.

This system also give us the ALE map $\mathcal{A}_f^t$.

Discretization

To realize these tests, we made the choice to used $P_N~-~P_{N-1}$ Taylor-Hood finite elements to discretize the space.

For the fluid time discretization, BDF, at order $q$, is the method we have chosen.

And finally Newmark-beta 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.

Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

 KSP case fluid solid type gmres relative tolerance $1e-13$ max iteration $1000$ reuse preconditioner true
 SNES case fluid solid relative tolerance $1e-8$ steps tolerance $1e-8$ 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 $1e-5$ 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 $1e-6$ max iterations $1000$ initial $\theta$ $0.98$ minimum $\theta$ $1e-12$

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.

`    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/P2P1G1-P1G1/np_4`

Results

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^{N-1}_c(\Omega_{f,\delta}) \times V^{N-1}_{s,\delta}$ $\Delta t$ Turek and Hron 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 [TurekHron] -2.69 ± 2.53 [10.9] 1.48 ± 34.38 [5.3] 457.3 ± 22.66 [10.9] 2.22 ± 149.78 [5.3] [Breuer] 464.5 ± 40.50 6.00 ± 166.00 [5.5] [TurekHron2] -2.88 ± 2.72 [10.9] 1.47 ± 34.99 [5.5] 460.5 ± 27.74 [10.9] 2.50 ± 153.91 [5.5] [MunschBreuer] -4.54 ± 4.34 [10.1] 1.50 ± 42.50 [5.1] 467.5 ± 39.50 [10.1] 16.2 ± 188.70 [5.1] [Gallinger] 474.9 ± 28.10 3.90 ± 165.90 [5.5] [Sandboge] -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].

Conclusion

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

Bibliography

References for this benchmark
• [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure 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. Münsch, T. Gallinger, and R. Wüchner, Fluid–structure interaction using a partitioned semi-implicit predictor–corrector coupling scheme for the application of clarge-eddy 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 fluid-structure 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 Lösung stark gekoppelter Probleme der Fluid-Struktur-Wechselwirkung, Shaker, 2010.

• [Sandboge] R. Sandboge, Fluid-structure interaction with openfsitm and md nastrantm structural solver, Ann Arbor, 1001 :48105, 2010.