# Turek & Hron CFD Benchmark

## Introduction

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 [img-geometry1].

The software and the numerical results were initially obtained from Vincent Chabannes.

 This benchmark is linked to the Turek Hron CSM and Turek Hron FSI benchmarks.

## Problem Description

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.

Figure 1. Geometry of the Turek & Hron CFD Benchmark

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.

### Boundary conditions

We set

• on \Gamma_{in}, an inflow Dirichlet condition :  \boldsymbol{u}_f=(v_{in},0)

• on \Gamma_{wall} and \Gamma_{obst}, a homogeneous Dirichlet condition :  \boldsymbol{u}_f=\boldsymbol{0}

• on \Gamma_{out}, a Neumann condition :  \boldsymbol{\sigma}_f\boldsymbol{ n }_f=\boldsymbol{0}

### Initial conditions

We use a parabolic velocity profile, in order to describe the flow inlet by  \Gamma_{in} , which can be express 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, we want to impose a progressive increase of this velocity profile. That’s why 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.

Moreover, in this case, there is no 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

\nu_f

kinematic viscosity

1\times 10^{-3}

m^2/s

\mu_f

dynamic viscosity

1

kg/(m \times s)

\rho_f

density

1000

kg/m^3

f_f

source term

0

kg/(m^3 \times s)

\bar{U}

characteristic inflow velocity

CFD1 CFD2 CFD3

0.2

1

2

m/s

## Outputs

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

(F_D,F_L)=\int_{\Gamma_{obst}}\boldsymbol{\sigma}_f \boldsymbol{ n }_f

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

## Discretization

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.

### Solvers

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

 type gmres relative tolerance 1e-13 max iteration 1000 reuse preconditioner false
 relative tolerance 1e-8 steps tolerance 1e-8 max iteration CFD1/CFD2 : 100 | CFD3 : 50 max iteration with reuse CFD1/CFD2 : 100 | CFD3 : 50 reuse jacobian false reuse jacobian rebuild at first Newton step true
 relative tolerance 1e-5 max iteration 1000 max iteration with reuse CFD1/CFD2 : 100 | CFD3 : 1000 reuse preconditioner false reuse preconditioner rebuild at first Newton step false
 type lu package mumps

## Running the model

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

Command line to execute CFD1 testcase
``\$ mpirun -np 4 /usr/local/bin/feelpp_toolbox_fluid_2d --config-file cfd1.cfg``

The result files are then stored by default in

Results Directory
``feel/applications/models/fluid/TurekHron/"case_name"/"velocity_space""pression_space""Geometric_order"/"processor_used"``

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

``feel/toolboxes/fluid/TurekHron/cfd2/P2P1G1/np_12``

## Results

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

### CFD1

Table 6. Results for CFD1
\mathbf{N_{geo}} \mathbf{N_{elt}} \mathbf{N_{dof}} Drag Lift

Reference Turek and Hron

14.29

1.119

1

9874

45533 (P_2/P_1)

14.217

1.116

1

38094

173608 (P_2/P_1)

14.253

1.120

1

59586

270867 (P_2/P_1)

14.262

1.119

2

7026

78758 (P_3/P_2)

14.263

1.121

2

59650

660518 (P_3/P_2)

14.278

1.119

3

7026

146057 (P_4/P_3)

14.270

1.120

3

59650

1228831 (P_4/P_3)

14.280

1.119

All the files used for this case can be found in this rep [geo file, config file, json file]

### CFD2

Table 7. Results for CFD2
\mathbf{N_{geo}} \mathbf{N_{elt}} \mathbf{N_{dof}} Drag Lift

Reference Turek and Hron

136.7

10.53

1

7020

32510 (P_2/P_1)

135.33

10.364

1

38094

173608 (P_2/P_1)

136.39

10.537

1

59586

270867 (P_2/P_1)

136.49

10.531

2

7026

78758 (P_3/P_2)

136.67

10.548

2

59650

660518 (P_3/P_2)

136.66

10.532

3

7026

146057 (P_4/P_3)

136.65

10.539

3

59650

1228831 (P_4/P_3)

136.66

10.533

All the files used for this case can be found in this rep [geo file, config file, json file]

### CFD3

As CFD3 is time-dependent ( from BDF use ), results will be expressed as

mean ± amplitude [frequency]

where

• mean is the average of the min and max values at the last period of oscillations.

mean=\frac{1}{2}(max+min)

• amplitude is the difference of the max and the min at the last oscillation.

amplitude=\frac{1}{2}(max-min)

