MONOLITHIC NEWTON-MULTIGRID SOLVER FOR MULTIPHASE FLOW PROBLEMS WITH SURFACE TENSION

: We have developed a monolithic Newton-multigrid solver for multiphase ﬂow problems which solves velocity, pressure and interface position simultaneously. The main idea of our work is based on the formulations discussed in [14], where it points out the feasibility of a fully implicit monolithic solver for multiphase ﬂow problems via two formulations, a curvature free level set approach and a curvature free cut-oﬀ material function approach. Both formulations are fully implicit and have the advantages of requiring less regularity, since neither normals nor curvature are explicitly calculated, and no capillary time restriction has to be respected. Fur-thermore, standard Navier-Stokes solvers might be used, which do not have to take into account inhomogeneous force terms. The reinitialization issue is integrated within the formulations. The nonlinearity is treated with a Newton-type solver with divided diﬀerence evaluation of the Jacobian matrices. The resulting linearized system inside of the outer Newton solver is a typi-cal saddle point problem which is solved using a geometrical multigrid method with Vanka-like smoother using higher order stable Q 2 /P disc 1 FEM for velocity and pressure and Q 2 for all other variables. The method is implemented into an existing software package for the numerical simulation of multiphase ﬂows (FeatFlow). The robustness and accuracy of this solver is tested for two diﬀerent test cases, static bubble and oscillating bubble, respectively.


INTRODUCTION
Multiphase flows are of great interest in different industrial and engineering applications. The simplest example of multiphase flow is two-phase flow [6], two fluids/phases are separated by an interface, where surface tension forces are applied. If the fluids have different densities and viscosities, then a discontinuous pressure jump is observed near the interface. The interface moves/deforms due to the flow movement, and to capture this behaviour, an efficient tracking/capturing method should be applied. There are two main methods for interface modeling in the multiphase flow problems, i.e. Lagrangian and Eulerian methods. The volume of fluid (VOF) [8], phase field [1,2] and level set [13,15] are among the most famous Eulerian interface capturing methods, which are very favourable for the computational and implementation point of view.
In the present work, numerical methods based on the level set and material cut-off function for two-dimensional incompressible two phase flow are implemented. The approximation of the surface tension force does not require the calculation of the curvature, normals and the delta function. However, the level set method requires some sort of redistancing [17]. In our improved fully implicit level set method, the explicit redistancing is removed by integrating the reinitialization term within the formulations [14]. Moreover, the advantage of this approach is that there is no capillary time restriction and the standard Navier-Stokes solver can simulate Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain the multiphase problems with homogeneous force terms. The numerical studies are carried out for two different test cases, static bubble and oscillating bubble, which show the accuracy and robustness of these formulations in the context of FEM. The system of equations in each formulation is solved monolithically (in a fully coupled manner).

GOVERNING EQUATIONS
In the methodology of our work, first the Continuum Surface Force (CSF) [5] is introduced and then it is linked to the classical Continuum Surface Stress (CSS) [10]. In the CSF approach, the interface between the fluids is not considered as a sharp discontinuity but as a smooth transition. As a result, the surface tension is also assumed to be continuous everywhere in the transition regime. A detailed discussion of this approach can be found in [4]. In the CSS approach, a stress tensor is introduced and the surface force term is written as the divergence of the stress tensor [4]. The following formulations in section (2.1) and (2.2) are based on the CSS approach.

Curvature free level set approach
The curvature free level set formulation [14] is introduced by adding a tensor field in the Navier-Stokes equations. The system of equations is defined as: Here, ρ is the density, u is the velocity, τ = (τ s + τ m ) is the full stress tensor, p is the pressure, φ is the level set function, ψ is the cut-off function and ψ is the parameter for the interface thickness. The standard stress tensor τ s and the modified stress τ m (derived in [14]) are defined as where µ is the viscosity, D(u) = 1 2 (∇u + ∇u T ) is the deformation stress tensor and σ is the surface tension coefficient. The modified pressure p m balances the pressure peaks at the interface. The mathematical expression is defined as In order to circumvent the explicit reinitialization, the additional normal diffusion term can be integrated into the level set equation (1). The complete system of equations is defined as follows: Here, γ nd is the relaxation parameter for normal diffusion and this term has the forward and backward diffusion property [11]. The system of equations (4) is a four field system with the unknowns (u, φ, ψ, p) T . The two main advantages of this formulation are that neither normals nor curvature have to be explicitly computed, which are the sources of numerical errors.

Curvature free cut-off material function approach
In this approach, we are no longer in the need of the level set function, instead we use a new equation for the cut-off material function. Olsson and Kreiss [12] have introduced the equation for the material cut-off function with fictitious time as where γ nc and γ nd are the relaxation parameters for the nonlinear convection in the normal direction and the normal diffusion, respectively. The nonlinear convection in the direction of the normal has the tendency to build the Heaviside step function, without depending on the convective parameter. Whereas, the sharpness of the interface is controlled by the normal diffusion. The full set of equations including the material cut-off function in physical time is defined as follows: The system of equations (6) is a three field system with the unknowns (u, ψ, p) T . The momentum equation has homogeneous force terms. We are solving the systems in a fully coupled manner with our monolithic multiphase flow solver.

