A balanced-force control volume finite element method for interfacial flows with surface tension using adaptive anisotropic unstructured meshes

A balanced-force control volume finite element method is presented for threedimensional interfacial flows with surface tension on adaptive anisotropic unstructured meshes. A new balanced-force algorithm for the continuum surface tension model on unstructured meshes is proposed within an interface capturing framework based on the volume of fluid method, which ensures that the surface tension force and the resulting pressure gradient are exactly balanced. Two approaches are developed for accurate curvature approximation based on the volume fraction on unstructured meshes. The numerical framework also features an anisotropic adaptive mesh algorithm, which can modify unstructured meshes to better represent the underlying physics of interfacial problems and reduce computational effort without sacrificing accuracy. The numerical framework is validated with several benchmark problems for interface advection, surface tension test for equilibrium droplet, and ∗Corresponding author. Email: z.xie@imperial.ac.uk; Tel.: +44(0)20-75943067. Preprint submitted to Computers and Fluids July 27, 2016 dynamic fluid flow problems (fluid films, bubbles and droplets) in two and three dimensions.

transport equation of volume fraction with a differencing scheme without re-adaptivity can be found in Piggott et al. [19].
In this paper, we focus on the surface tension force model in three dimensions. Many different types of surface tension force model have been For each fluid component i, the conservation of mass may be defined as, and the equations of motion of an incompressible fluid may be written as: where t is the time, u is velocity vector, p is the pressure, the bulk density tions of the degrees of freedom for the P 1 DG-P 1 and P 1 DG-P 2 elements and 105 the boundaries of the control volumes in two dimension (2D).  In order to discretise the above governing equations, a finite element representation for u and p is assumed, expressed in terms of their finite element basis functions Q j and P j , respectively, as: Here N u and N p are the total number of degrees of freedom for the velocity 128 and pressure representations. Q j = Q j I, where I is the identical matrix.

129
When multiplying the global continuity equation with a quadratic con-130 tinuous Galerkin (CG) basis function P i and applying integration by parts 131 once, the discrete form of the global continuity equation can be obtained as: where n is the current time level, n is the outward-pointing unit normal can be obtained as: where n is the outward-pointing unit normal vector to the surface Γ E of the wher u in is the upwind velocity calculated from the neighbouring element or boundary and the subscript ( * ) represents the latest value during the iteration 145 in one time step. The advection term for the time level n can be obtained in 146 a similar way.

147
For the viscous terms τ = µ(∇u + ∇ T u), we use a high order linear 148 scheme which results in a compact stencil with an element coupling only to 149 its surrounding elements. For example, where τ nb is the value of τ in the neighbouring element along the face and 151 the viscous term for the time level n can be obtained in a similar way.

152
In order to evalute the viscous stress tensor τ on the boundary of element Γ E in Eq. (8), we integrate over the volume of two neighbouring elements in order to calculate the derivatives on the element face between the two elements. For example, the derivative in x for the x component of u is obtained as: in which Γ E1∩E2 is the the shared face between element 1 (E1) and element is not only validated for discontinuous elements but also can be applied for 156 continuous elements.

Projection method 158
The discretised form of the momentum (Eq. (6)) and global mass balance 159 (Eq. (5)) equations are solved using a pressure projection method. This effec-160 tively eliminates the unknown velocity and solves a system of equations for 161 pressure or pressure correction. The discretised momemtum and continuity 162 equations, at time level n + 1, can be written in matrix form respectively, as: where u n+1 and p n+1 are the FEM solution fields for velocity and pressure, The matrix equation for velocity to be satisfied is: Subtracting these two equations the velocity correction equation is obtained:  It is worth mentioning that in the first time step after the mesh adaption, 181 the interpolated velocity field on the new mesh is projected to a continuity-182 satisfied space, otherwise it will not satisfy the divergence free condition.
where σ is the surface tension coefficient, κ is the interfacial curvature,ñ is the interface unit normal, δ is the Dirac delta function. Here, we use δ = |∇α| andñ = ∇α |∇α| to reformulate the CSF based on the component volume fraction as: In this approach, the curvature is estimated as: whereñ df = ∇ϕ |∇ϕ| is the interface unit normal calculated from the signed 222 distance function ϕ from the interface (similar to the level set function).

223
Contrary to the standard CLSVOF method, only one function (volume frac-224 tion α) is advected to capture the interface here, and then a distance function 225 ϕ is calculated based on the volume fraction α by the following three steps.

226
In the first step, the distance function is initialised as: where h min is the minimum mesh size around the interface. As we assume 227 the contour of α = 0.5 is the interface, it can be seen that Eq. (20) provides 228 a good estimate for the initial distance function, where it is zero at the 229 interface (required by the definition) and has a different sign on either side of 230 the interface with a good guess for the distance in the vicinity of the interface.

231
In the second step, we follow the re-initialisation process in the level set method to obtain the actual sign distance function ϕ by solving the following equation with the initial value ϕ 0 : where S(ϕ 0 ) = ϕ 0 / ϕ 2 0 + ε is the sign function and ε = h min is used here In addition, we use the pressure basis functions as it is high order.