• frequency can be obtain by Fourier analysis on periodic data and retrieve the lowest frequency or by the following formula, if we know the period time T.

frequency=\frac{1}{T}

Table 8. Results for CFD3
\mathbf{\Delta t} \mathbf{N_{geo}} \mathbf{N_{elt}} \mathbf{N_{dof}} \mathbf{N_{bdf}} Drag Lift

0.005

Reference Turek and Hron

439.45 ± 5.6183[4.3956]

−11.893 ± 437.81[4.3956]

 0.01 1 8042 37514 (P_2/P_1) 2 437.47 ± 5.3750[4.3457] -9.786 ± 437.54[4.3457] 2 2334 26706 (P_3/P_2) 2 439.27 ± 5.1620[4.3457] -8.887 ± 429.06[4.3457] 2 7970 89790 (P_2/P_2) 2 439.56 ± 5.2335[4.3457] -11.719 ± 425.81[4.3457]
 0.005 1 3509 39843(P_3/P_2) 2 438.24 ± 5.5375[4.3945] -11.024 ± 433.90[4.3945] 1 8042 90582 (P_3/P_2) 2 439.25 ± 5.6130[4.3945] -10.988 ± 437.70[4.3945] 2 2334 26706 (P_3/P_2) 2 439.49 ± 5.5985[4.3945] -10.534 ± 441.02[4.3945] 2 7970 89790 (P_3/P_2) 2 439.71 ± 5.6410[4.3945] -11.375 ± 438.37[4.3945] 3 3499 73440 (P_4/P_3) 3 439.93 ± 5.8072[4.3945] -14.511 ± 440.96[4.3945] 4 2314 78168 (P_5/P_4) 2 439.66 ± 5.6412[4.3945] -11.329 ± 438.93[4.3945]
 0.002 2 7942 89482 (P_3/P_2) 2 439.81 ± 5.7370[4.3945] -13.730 ± 439.30[4.3945] 3 2340 49389 (P_4/P_3) 2 440.03 ± 5.7321[4.3945] -13.250 ± 439.64[4.3945] 3 2334 49266 (P_4/P_3) 3 440.06 ± 5.7773[4.3945] -14.092 ± 440.07[4.3945]

All the files used for this case can be found in this rep [geo file, config file, json file].

Figure 2. Lift and drag forces

## Geometrical Order

 Add a section on geometrical order.

## Conclusion

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.

## 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ée de Grenoble, 2013.

# 2D Drops Benchmark

This benchmark has been proposed and realised by Hysing. It allows us to verify our level set code, our Navier-Stokes solver and how they couple together.

Computer codes, used for the acquisition of results, are from Vincent Doyeux, with the use of Chabannes's Navier-Stokes code.

## Problem Description

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 :

• The gravity force : \boldsymbol{f}_g=\rho_\phi\boldsymbol{g}

• The surface tension force : \boldsymbol{f}_{st}=\int_\Gamma\sigma\kappa\boldsymbol{ n }

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.

Figure 3. Geometry used in 2D Bubble Benchmark

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

### Boundary conditions

• On the lateral walls, we imposed slip conditions

\begin{eqnarray} \boldsymbol{u}\cdot\boldsymbol{n}&=&0 \\ t\cdot(\nabla\boldsymbol{u}+^t\nabla\boldsymbol{u})\cdot \boldsymbol{n}&=&0 \end{eqnarray}
• On the horizontal walls, no slip conditions are imposed : \boldsymbol{u}=0

### Initial conditions

In order to let the bubble rise, its density must be inferior to the density of the exterior fluid, so \rho_1>\rho_2

## Inputs

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

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

\boldsymbol{g}

gravity acceleration

(0,0.98)

m/s^2

l

length domain

1

m

h

height domain

2

m

r

0.25

m

B_c

bubble center

(0.5,0.5)

m

## Outputs

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