NUMERICAL METHOD
For solving the system of equations (4) and (6), first we discretize in time with a fully implicit 2nd order time stepping scheme, i.e. Crank Nicolson. For the space discretization, the velocity and pressure fields are discretized using higher order stable Q 2 /P disc 1 FEM [3] and Q 2 for the level set function as well as for the cut-off material function, presented in Fig. 1. The nonlinearity in the system of equations (4) and (6) is treated by Newton solver and the resulting linear system is then solved using multigrid solver. Our Newton-multigrid solver is fully monolithic, which means it solves all the variables (u, φ, ψ, p) simultaneously.

NUMERICAL RESULTS
Two numerical studies are performed to assess the accuracy and robustness of the solver. The static and the oscillating bubble are the prototypical configurations and present the behaviour of two phase flows.

Static bubble
A static bubble in a two-dimensional incompressible two phase flow is considered [18,7,16,9]. This is a simple example for demonstrating the prototypical behaviour of multiphase flows. For simplicity, we consider a stationary bubble at equilibrium. Since the bubble is at rest inside the domain so it should show a zero velocity field but unfortunately spurious velocity/currents are observed near the interface [18,7,16,9]. Moreover, flows involving interfaces lead to large pressure jumps. Approximation of the incompressibility constraints, the interface and the local external force are three different responsible sources for these phenomena [7]. The flow dynamics depends strongly on the magnitude of the spurious velocities. A non-physical movement of the interface might also be observed due to this spurious velocity.

Geometrical configuration
Both fluids are immiscible and separated by an interface Γ. The first fluid Ω 1 is completely inside of the second fluid Ω 2 as shown in Fig. 2 (lef t). A circular static bubble of radius r = 0.25 is placed at the center [0.5, 0.5] of a unit square Ω = [0, 1] 2 . The surface tension coefficient, viscosities and densities are set to unity in the absence of the gravitational force. The relation between the pressure (inside and outside of the static bubble) should satisfies Laplace Young law: Here, p i is the pressure inside the bubble, p o is the pressure outside the bubble and σ is the surface tension coefficient. The numerical studies for the systems of equations (4) and (6) are performed with a fixed time step ∆t = 10 −2 until t = 10. The interface thickness is controlled by the parameter ψ in both formulations. The spurious velocity/currents are observed and visually represented in Fig. 3 (a,b) and Fig. 4 (a,b). It can be seen in Fig. 3 (c,d) and Fig. 4 (c,d) that the pressure difference inside and outside of the bubble converges to the magnitude 4. Hence, the pressure jump across the interface successfully satisfies the Laplace Young law. By decreasing the value of interface thickness parameter ψ , the surface tension force becomes sharp, resulting in a sharp pressure jump. The graphical representation of the pressure can be seen by a cross section through y = 0 in Fig. 3 (c,d) and Fig. 4 (c,d). Moreover, the cut-off function exhibits the same behaviour w.r.t. the interface thickness parameter ψ in Fig. 3 (e,f) and 4 (e,f).

Oscillating bubble
A two-dimensional unsteady incompressible two phase flow is considered. These fluids are immiscible and separated by an interface Γ. The first fluid Ω 1 is completely inside the second fluid Ω 2 as shown in Fig. 2 (right).

Geometrical configuration
The circle is initially perturbed to an elliptical shape with the semi-axis 0.25 in the xdirection and 0.125 in the y-direction. The ellipse is placed at the center [0.5, 0.5] of a unit square Ω = [0, 1] 2 . The surface tension coefficient, viscosities and densities are set to unity. The gravitational force is neglected.
The numerical studies for the systems of equations (4) and (6) are performed with a fixed time step ∆t = 10 −2 until t = 100. The transition of the radius (r x , r y ) to a steady state with respect to time, is presented in Fig. 5. This study is performed for three different mesh refinement levels, i.e. h = 1/16, 1/32 and 1/64. As the mesh is refined, the smooth transition from oscillating to steady state is illustrated in Fig. 5. It is observed that after certain oscillations the bubble reached the steady state, with no mass loss. It can be observed that the numerical instability arising from the interface capturing in the level set method vanishes in formulation (6), so the oscillating bubble reaches the steady state with much less oscillations. To analyse the temporal development of the interface, the bubble shapes are extracted at different time intervals, for mesh refinement level 6 (h = 1/64) in Fig. 6. At the final time, the ellipse is expected to reach an equilibrium state, that is a stable circular shape. As expected, the oscillating bubble is transformed into a stable circular shape.  : Results for the system of equations (4(lef t), 6(right)): The radius (r x , r y ) of oscillating bubble using Crank Nicolson time stepping scheme for three different mesh refinement levels, ∆t = 10 −2 . Figure 6: Results for the system of equations (4(lef t), 6(right)): Oscillating bubble shape with respect to time t, ∆t = 10 −2 .

CONCLUSIONS
In this work, we have developed a monolithic Newton-multigrid solver for multiphase flow problems which solves velocity, pressure and interface position simultaneously. There are three main advantages of these formulations. Firstly, no explicit computation of the curvature and normals are required. Secondly, the explicit redistancing is removed by integrating the reinitialization term within the formulations. Thirdly, there is no capillary time restriction. In order to investigate the accuracy and robustness of the solver, numerical studies are performed for two-dimensional static and oscillating bubble. The results expectedly confirms the accuracy of the solution approximation as the magnitude of the pressure across the interface satisfy the Laplace Young's law. In both test cases, the bubble reached its steady state satisfying the theoretical predictions.