c
TU Delft, Delft The Netherland, 2006
LEVEL SET BASED PARALLEL COMPUTATIONS OF
UNSTEADY FREE SURFACE FLOWS
E. Marchandise∗, P. Geuzaine† and J.F Remacle∗ ∗ Universit´e catholique de Louvain
Department of Civil Engineering Place du Levant 1
1348 Louvain-la-Neuve, Belgium e-mail: marchandise@gce.ucl.ac.be
† Cenaero
Avenue Jean Mermoz 30 6041 Gosselies, Belgium. e-mail:philippe.geuzaine@cenaero.be
Key words: Free surface, Parallel Computing, Level Set, Stabilised Finite Element Method
Abstract. The simulation of free surface flow is a computationally demanding task. Real life simulations require often the treatment of unstructured grids to capture the complex geometries and the use of parallel supercomputers for the efficient solution.
For that aim, we have developed an efficient parallel free surface flow solver based on a full tetrahedral discretization. An implicit stabilized finite element method is used for solv-ing the unsteady incompressible two-phase flow in three-dimensions. A pressure stabilized Petrov Galerkin technique is used to avoid spurious pressure modes, while an upwind finite volume discretization is used to discretize the advective fluxes in a stable manner. The interface between the fluid phases is captured with the level set method implemented with a quadrature-free discontinuous Galerkin method. The parallel implementation is based on the MPI message-passing standard and is fully portable.
We show the effectiveness of the method in the simulation of complex 3D flows, such as the flow past a cylindrical and the flow in a partially-filled tank of a car that suddenly brakes.
1 INTRODUCTION
The study of free surface flows covers a wide range of engineering and environmental flows, including molten metal flow, sloshing in containers, wave mechanics, open channel flows, flows around a ship or structure.
different fluids. This involves the tracking of a discontinuity in the material properties like density and viscosity.
The principal computational methods used to solve incompressible two-phase flows are the front tracking methods [10, 12], and the front capturing methods (volume of fluid [11], [24] and level set [22, 30]).
A successful approach to deal with two phase flows, especially in the presence of topo-logical changes, is the level set method (Osher and Sethian [22]). Application of level sets in two-phase flow calculations have been extensively described by Sussman, Smereka and Osher in [29–31]. The level set function is able to represent an arbitrary number of bubbles or drops interfaces and complex changes of topology are naturally taken into account by the method. The level set function φ(~x, t) is defined to be a smooth function that is positive in one region and negative in the other. The implicit surface φ(~x, t) = 0 represents the current position of an interface. This interface is advected by a velocity field u(~x, t) that is, in case of two-phase flows, the solution of the Navier-Stokes equations. The elementary advection equation for interface evolution is:
∂tφ + u· ∇φ = 0 . (1)
In [19], we have developed a high order quadrature free Runge-Kutta discontinuous Galerkin (DG) method to solve the level set equation (1) in space and time. The method was compared with classical Hamilton-Jacobi ENO/WENO methods [6, 13, 23, 29] and showed to be computationally effective and mass conservative. Besides, we showed that there was no need to reinitialize the level set.
Level sets represent a fluid interface in an implicit manner. The main advantage of this approach is that the underlying computational mesh does not conform to the interface. Hence, discontinuous integrals have to be computed in the fluid formulation because both viscosities and densities are discontinuous in all the elements crossed by the interface. The most common approach is to define a zone of thickness 2 in the vicinity of the interface (|φ| < ) and to smooth the discontinuous density and viscosity over this thickness [20,21,30,35,36]. Smoothing physical parameters in the interfacial zone may be the cause of two problems. The first one is the introduction of non physical densities and viscosities in the smoothed region, leading to possible thermodynamically aberrations [7]. The second problem is the obligation to keep the interface thickness constant in time. To ensure that the smoothed region has a constant thickness, one has to reinitialize the level set so that it remains a distance function. In this work, we rather adopt a discontinuous approach [28,34] to compute the discontinuous integrals. The use of a recursive contouring algorithm [25] allows to localize the interface accurately. Consequently, we are able to compute the discontinuous integrals with a very high level of accuracy.
A key feature of stabilized methods is that they have proved to be LBB stable and to have good convergence properties [15] [8].
In this work, we present a stabilized finite element method for computing flows in both phases [18] and combine it with a discontinuous Galerkin level set method [19] for com-puting the interface motion. The overall algorithm avoids the cost of the renormalization of the level set as well as the introduction of a non physical interface thickness and ex-hibits good mass conservation properties. This numerical technology is implemented in the parallel mode using mesh partitioning and the message passing interface (MPI).
The outline of this paper is as follows: we first present the governing equations in Section 2. Section 3 is devoted to the description of our computational method. We present the Navier-Stokes solver and the coupling with the discontinuous Galerkin method for the level set equation. In section 4, we show the effectiveness of the method in the simulation of complex 3D problems, such as the flow past a cylindrical and the flow in a partially-filled tank of a car that brakes suddenly.
2 GOVERNING EQUATIONS
We consider that the interface is represented by the iso-zero of the level set function φ. The level set function φ is positive in the liquid (denoted fluid 1) and negative in the gas (denoted fluid 2).
The governing equations describing the motion are then given by the non dimensional Navier Stokes equations:
∂tu+ u· ∇u = − ∇p
ρ(φ) + 1
ρ(φ)Re∇ · µ(φ)(∇u + ∇u
T) + g
F r2 (2)
∇ · u = 0. (3)
and by the level set equation :
∂tφ + ∇ · (uφ) = 0. (4)
From the level set function one can derive the density ρ(φ), viscosity µ(φ), the normal ~n and the curvature κ. Density and viscosity are written as
ρ(φ) = ρ2(1− H(φ)) + ρ1H(φ) (5)
µ(φ) = µ2(1− H(φ)) + µ1H(φ) (6)
where H(φ) is the Heaviside function that is equal to 1 in the liquid and 0 in the gas. Here u is the velocity field, p the pressure field, ρ and µ are the density and dynamic viscosity of the fluid all in non-dimensional forms; these are variables in whole domain but constant and in general different in each phase. g denotes the gravitational field.
The key flow parameters are then the ratios of density and viscosity (ρ1/ρ2, µ1/µ2),
the Reynolds number Re = ρRURLR/µR , the Froude number F r = UR/√gLR . The
3 COMPUTATIONAL METHOD
To derive the finite element discretization, we assume that we have some appropriate finite dimensional function spaces for the trial (Sh
u,S h
p,Sφh) and weighting (Vuh,V h p,Vφh)
function spaces corresponding to the velocity, pressure and the level set function.
Concerning the function spaces, we use linear continuous approximations for the ve-locity and pressure and piecewise continuous approximations of order p on each element for the level set. The physical domain Ω of boundary Γ is discretized into a collection of Ne elements (Ωe) called a mesh. The same computational mesh is used for the Navier
Stokes solver and the Interface solver. 3.1 The Navier-Stokes solver
We employ a pressure stabilized pressure finite element (PSPG) [32] method for the discretization of the Navier Stokes equations. For the stabilization of the upwind term, we follow the work of Barth and Selmin [2] and use a finite volume stabilization with the finites volumes being the median cells of the mesh .
The stabilized finite element formulation of Eqs. (2) and (3) can be written as follows: find uh ∈ Sh
u and p h ∈ Sh
p such that ∀vh ∈ Vuh and ∀q h ∈ Vh p : R Ωev h· (uh t + uh· ∇uh+ ∇ph ρ − g F r2) dv + R Ωe ∇v h· µ
ρRe(∇u + ∇u
T) dv (7) +R Ωeq h∇ · uh dv + ST = R Γv h· hh ds
where hh denotes the Neumann-type boundary condition associated with the momentum
equations hh = σ· n on the interface Γ, where σ = −pI + µ(∇u + ∇uT) is the stress
tensor.
The stabilization term ST
ST = X
e
τe∇qh, R(ph, uh)
Ωe (8)
contains the residual of the momentum equation R(ph, uh) = uh
t + (u h
· ∇)uh
− 1
ρRe∇ · µ(∇u + ∇u
T) + ∇ph
ρ −
g
F r2. (9)
The stabilization parameter τ is of orderO (h2e/ν) in the diffusion dominated case and
of order of O (he/kuk) in the advection dominated case [32, 33].
3.1.1 Discontinuous density and viscosity
As the density and the viscosity are constant in each phase but may be discontinu-ous across the interface (the interface does not conform with the mesh), some integrals of Eq. (7) are discontinuous. Those discontinuous integrals are computed exactly by numerical integration: Z Ωe vh ∇p h ρ dv = Z V1 vh ∇p h ρ1 dv + Z V2 vh ∇p h ρ2 dv (10) Z Ωe ∇vh µ ρRe∇u h dv = Z V1 ∇vh µ1 ρ1Re∇u h dv + Z V2 ∇vh µ2 ρ2Re∇u h dv (11)
where the subscripts 1 and 2 denote the fluid 1 and the fluid 2. 3.2 The Interface solver
We employ a quadrature free discontinuous Galerkin level set method [19] for the discretization of the level set equation.
The discontinuous finite element formulation of Eqs.˜refeq:levelset2 can be written as follows: find φh ∈ Sh φ such that ∀ ˆφh ∈ Vφh: Z Ωe ˆ φh· ∂tφ dv + Z Ωe ∇ ˆφh· (uh φ) dv− Z ∂Ωe ˆ φh· f ds (12)
where f = (φu· n) is the normal trace of the fluxes and is chosen to be the upwind flux. The resulting system is solved in time using an explicit Runge-Kutta method.
More details of implementation of the discontinuous Galerkin method as well as super-convergence analysis can be found in [19].
3.3 Summary of the algorithm
Our computational approach for numerical modeling of two-phase flows can be sum-marized as follows. First we initialize all the variables (velocity u, pressure p, level set φ). Then,
For each time tn, n = 0, 1, 2...:
• First advance implicitly the Navier-Stokes equations from tn−1 until tn and find a
velocity un. This equation is solved implicitly for the velocity and pressure variables
but uses the previous position of the level set φn−1 to identify the fluid variable such as viscosity and density. To go to step 2, project the velocity field u onto the dofs of the interface solver.
• Second advance explicitly the level set (Eq. (4)) from tn−1 until tn (using several
smaller time steps) with a linear velocity varying in time between un−1 until un and
find a new position for the level set φn. Then project the level set function φ onto
• Third go back to step 2, with this new position φn and repeat those steps until
kφn−1− φnk < eps.
4 NUMERICAL RESULTS
In this section, we show the effectiveness of the method in the simulation of complex 3D problems, such as the flow past a cylindrical and the sloshing in a tanker during braking. Our methods were implemented in Standard C++, compiled with INTEL icc v8.0, and run on the NEC/Clustervision Linux cluster installed at CENAERO.
The iterations are carried out with a GMRES update technique with a Krylov space of 60.
4.1 3D Dam break with a cylindrical obstacle
To show the ability to simulate 3D free surface flows, we consider the breaking of a cubic water column in a domain containing a cylindrical obstacle .
The computational domain is described by the length L = 6 m and the height H = 1.5 m . The height and the width of the water column are hl = 1 m and hw = 1.5 m. The
cylinder is located 1.3 m downstream the water column and has a diameter of 0.4 m. The mesh depicted in fig. 1 is unstructured and made of 30, 820 nodes.
Figure 1: Computational mesh for the dam break problem made of 30.820 nodes.
Figure 2 shows snapshots of the water surface position at selected times. We see that from the time the water front reaches the cylinder, the flow shows clearly three dimensionality.
Figure 2: Perspective view of the free surface for the dam break with a cylindrical obstacle at the selected times t = 1s and t = 2s. calculated by: V = V (tf)− V (0) V (0) (13)
where V (tf) is the total volume of the liquid at the final time tf.
Table 1: Mass fluctuation of the method for the 3D dam break with cylindrical obstacle computed at time t = 1s and t = 2s on different computational meshes and with different polynomial orders p to approximate the level set
Mesh order Mass loss Mass loss #CP U Elapsed Time
# nodes p (t = 1s) (t = 2s) 12,500 1 7.4 % 14.3 % 2 1h30 12,500 2 3.4% 6.2 % 2 3h10 12,500 3 2.1% 3.1 % 4 15h10 30,820 2 1.16% 0.86 % 8 13h10 74,016 2 0.62% 0.11 % 24 6h45
mass conservation. We show that those results are obtained within reasonable computa-tional time. The computacomputa-tional time may however be dramatically reduced by performing mesh adaptation near the interface. We are currently working on this topic.
4.2 Sloshing in a partially filled tank that brakes suddenly.
Sloshing is known as the oscillation of the free surface of a liquid in an externally excited container. Sloshing of liquids in moving containers is of practical concern in many engineering applications, such as fuel sloshing in aircraft tanks as a result of sharp flight maneuvers, liquid sloshing in large storage tanks due to earthquakes and fluid oscillation in tanker trucks traveling on highways.
Partially filled liquid tanks undergoing accelerated motions are susceptible to slosh-induced loads, which may affect the dynamic behavior of the liquid-retaining structure, and may even be severe enough to cause structural damage.
Accurate numerical simulation of the sloshing phenomena is very important in the design of the tank. It allows fuel system validation before actual prototyping and indus-trialization and therefore, reducing time and costs for their development.
Large amplitude liquid sloshing in partially filled tank vehicles has been studied in [1, 16, 27].
The geometry of the tank considered here is a simplified geometry of a vehicle tank. The computational mesh (fig. 3 is made of 605, 923 tetrahedral elements and 106, 384 nodes).
It is assumed that the vehicle is moving at a constant velocity and it suddenly brakes with a constant deceleration equal to 5 m/s2 for 0.5 seconds. Afterwards, it moves again
with a constant velocity.
Since the inertial force due to the acceleration of the tank is very large compared to the viscous forces in the fluids, a slip-wall boundary condition was imposed on the tank walls.
The level of water is 12 cm measured from the bottom of the fuel tank. This represents about 45 liters of fuel.
Figure 3: Computational mesh for the vehicle tank made of 106, 384 nodes (left) and cut in the mesh (right).
The computation was carried out with 32 processors and took about 19h to simulate 2s.
5 CONCLUSIONS
A unified approach for the numerical simulation of three dimensional two-phase flows has been presented. The approach relies on an implicit stabilized finite element approxi-mation for the Navier-Stokes equations and discontinuous Galerkin method for the level set method.
Such a combination of those two numerical methods results in a simple and effective al-gorithm that allows to simulate diverse flow regimes presenting large density and viscosity ratios (ratio of 1000).
Two advantages of the method are:
• Simplicity and flexibility: The stabilized method utilizes linear elements for the unknowns of velocity and pressure and a discontinuous Galerkin method of higher order for the level set unknown. The tetrahedral meshes can be structured or unstructured;
• Accuracy: The overall scheme is second order accurate in space and time. Besides we do not need to introduce any artificial parameter as interface thickness. The interface can be localized precisely which enables us to compute accurately the integrals with a discontinuous density and/or viscosity. Furthermore, the method exhibits excellent conservation properties using high order polynomials for the level set.
X Y Z X Y Z
Figure 4: Sloshing during a vehicle tank during braking. The images show at different instants (t = 0.1 s, 0.2 s, 0.4 s) , the free surface position (left in red) and the pressure distribution (right).
REFERENCES
[1] S. Aliabadi and T.E. Tezduyar. Stabilized-finite-element/interface-capturing tech-nique for parallel computation of unsteady flows with interfaces. Computer Methods in Applied Mechanics and Engineering, 190:243–261, 2000.
[3] X. C. Cai, C. Farhat, and M. Sarkis. A minimum overlap restricted additive Schwarz preconditioner and applications in 3D flow simulations. Contemporary Mathematics, 218:478–484, 1998.
[4] A.J Chorin. Numerical solution of the Navier-Stokes equations. Math Comput, 22:745–762, 1968.
[5] A.J Chorin. A numerical method for solving viscous floxw problems. Journal of Computational Physics, 135:118–125, 1997.
[6] D. Enright, R. Fedkiw, J. Ferziger, and I. Mitchell. A Hybrid Particle Level Set Method for Improved Interface Capturing. Journal of Computational Physics, 183:83–116, 2002.
[7] J.H. Ferziger. Interfacial transfer in Tryggvason’s method. International Journal for Numerical Methods in Fluids, 41(5):551–560, 2003.
[8] L.P. Franca and S. Frey. Stabilized finite element methods: Ii. the incompressible navier stokes equations. Computer Methods in Applied Mechanics and Engineering, 59:209–233, 1992.
[9] P. Geuzaine. Newton-Krylov strategy for Compressible Turbulent Flows on Unstruc-tured Meshes. AIAA Journal, 39:528–531, 2001.
[10] F.H. Harlow and J.E. Welch. Volume tracking methods for interfacial flow calcula-tions. Physics of fluids, 8:21–82, 1965.
[11] C Hirt and B. Nichols. Volume of Fluid Method (VOF) for the dynamics of free boundaries. Journal of Computational Physics, 39:201–225, 1981.
[12] C. W. Hirt, A. A. Amsden, and J. L. Cook. An arbitrary Lagrangian-Eulerian computing method for all flow speeds. J. Comput. Phys., 135(2), 1997.
[13] C Hu and C.-W. Shu. Weighted ENO Schemes for Hamilton–Jacobi Equations. SIAM Journal on Scientific Computing, 21(6):2126–2143, 1999.
[14] T. J. R. Hughes and L. P. Franca. A new finite element formulation for computational fluid dynamics: V. circumventing the Babuska–Brezzi condition: A stable Petrov– Galerkin formulation of the stokes problem accomodating equal-order interpolations. Computer Methods in Applied Mechanics and Engineering, 59:85–99, 1986.
[16] R.A. Ibrahim. Liquid Sloshing Dynamics : Theory and Applications. Cambridge University Press, 2005.
[17] Y. Kallinderis and K. Nakajima. Finite Element Method for incompressible viscous flow with adaptative grids. AIAA Journal, 8(32), 1994.
[18] E. Marchandise and J.-F. Remacle. A stabilized finite element method using a dis-continuous level set approach for solving two phase incompressible flows. Journal of Computational Physics, 2006.
[19] E. Marchandise, J.-F. Remacle, and Nicolas Chevaugeon. A Quadrature free dis-continuous Galerkin method for the level set equation. Journal of Computational Physics, 212:338–357, 2006.
[20] S Nagrath and K.E. Jansen. Three dimensional simulation of incompressible two-phase flows using the finite element method and level set approach. Computer Meth-ods in Applied Mechanics and Engineering, page in press, 2004.
[21] S. Nagrath, K.E. Jansen, and R.T. Lahey. Computation of incompressible bubble dynamics with a stabilized finite element level set method. Computer Methods in Applied Mechanics and Engineering, 194:4565–4587, 2005.
[22] S Osher and J. A. Sethian. Fronts propagating with curvature dependent speed: al-gorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79:12–49, 1988.
[23] S Osher and C.-W. Shu. High order essentially non-oscillatory shemes for Hamilton-Jacobi equations. SIAM Journal on Numerical Analysis, 728:902–921, 1991.
[24] J.E Pilliot and E.G. Puckett. Second order accurate Volume-of-Fluid algorithms for tracking material interfaces. Journal of Computational Physics, 199:465–502, 2004. [25] J.-F. Remacle, E. Marchandise, N. Chevaugeon, and C. Geuzaine. Efficient
visualiza-tion of high order finite elements. SIAM Journal of Scientific Computing, Accepted, 2006.
[26] S. Rogers, D. Kwak, and C. Kiris. Steady and unsteady solutions of the incompress-ible Navier-Sstokes equations. AIAA Journal, 29:306–610, 1991.
[27] W. Rumold. Modeling and simulation of vehicles carrying liquid cargo. Multibody System Dynamics, 5(4):351–374, 2001.
[29] M. Sussman, E. Fatemi, and P. Smereka. An improved level set method for incom-pressible two-fluid flows. Computers and Fluids, 27:663–680, 1998.
[30] M. Sussman, E. Fatemi, P. Smereka, and S. Osher. An efficient interface preserving level set re-distancing algorithm and its application to interfacial incompressible fluid flow. SIAM Journal on Scientific Computing, 20:1165–1191, 1999.
[31] M. Sussman, P. Smereka, and S. P. Osher. A level set approach for computing solu-tions to incompressible two-phase flows. Journal of Computational Physics, 114:146– 159, 1994.
[32] S.E. Ray T.E. Tezduyar, S. Mittal and R. Shih. Incompressible flow computa-tions with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements. Computer Methods in Applied Mechanics and Engineering, 95:221–242, 1992.
[33] T.E. Tezduyar and Y. Osawa. Finite element stabilization parameters computed from element matrices and vectors. Computer Methods in Applied Mechanics and Engineering, 190:411–430, 2000.
[34] A.-K Tornberg. Interface Tracking Methods with Applications to Multiphase Flows. PhD thesis, Royal Institute of Technology (KTH), Finland, 2000.
[35] A.-K. Tornberg and B. Engquist. A finite element based level set method for multi-phase flow applications. Computing Visualization in Science, 3:93–101, 2000.
