On numerical approximation of incompresible fluid flow with free surface influenced by surface tension

This paper deals with the numerical approximation of two immiscible incompressible fluids flow. The two-dimensional problem is mathematically described with the aid of the incompressible Navier-Stokes equations. The motion of the two phases is treated with the aid of level set method. The surface tension effects are taken into account as well as the contact angle influence. The problem is formulated in a weak sense and discretized with the aid of the finite element method. The details of the numerical methods are given and several possibilities of the realization of the surface tension term are discussed. Numerical results are shown.


Introduction
The numerical approximation of two-phase flows with free surface is very important in various scientific as well as in technical applications or industrial processes, see e.g. [1], [2], [3]. The kinematics and dynamics of the fluid interface plays an important role in many of these processes. The approximation of free surface flow influenced by the high surface tension and possibly contact angles is very difficult problem, see e.g. [4], [5], [6], [7], [8].
In order to approximate free surface flows there are several fast, robust and accurate numerical methods available dealing with fixed computing grids. Particularly, these methods are the marker and cell method [9], volume of fluid method [10] and level set method [11], [12], [13]. There are many papers describing the application of the mentioned methods, but its reliability is usually not tested.
In this paper, the two-dimensional flow of two immiscible fluids is considered, mathematically described. In order to take into account the surface tension effects the variational formulation is introduced and the discretization with the the finite element method is applied. The motion of the interface between these two fluids (free surface) is realized using the level set method, cf. [11]. A modification of the standard finite element method is used in order to properly address the arising spurious currents nearby the free surface caused by the discontinuous pressure, see [7] or [5]. The implemented method is verified by solution of a benchmark problem, cf. [4].

Mathematical model
In this section two immiscible fluids are assumed to occupy a domain Ω ⊂ R 2 . The domain Ω is a bounded coma Corresponding author: petr.svacek@fs.cvut.cz putational domain with the Lipschitz continuous boundary decomposed into three mutually disjoint parts Γ W (representing fixed wall) , Γ S (representing part of boundary where the slip boundary condition is prescribed) and Γ O (outlet part of the boundary). Each fluid is assumed to occupy a time dependent subdomain (Ω 1 (t) or Ω 2 (t) ) with the moving interfaceΓ t (free surface), i.e. the interface is the common boundary of both subdomainŝ Γ t = ∂Ω 1 (t) ∩ ∂Ω 2 (t) . Both incompressible fluids are characterized by the density (ρ 1 and ρ 2 ) and its viscosity (μ 1 and μ 2 ).
In order to describe the mathematical model in a single equation, the Heaviside function H(x, t) is introduced Further, the density ρ(x, t) and viscosity μ(x, t) functions are defined using the Heaviside functions as functions giving either the fluid local density or viscosity at the point x ∈ Ω at time t ∈ (0, T ). Thus the density function reads and the viscosity function is given by formula Similarly, the pressure and the flow velocity functions defined on the domain Ω are denoted by p(x, t) and u = u(x, t), respectively. The flow motion in Ω i (t) , i = 1, 2 is then described by the incompressible Navier-Stokes equations ρ ∂u ∂t + ρ(u · ∇)u − ∇ · σ = ρf . On the interfaceΓ t the following boundary conditions including the surface tension effects needs to be specified where γ is the surface tension coefficient, κ denotes the curvature of the interface Γ I and n here denotes the normal to the Γ I pointing into Ω 2 (t) . The Navier-Stokes equations (4) together with the interface condition (5) can be written in where δΓ t is the delta function of the interfaceΓ t and n is the unit normal vector to the interfaceΓ t (oriented into Ω 2 (t) ). The meaning of Eq. (6) is clear in the weak formulation: we seek u such that for any v ∈ H 1 (Ω) holds Furthermore, system (6) is equipped with an initial condition u = 0 and with suitable boundary conditions. Further, on the boundary ∂Ω the following boundary conditions are prescribed Level set method. The free surface motion is numerically treated with the aid of the level set method, see e.g. [11]. The level set function generally can be considered as the signed distance function being zero on the interfacê Γ t at any time instant t, positive in Ω 1 (t) and negative in Ω 2 (t) . In practical realization the level set function is used only to implicitly define the interfaceΓ t as the set of points satisfying the condition φ(x, t) = 0. With the aid of a level set function φ(x, t) the Heaviside function H(x, t) satisfy The motion of the interfaceΓ t is then realized by solving the transport equation for the function φ(x, t), i.e.
Equation (9) means that the interfaceΓ t moves in time with the flow velocity u.
Surface tension. The surface tension realized with the aid of the weak reformulation, see [4], [14]. To this end the tangent derivative ∇ Γ is defined as and similarly the second order tangent derivative, i.e. the Laplace-Beltrami operator, (see [14]) as Now, using the relation and applying the integration by parts onΓ t we get where on the right hand side the integral over ∂Γ t denotes the integral over the contact lines (in 3D). For the twodimensional problem the integral over ∂Γ t is the sum over the contact points, which -in the case ofΓ t being a closed curve -is the empty set and equals zero. In the other case this reformulation naturally includes the contact angles into the weak formulation of the problem.
The base functions of the extended finite element space shown for 1D case on the grid given by points xi−2, xi−1, xi and xi+1. The interface at time t is located between the nodes xi and xi+1 at x = c and the interval where the levelset function is positive (φ(x, t) > 0) is marked by the dashed part of the x-axis.

