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 [Chabannes]. In the current section, we present a 3D extension of the 2D benchmark introduced in [Hysing] and realised using Feel++ in [Doyeux].

1. 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 \} 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ötvös 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 [Hysing].

Table 1. Numerical parameters taken for the benchmarks.

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 [Hysing] 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}\).

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

Table 2. Numerical parameters used for the test 1 simulations: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation.

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

Table 3. 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

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.

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

4. Bibliography

.

References for this benchmark
  • [cottet]

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

  • [osher] Osher1988, book_Sethian, book_Osher

  • [Franca1992] Franca 1992

  • [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.

  • [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.