238
In the last step, the calculated discontinuous distance function (ϕ DG ) is projected to a continuous space (ϕ CG ) by a volume-weighted interpolation as: where N cv is the number of control volumes connected to the continuous 239 function ϕ CG and V i is the volume of the control volume (see Fig. 1 for 240 example). Then the continuous function ϕ CG is used to calculate the gradient, 241 which helps to accurately estimate the unit normalñ df and the curvature κ.

b. Diffused interface approach
In this approach, we use a diffused interface based on the original volume fraction α to estimate the curvature as: whereñ di = ∇ψ |∇ψ| is the interface unit normal calculated from the diffused volume fraction ψ from the interface. In contrast to some smoothing techniques for the volume fraction to calculate the curvature, we obtain the diffused volume fraction by solving a diffusion equation as: where τ is also an artificial time and D is an artificial diffusion coefficient.
Here, we initialise the diffusion value based on the original volume fraction α as:  Finally, similar to the distance function approach, the calculated discon-251 tinuous diffused volume fraction is projected to a continuous space, which is 252 used to approximate the unit normalñ di and the curvature κ.
where u and v are the horizontal and vertical components of the velocity 281 field, respectively.

282
The initial interface shape is deformed by the velocity field and return to very elongated filament which is very thin at the tail. The interface has been 291 efficiently captured in the computation by using the adaptive unstructured mesh, which provides fine resolution equivalent to that of a 512×512 uniform 293 mesh. After t = 2, the velocity field is reversed and the interface shape is 294 returned to its initial shape which is well captured during the computation.

295
This test demonstrates the power of the adaptive mesh approach, which can 296 refine the mesh in the vicinity of the interface or an area of interest, and 297 reduce computational effort without sacrificing accuracy.

298
As the location of the interface is known (which is represented as the contour of the volume fraction at α = 0.5), the deviation of the interface position after one rotation can be calculated as: where N is the number of points along the interface. Fig. 3 shows the conver-299 gence for the computations with three different simulations. It can be seen 300 that the present method is close to second order accurate, which is consistent 301 with the quadratic polynomial function used for the P 2 finite element type.

Static drop in equilibrium 303
In order to validate the proposed framework for surface tension, we con-  putational mesh with 40 layers in each direction is used for the computations with the P 1 DG-P 2 element pair.

312
In the first set of calculations, we test the coupling of surface tension force with the pressure gradient by specifying the exact curvature κ exact = 1/R = 0.5. Table 1 shows the comparison for the maximum velocity and pressure jump errors after one time step for three different density ratios, and also with the results obtained for the structured Cartesian mesh in Francois et al.
[22]. These errors are defined as: where ∆P exact = σκ = 36.5, and w denotes different evaluation ways by 313 using pressure points in the areas of r < R and r > R (total) and the areas 314 of r < R/2 and r > 3R/2 (partial), where 'total' considers the whole region 315 inside and outside the drop and 'partial' considers some parts of the region 316 by avoiding the transition zone. It can be seen from Table 1 that the spurious   317 currents are very small, related to the machine accuracy. For the pressure 318 drop, the total error E(∆P ) total is of the order of 10 −2 , and is independent 319 of the density ratio. The partial error E(∆P ) partial is much smaller because 320 the error measurement does not include the transition region.  same setup above for density ratio ρ 1 /ρ 2 = 10 3 . Table 2 shows the results for 327 four different techniques: the first is for the given exact curvature; the second 328 is for the curvature calculated from the analytically given distance function 329 (the distance function is given rather than solving from Eq. (21)); the third 330 and fourth are for the distance function approach and diffused interface ap-331 proach for the curvature calculation proposed in this study, respectively. It 332 can be seen from Table 2 that the two approaches for the curvature cal-333 culation are in the same order of accuracy, whereas the diffused interface 334 approach performs slightly better than the distance function approach, and 335 even better or close to the case when the distance function is specified. The  Table 2: Error in velocity and pressure drop after one time step with ∆t = 10 −6 for the inviscid static drop in equilibrium for density ratio ρ 1 /ρ 2 = 10 3 with different curvature calculation techniques for the structured mesh.

347
As shown in Fig. 4 and     no-slip boundary condition is applied at both walls in the water and air side, and a zero-gradient boundary condition is applied at the outlet. At the inlet x = 0, a parabolic velocity profile is imposed for the liquid phase as: where u av is the average velocity containing the forcing perturbation: in which u N is the Nusselt velocity, ξ is the disturbance magnitude and F is the forcing frequency. The velocity for the air phase at the inlet is set as: The simulation is initialised with a flat film with a fully developed velocity 353 field as prescribed at the inlet. The computation is carried out using the     Buckmaster [34] (case A in Table 1 in [34]), which have been used to validate

Coalescence of two bubbles
Here we consider merging of two bubbles with co-axial and oblique coales-  For the oblique coalescence case, Fig. 10      The results presented here established with sufficient confidence that this 483 method can be used to successfully model multiphase flows in a wide range 484 of applications. This approach has the potential to be used for an arbi-485 trary number of components although that has not been demonstrated here.

486
Future work will include parallel computing and surfactant modelling. Nuclear Energy project (EP/M012794/1) and H2020 IVMR. Funding for Dr 492 P. Salinas from ExxonMobil is gratefully acknowledged.