\boldsymbol{X}_c = \dfrac{ \displaystyle \int_{\Omega_2} \boldsymbol{x}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{ \displaystyle \int_\Omega \boldsymbol{x} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }
\boldsymbol{U}_c = \dfrac{\displaystyle \int_{\Omega_2} \boldsymbol{u}}{ \displaystyle \int_{\Omega_2} 1 } = \dfrac{\displaystyle \int_\Omega \boldsymbol{u} (1-H_\varepsilon(\phi))}{ \displaystyle \int_\Omega (1-H_\varepsilon(\phi)) }
c = \dfrac{\left(4 \pi \displaystyle \int_{\Omega_2} 1 \right)^{\frac{1}{2}}}{ \displaystyle \int_{\Gamma} 1} = \dfrac{ \left(4 \pi \displaystyle \int_{\Omega} (1 - H_\varepsilon(\phi)) \right) ^{\frac{1}{2}}}{ \displaystyle \int_{\Omega} \delta_\varepsilon(\phi)}

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.

## Discretization

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

 Case \rho_1 \rho_2 \mu_1 \mu_2 \sigma Re E_0 ellipsoidal bubble (1) 1000 100 10 1 24.5 35 10 squirted bubble (2) 1000 1 10 0.1 1.96 35 125

## Results

### Test 2

We describe the different quantitative results for the two studied cases.

 h c_{min} t_{c_{min}} u_{c_{max}} t_{u_{c_{max}}} y_c(t=3) lower bound 0.9011 1.8750 0.2417 0.9213 1.0799 upper bound 0.9013 1.9041 0.2421 0.9313 1.0817 0.02 0.8981 1.925 0.2400 0.9280 1.0787 0.01 0.8999 1.9 0.2410 0.9252 1.0812 0.00875 0.89998 1.9 0.2410 0.9259 1.0814 0.0075 0.9001 1.9 0.2412 0.9251 1.0812 0.00625 0.8981 1.9 0.2412 0.9248 1.0815
 h c_{min} t_{c_{min}} u_{c_{max_1}} t_{u_{c_{max_1}}} u_{c_{max_2}} t_{u_{c_{max_2}}} y_c(t=3) lower bound 0.4647 2.4004 0.2502 0.7281 0.2393 1.9844 1.1249 upper bound 0.5869 3.0000 0.2524 0.7332 0.2440 2.0705 1.1380 0.02 0.4744 2.995 0.2464 0.7529 0.2207 1.8319 1.0810 0.01 0.4642 2.995 0.2493 0.7559 0.2315 1.8522 1.1012 0.00875 0.4629 2.995 0.2494 0.7565 0.2324 1.8622 1.1047 0.0075 0.4646 2.995 0.2495 0.7574 0.2333 1.8739 1.1111 0.00625 0.4616 2.995 0.2496 0.7574 0.2341 1.8828 1.1186

## Bibliography

References for this benchmark
• [Hysing] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska, Quantitative benchmark computations of two-dimensional bubble dynamics, International Journal for Numerical Methods in Fluids, 2009.

• [Chabannes] V. Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles. PhD thesis, Université de Grenoble, 2013.

• [Doyeux] V. Doyeux, Modélisation et simulation de systèmes multi-fluides, Application aux écoulements sanguins, PhD thesis, Université de Grenoble, 2014.

# 3D Drop benchmark

The previous section described the strategy we used to track the interface. We couple it now to the Navier Stokes equation solver described in \cite{chabannes11:_high}. In the current section, we present a 3D extension of the 2D benchmark introduced in \cite{Hysing2009} and realised using Feel++ in \cite{Doyeux2012}.

## Benchmark problem

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:

\begin{array}[lll] \rho\rho(\phi(\mathbf{x}) ) \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) + \nabla p - \nabla \cdot \left( \nu(\phi(\mathbf{x})) (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) \right) &=& \rho ( \phi(\mathbf{x}) ) \mathbf{g}, \\ \nabla \cdot \mathbf{u} &=& 0, \\ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi &=& 0, \end{array}

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

 Tests \rho_1 \rho_2 \nu_1 \nu_2 \sigma Re E_0 Test 1 (ellipsoidal bubble) 1000 100 10 1 24.5 35 10 Test 2 (skirted bubble) 1000 1 10 0.1 1.96 35 125

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

## Simulations parameters

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.

 h Number of processors \Delta t Time per iteration (s) Total Time (h) 0.025 360 0.0125 18.7 1.25 0.02 360 0.01 36.1 3.0 0.0175 180 0.00875 93.5 8.9 0.015 180 0.0075 163.1 18.4 0.0125 180 0.00625 339.7 45.3
 h Tetrahedra Points Order 1 Order 2 0.025 73010 14846 67770 1522578 0.02 121919 23291 128969 2928813 0.0175 154646 30338 187526 4468382 0.015 217344 41353 292548 6714918 0.0125 333527 59597 494484 11416557

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.

## Results Test 1: Ellipsoidal bubble

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.

## Bibliography

.

References for this benchmark
• [cottet]

• [Feelpp] C. Prud’homme et al.

• [osher] Osher1988, book_Sethian, book_Osher

• [Franca1992] Franca 1992