Numerical approximation
Flow step. In order to discritize the problem in time an equidistant division t n = nΔt of the time interval [0, T ) is considered with a constant time step Δt = T /N. Furthermore, the approximations of the flow velocity u, the presure p, the density ρ and the viscosity μ at the time instant t n are denoted by u (n) , p (n) , φ (n) , ρ n and μ n , respectively. The time derivatives are then approximated by the second order backward difference formula as The time discretized weak formulation of Eq. (6) then reads: Find u = u n+1 ∈ V and p = p n+1 ∈ Q such that holds for all v ∈ V and q ∈ Q.
In the practical computations the polygonal computational domain Ω is used and the spaces V and Q are approximated by the FE subspaces V h and Q h , respectively, defined over an admissible triangulation T h of the domain Ω. The Taylor-Hood finite elements are used for approximation of velocity and pressure, i.e. the velocity components are sought in the space and P k (K) denotes the space of all polynomials of degree less or equal to k on K. Further, the pressure (as well as the level set function) is approximated in the space The discrete flow problem then reads: Find u h = u n+1 h ∈ V h and p h = p n+1 h such that Eq. (11) holds for any test function v := v h ∈ V h and q := q h ∈ Q h . In order to treat the discontinuity of the pressure due to the presence of the surface tension the extended finite element method (XFEM) is applied, see e.g. [7].
Extended finite element method. The XFEM method is used to treat the possible numerical instability caused by the surface tension: the pressure is discontinuous along the interface due to the applied surface tension term. The numerical approximation with the aid of Taylor-Hood finite element means that the presure is continous. In order to capture the pressure discontinuity the pressure finite element space Q h is enriched using the localization of a enrichment function, in this case the enrichment of the Heaviside function H Γ (x) given as the Heaviside function H Γ (x) = H(x, t n+1 ) at the time instant t n+1 .
In order to enrich the finite element space, the original base functions of Q h are used, i.e. we denote the index set J = {1, . . . , n}, n = dim Q h and the mesh nodes by x j , j ∈ J . The nodal base functions are then denoted by q i ∈ Q h , i ∈ J and satisfy q i (x j ) = δ ij . The J is the subset of all the neighbours of the interfaceΓ t , i.e. J = {j ∈ J : supp q j ∩Γ t = ∅}.
Finally, the enrichment of the space Q h is made using the discontinuous base functions q xf em j defined by where the term −q j (x)H Γ (x j ) makes the support of the function q xf em j (x) localized to the single element containing the interface, see Fig. 2.
Here, H Γ (x j ) can be left out from the right hand side as this only adds a constant multiple of the continuous base function q j (x). On the other hand, this term makes the function q xf em j (x) being zero at every node x i , i ∈ J and also makes the support of q xf em j (x) localized only to the elements containing the interfaceΓ t , which simplifies the practical discretization of the problem. The FE space Q h is then replaced by the extended FE space   Level set step and coupled problem. The level set equation (9) is numerically approximated with the aid of finite element method. In order to get the numerically stability of the scheme, the algebraic flux corrections are applied, see [15]. Moreover, in order to keep the level set function as the signed distance function, a re-initialization step is applied, see also [4]. The solution of the coupled problem is then performed by a de-coupled algorithm with strong coupling. This means, that both problems are solved several times in order to get a convergence criteria satisfied. In practical computations the level set step is solved first followed by solution of equation (11) for approximation of the flow velocity u n+1 and the pressure p n+1 .

Numerical results
In order to verify the described numerical method, the two benchmarks for the case of a rising bubble considered in [4] were numerically approximated. The first test case (case I) used the following values of the fluid's densities Due to the gravitational force , the second fluid with the lower density starts to rise and deform its shape, see Fig. 3 However, due to high surface tension a stable shape is developed after a short time period. Since this time the fluid keeps rising in an undeformed shape, see Fig. 5.
Following [4] also the second benchmark test was approximated (test case II). In this case the fluid's densities were ρ 1 = 1000 kg m −3 and ρ 2 = 1 kg m −3 , fluid's viscosities were μ 1 = 10 Pa s and μ 2 = 0.1 Pa s. The same gravitational acceleration was used f = (0, −0.98) m s −2 but lower surface tenstion γ = 1.96 N/m was considered. The development of the free surface is shown in Fig. 4. The development of the shape as well as the vertical component of the flow velocity is shown in Fig. 6.
The computations were performed on triangular meshes with the equidistant partition. The spatial steps h = 1/20, 1/40 or 1/80 were used. The step h = 1/40 was the coarsest mesh used in [4]. The time step used in the computation was Δt = 0.002. The motion of the domain Ω 2 (t) with the area A(t) was tracked in terms of the y−coordinate of the center of mass T y (t) = In order to verify the presented numerical method the values of T y , C and V were computed at every time instant for both the test cases. The graphs of T y , C and V in dependence on time shown in Fig. 7 for the test case I agrees well with the results in [4]. The computations were performed on three grids and the differences were not significant. Similarly, the graphs of T y , C and V in dependence on time shown in Fig. 8 for the test case II is similar to the results presented in [4]. Nevertheless for the benchmark II the results are much more sensitive (particularly for the circularity) to the mesh refinement. See also the quantitative comparison for both cases in Tables  1, 2.

Conclusion
The motion of two immiscible fluid was mathematically modelled with the aid of a variational formulation which includes the surface tension and possibly also the contact angles effects. The time discretization of the problem was performed with the aid of the second order backward difference formula. The spatial discretization was done with the aid of the finite element method, where the XFEM was applied in order to capture the disconti-  nuity of the pressure along the free surface. The motion of the free surface was treated with the aid of the level set method. The presented numerical method was tested on two benchmark problems from [4] and the obtained numerical results were compared to the reference values.