A Low Mach Number IMEX Flux Splitting for the Level Set Ghost Fluid Method

Considering droplet phenomena at low Mach numbers, large differences in the magnitude of the occurring characteristic waves are presented. As acoustic phenomena often play a minor role in such applications, classical explicit schemes which resolve these waves suffer from a very restrictive timestep restriction. In this work, a novel scheme based on a specific level set ghost fluid method and an implicit-explicit (IMEX) flux splitting is proposed to overcome this timestep restriction. A fully implicit narrow band around the sharp phase interface is combined with a splitting of the convective and acoustic phenomena away from the interface. In this part of the domain, the IMEX Runge-Kutta time discretization and the high order discontinuous Galerkin spectral element method are applied to achieve high accuracies in the bulk phases. It is shown that for low Mach numbers a significant gain in computational time can be achieved compared to a fully explicit method. Applications to typical droplet dynamic phenomena validate the proposed method and illustrate its capabilities.


Introduction
In many situations in nature or technical applications, more than one material phase is present. Often, two phases are separated by a distinct interface as it is the case for, e.g., rain drops or spray related processes. In such applications, the characteristic Mach number Ma = ‖u‖ 2 c , relating the velocity of the convective phenomena u to the speed of sound c, often covers a wide range of values. Especially in the liquid phase, the Mach numbers are often low, whereas in the gaseous phase higher Mach numbers are present. The simulation of such phenomena requires numerical methods which are capable of treating flows with a wide range of Mach numbers accurately (including the correct asymptotic scaling) and efficiency. Considering the characteristics of the Euler equations shows the difficulties numerical methods are facing. For large Mach numbers ( Ma = O(1) ) the Euler equations are hyperbolic. In the limit ( Ma → 0 ) they change their type to hyperbolic-elliptic. For small Mach numbers ( Ma ≪ 1 ), the speeds of the present characteristic waves differ tremendously: the acoustic waves are much faster than the convective phenomena. If acoustics plays a minor role, standard explicit methods are not efficient as they have to resolve the fast characteristic waves for stability reasons, leading to a prohibitively small timestep.
To overcome this timestep restriction, different approaches exist, mostly relying on the implicit time discretization: fully implicit time discretization for the modeling of two-phase flow are, e.g., proposed in [25,31,47]. Another class of schemes are pressure-based methods, relying on the implicit solution of an elliptic equation for the pressure. They can be seen as an extension of incompressible schemes and require a reformulation of the equation of state (EOS), which is a non-trivial task for complex EOS. Pressure-based methods have been derived for multiphase flows, e.g., in [4,15,18,26,35,45]. Another option are implicit-explicit (IMEX) flux splitting schemes which are based on the separation of the fast acoustic and slow convective waves of the flux. The idea is to treat the fast waves implicitly, while the slow waves are treated explicitly. This splitting often serves as a starting point for the derivation of pressure-based methods, see, e.g., [13,30,51,64]. The application of flux splittings with general EOS is often straightforward as the EOS does not have to be reformulated. Nevertheless, IMEX schemes have only rarely been applied to multiphase problems, see [7,8,55].
Low Mach number two-phase flows not only necessitate a specific design to tackle the different orders of magnitude of the characteristic waves. Additionally, the presence of two phases demands a modeling of the different states of matter and their coupling. For the modeling of multiphase flows mainly two concepts exist: diffuse and sharp interface approaches. Common approaches to capture the interface in a sharp interface scheme are the volume-of-fluid [33] and the level set method [54,62]. The level set method has the advantage that the normals and the curvature of the phase interface are directly given by the derivatives of the level set function and hence the interface geometry does not have to be reconstructed, as it would be the case for the volume-of-fluid method. For the modeling of the physical phenomena at the phase boundary, the ghost fluid method (GFM) has been introduced in [23]. In this work, we focus on a specific type of a sharp interface GFM based on the ideas of [49]. Here, a multiphase Riemann problem is solved at the phase boundaries. The idea has been adapted in [20,21] and slightly modified in [37,50]. A discontinuous Galerkin scheme [32] is used for the spatial discretization of the bulk phases and at the phase boundary a finite volume sub-cell refinement [60] is applied. This framework will serve as a foundation of this work.
Since IMEX flux splitting schemes have been found to be well suited for low Mach number computations within a discontinuous Galerkin framework [72,73], the idea in this work is to combine an IMEX flux splitting scheme with the sharp interface level set ghost fluid method (LSGFM). Although the fully implicit time discretization has been used within the LSGFM [31,37,47], IMEX flux splitting schemes have not been investigated in this context. As a flux splitting, the idea of Toro and Vázquez-Cendón (TV) [67] and its extension to general EOS [66] are used. Hence, the main novelty of this paper is the analysis and presentation of how an IMEX flux splitting can be combined with a discontinuous Galerkin sharp interface LSGFM to overcome the restrictive timestep restriction of weakly compressible multiphase flows.
We will proceed as follows. In Sect. 2, the considered governing equations are summarized. Following in Sect. 3, the IMEX LSGFM is described. Starting with a short summary on the level set ghost fluid framework in Sect. 3.1, it is described how IMEX flux splitting schemes can be applied in Sect. 3.2. Next, novel scheme is validated and its efficiency is investigated in Sect. 4. Following, some illustrative applications are shown in Sect. 5. Finally, in Sect. 6 conclusion and outlook are given.

Governing Equations
The idea of the sharp interface approach is to model two-phase flows as two pure bulk phases without a mixing zone. Inside the pure phases, the Euler equations are chosen to model the physical behavior of the fluids. In three dimensions they are given by with density , velocity u = (u, v, w) T , total energy E, pressure p and the three dimensional identity matrix . The primitive state vector is given as w prim = ( , u, p) T . The pressure is linked to the density and the specific internal energy ∶= E − ‖u‖ 2 2 2 via an EOS. For gases, the ideal gas EOS with the isentropic coefficient is chosen. Liquids are modeled with the stiffened gas EOS: where the material-dependent parameter p ∞ is used to model the increased stiffness of liquids compared to gases. Note that the algorithm introduced in this work is not restricted to the use of such simple EOS. It is also possible to use more complex or tabulated EOS [24].
The phase interface separating the two phases is tracked by a level set function Φ , describing the signed distance to the phase interface. Hence its root describes the position of the phase interface. The level set function is transported with a velocity field according to

3
The geometrical quantities of the phase interface, i.e., the normal vector Φ and the curvature Φ , are directly given by the derivatives of the level set field [17].

IMEX Flux Splitting for the LSGFM
In this section, first the compressible LSGFM according to [21,37,50] is briefly recalled. Following, it is outlined how IMEX flux splitting schemes can be used to enable efficient simulations at low Mach numbers.

Domain Decomposition
The position of the phase interface is described by the zero position of the level set function. Hence, the sign of the level set field can be used as an identifier to divide the computational domain into a liquid l and a gaseous domain g . Discretely, the interface position is shifted to the next boundary of the discretization elements, in the following called the surrogate surface. To obtain an accurate scheme, the high order discontinuous Galerkin spectral element method (DGSEM) [32,41] is applied for the spatial discretization. The DGSEM is based on the polynomial approximation of the solution as a tensor-product of nodal one-dimensional Lagrange polynomials of degree N inside each element C of the computational mesh C ∈ . For this piecewise smooth function w h the weak formulation is then given by for every polynomial test function (x) of degree N. Across the element boundaries discontinuities are allowed and the elements are coupled via a numerical flux function f * . Choosing the same nodal Lagrange polynomials of degree N as test functions and ansatz functions, leads to the DGSEM which is (N + 1)th-order accurate in space. A detailed description of the used DGSEM can be found in [32] and an overview on the actual implementation is provided in [42]. The core method is available as open source software FLEXI 1 . Due to its solution representation, each grid element contains multiple (nodal) degrees of freedom (DOFs). In the LSGFM context this leads to a poor relation of discrete interface resolution per DOF. Therefore, a refinement of the elements adjacent to the phase interface is applied. For that purpose, the mixed finite volume (FV)/DG formulation from [60,61] is used. Figure 1 illustrates this sub-cell refinement in the cells neighboring the phase interface and the domain decomposition in a liquid and a gaseous domain l and g , respectively which are separated by the surrogate surface .

Interface Treatment
Two-phase Riemann solver While in the pure bulk phases, no multiphase effects have to be considered, directly at the discrete phase interface a special treatment of the fluxes is necessary. The procedure at the phase boundary is based on the idea outlined in [49], to use the solution of a two-phase Riemann problem as a boundary condition for the pure phases. This approach has been modified in [20,21] to use the approximate solution of the two-phase Riemann problem instead of the exact solution. In Fig. 1, positions on the surrogate surface, where two-phase Riemann problems are solved are marked with a small white square. In this work, the HLLP Riemann solver from [20,59] is used at the phase boundary. This Riemann solver is based on the well-known HLLC Riemann solver, where the contact discontinuity "C" is considered as the phase boundary "P". It enables the treatment of surface tension effects by considering the pressure jump of the Young-Laplace law in the jump conditions of the Riemann solver. From the solution of the HLLP Riemann solver three quantities can be extracted: fluxes for the liquid and the gaseous phase f * l , f * g , the velocity of the phase interface s # and so-called ghost states for the liquid and the gaseous phase. For a liquid cell adjacent to the interface, its ghost state is a gaseous state and vice versa. Here, the ghost states are chosen as the interchanged initial conditions of the two-phase Riemann problem. They are used if the interface position moves and a sub-cell changes its assignment to the liquid or gaseous domain.
Note that the GFM is inherently non-conservative [23,49]. This stems from the nonequality of the fluxes at the interface f * g ≠ f * l and from the definition of the new states if a sub-cell element changes its assignment from the liquid to the gaseous domain or vice versa. Attempts to overcome the non-conservativity of the GFM for incompressible flows relying the divergence-free condition of the incompressible velocity field have been made in, e.g., [16,52,53]. For compressible flows, an attempt to improve the conservativity has been presented in [46]. The idea in [46] is to modify the novel states which are set when the interface moves. Nevertheless, it is not considered here as this modification can have an unfavorable influence on the stability. Ω Ω g l Velocity extrapolation From the solution of the two-phase Riemann problem, the velocity at the phase boundary s # can be extracted. This pointwise available information is then used to construct a velocity field for the transport of the level set field. This is done by solving a Hamilton-Jacobi-type equation [56] with the numerical Hamiltonian according to [19], with the componentwise absolute value of the level set normal | Φ | and a numerical sign function sgn . The discretization of the up-and downwind gradients ∇ + x and ∇ − x is done with the HJ-WENO scheme from [36] and the pseudo time iteration is performed with a 3rd order explicit Runge-Kutta method. More details on this can be found in [19].
Level set transport: The resulting velocity field is then used for the solution of Eq. (2). Differently to the weak form of Eq. (1) given in Eq. (3), the variational form for the level set transport reads As this is a non-conservative equation, the discretization has to be slightly adapted. Similarly to the mixed DGSEM/FV formulation used for the bulk phases, the solution Φ h is either represented with the tensor-product of nodal one-dimensional Lagrange polynomials of degree N (DGSEM) or with the piecewise constant refined FV representation. When choosing a path-conservative formulation for the discretization of the weak form introduces now modified numerical fluxes at the cell boundaries to account for the non-conservativity of the equation. Details on the derivation of this specific fluxes can be found in [37].
Level set reinitialization Unfortunately, the level set field does not retain its signed distance property during the advection and hence it has to be reinitialized. Again, a Hamilton-Jacobi equation is solved with the WENO scheme from [36] and an explicit 3rd order Runge-Kutta scheme. The Hamiltonian H reinit (∇ x Φ) = sign(Φ) |∇ x Φ| − 1 is approximated with the numerical Hamiltonian according to [17] by

3
The reinitialization and the velocity extrapolation procedure are proceeded in pseudo time until information from the level set zero position has propagated two to three mesh elements away from the phase interface.
Interface derivatives calculation As stated before, the geometrical quantities of the phase interface, i.e., normal vector and interface curvature, are directly given by the derivatives of the level set field. It has been found to be beneficial to use second order finite differences for the calculation of the derivatives, although the level set field is discretized with higher order polynomials [48]. This approach is also pursued here, i.e., the gradients of the level set, required for the normals of the phase interface Φ = ∇ x Φ∕|∇ x Φ| , are discretely calculated at a DOF ijk with where Φ x 1 denotes the derivative of the level set in x 1 -direction. The level set normal vector is used to transform the states for the solution of the two-phase Riemann solver such that the Riemann problem is solved in the phase normal direction, i.e., in the normal direction of the physical phase interface. This ensures that the pressure jump introduced by the surface tension acts in correct direction. The pressure jump is calculated via Young-Laplace law using the curvature of the phase interface where Φ x 1 x 2 denotes the derivative 2 Φ x 1 x 2 and is obtained by applying the finite difference operator in Eq. (4) twice.
These procedures now serve as building blocks for the construction of the IMEX LSGFM, which is introduced in the next subsection.

The IMEX LSGFM
The idea of IMEX flux splittings is to separate the fast (acoustic) and the slow (convective) waves. To obtain a timestep restriction independent of the fast waves, they are then treated implicitly and the slow waves are treated explicitly. Fast waves are presented in Eq. (1) but not in Eq. (2). Therefore, the level set transport is treated purely explicitly and only for the time discretization of the Euler equations and an IMEX flux splitting scheme is applied.

Overview of the Algorithm
In Fig. 2, it is illustrated how the different building blocks can now be assembled to a consistent algorithm, where the level set transport is handled explicitly. In a first step fully implicit time discretization is applied to the flow field. The separation of the fast and slow waves using a flux splitting and leading to an IMEX discretization of the flow field is described thereafter.
(i) Update the DG/FV distribution of the level set field according to the information about the position of the interface. (ii) Advance the level set field explicitly in time with the velocity field to obtain Φ n+1 by solving Eq. (2). (iii) Define l and g at the new time instance t n+1 according to the sign of the level set field. Use the ghost states from the old time instance t n if an assignment of a sub-cell has changed. (iv) Repeat step (i) for the flow solver. (v) Reinitialize the level set field to recover its signed distance property. (vi) Calculate the geometrical properties n+1 Φ and n+1 Φ of the phase interface at the new time instance. (vii) Iteratively advance the flow field using the implicit/IMEX time discretization.
-Use the geometrical information n+1 Φ and n+1 Φ for the solution of the two-phase Riemann problems to obtain f * l , f * g , s # and the ghost states at the new time instance t n+1 .
Advance bulk phases (viii) Construct the velocity field for the level set transport from s # , when the iterative procedure from step (vii) has converged. (ix) Repeat steps (i)-(viii) until the final time t end is reached.
This quasi-serial advancement of the level set function and the flow field ensures that the interface position of the initial guess for the implicit method is already at the correct location for t n+1 . Note that this coupling is first order in time. Nevertheless, the errors of the temporal discretization can be reduced by applying Runge-Kutta methods for the level set transport and the advancement of the flow field.

TV Flux Splitting
If a proper splitting of the Euler equations (1) is chosen and coupled with the above described algorithm, this allows timesteps being independent of the fast acoustic waves. There are a lot of options of how to split the flux of the Euler equations into an implicitly F (w) and explicitly treated flux F (w) with F(w) =F(w) +F(w) , see, e.g., [5,12,29,67,74]. In this work, the TV splitting approach according to [67] is used due to its simplicity and its ability to be applied to general EOS [66]. It splits the flux into a pressure system F (w) and an advection system F (w) via According to [66], the eigenvalues of the implicit subsystem are given by with For the explicit subsystem the eigenvalues can be calculated as From this one can see, that only convective phenomena influence the timestep if only the advection system is treated explicitly as the timestep of the DGSEM is given by for a characteristic length scale l ref of the discretization elements and CFL ⪅ 1 . For the DGSEM, the CFL number is typically scaled with an empirically obtained constant which is specific for the chosen Runge-Kutta method and polynomial degree, see, e.g., [28]. In this work we do not use optimized scaling factors for the IMEX Runge-Kutta methods and simply choose the scaling factor 1. Therefore, the stability limit is not strictly CFL = 1 but can be slightly below or above. Additionally, a numerical flux function for the subsystems is required. We choose the standard TV flux proposed in [66] or the less dissipative TVHLL-LM Riemann solver which are given in Appendix A.

Spatial Discretization
With the TV flux splitting only the flux in the bulk phases can be split into an explicit and an implicit part. For the flux obtained by the two-phase Riemann solver ( f * g , f * l ), this is not possible for general two-phase Riemann solvers as physical modeling of two-phase effects is considered. This can include surface tension such as, e.g., in [20,59] and phase transition, see, e.g., [22,34]. Therefore, one can either treat this flux fully implicitly or fully explicitly. The latter would cause a timestep restriction depending on the acoustic waves and hence this flux is treated fully implicitly. Moreover, sub-cell elements that change their assignment from the liquid to the gaseous domain or vice versa cannot be treated partially implicit and partially explicit. Therefore, the cells located around the phase interface also have to be treated fully implicitly. This is realized with a narrow band of two elements adjacent to the interface in each direction, where no flux splitting is applied. This is illustrated in Fig. 3. Away from the interface, the TV splitting is used, indicated with the hatched area in Fig. 3. Inside the fully implicit narrow band, a finite volume sub-cell discretization is applied, whereas the DGSEM is used in the remaining IMEX flux splitting domain.
For the coupling of the IMEX flux splitting and the fully implicit domain, a special treatment is necessary. In analogy to domain-based IMEX methods, see, e.g., [38], the flux at the interface is calculated separately for both domains. That means that, during the flux calculation the neighboring domain is assumed to belong to the same type of the domain as the current element. With this, the numerical fluxes can be calculated as With this approach conservativity at the coupling interface cannot be guaranteed. However, since the ghost fluid method is not conservative by construction, this influence is supposed to be small. In Sec. 4.3.2 this is experimentally justified (see Fig. 9 right). Note that choosing the TV flux instead of the TVHLL-LM flux in the sub-cell directly neighboring the phase interface, has a favorable influence on the stability. Therefore, in those sub-cells always the TV flux is chosen as a Riemann solver. In the remaining domain one can either use the TV or the TVHLL-LM flux.

Temporal Discretization
In Fig. 2, the first order in time procedure has been shown. To increase the temporal accuracy, IMEX Runge-Kutta schemes [39] can be used for the advancement of the bulk phases and the level set field. Note that the coupling of both fields is still done only once per timestep and hence the global order of accuracy is not increased. Nevertheless, one can expect more accurate results by using higher order IMEX Runge-Kutta schemes.
In this work, globally stiffly accurate IMEX Runge-Kutta methods from [39] are used. The time integration procedure in semi-discrete form for a scheme with s stages, integrating from t n to t n+1 , can be written as (i) for i = 1, ⋯ , s , solve for w n,i (ii) set w n+1 ∶= w n,s .
The coefficients Ã ij , Â ij , c j and ĉ j are given in the Butcher table of the chosen method. Note that the second step only holds due to the globally stiffly accurate property of the chosen IMEX Runge-Kutta methods. The arising non-linear equation system with b n,i denoting the sum of all values known from previous stages, is solved with a Jacobian-free Newton-GMRES method, see, e.g., [40] for an overview. Newton's method for the solution of Eq. (8) requires an initial guess w 0 and a user defined tolerance Newton and is given by dg(w k ) dw k Δw k = −g(w k ), w k+1 = w k + Δw k ; As initial guess, the solution at the previous stage is chosen, i.e., w 0 ∶= w n,i−1 . The linear system in Eq. (9) is then solved with a Jacobian-free GMRES method to ensure a small memory footprint. An analysis on the memory requirements of a matrix-free and a matrix-based approach for the DGSEM can be found in [68]. The arising Jacobian-vector product in the GMRES method is approximated via a first or a second order finite difference thus requiring one or two evaluations of the spatial operator per iteration of the linear solver. Note that g(w k ) can be stored and only has to be evaluated once per Newton iteration. The step sizes of the finite differences are chosen along the ideas from [1] and are given by with the machine accuracy machine and the free parameter scale .
To accelerate the convergence of the GMRES method, a preconditioner is applied. As the purpose of the preconditioner is to approximate the Jacobian dg(w)∕dw but being much simpler to invert, an elementwise block-Jacobi preconditioner is applied. It accounts for inner-element dependencies but neglects inter-element dependencies. Hence, it consists of one independently invertible matrix for each element. More details on the derivation of the block-Jacobi preconditioner for the mixed DG/FV sharp interface formulation can be found in [71].
Timestep restriction If no surface tension effects are considered, the timestep of the method is given by Eq. (7). Otherwise, the capillary timestep restriction [14] has to be considered as well. Taking this constraint additionally into account and following [14], the timestep is given by where the capillary eigenvalue can be calculated as with the surface tension coefficient and the characteristic length of the FV sub-cells l * ref = l ref ∕ (N + 1) . Note that the capillary timestep constraint only has to be evaluated at FV sub-cell DOFs directly adjacent to the surrogate surface . Differently to [14], we consider the total velocity of the flow at the interface |u| and not only the velocity component tangential to the phase interface.

Validation and Efficiency of the IMEX LSGFM
In this section, the novel IMEX LSGFM is validated and the performance of the scheme is investigated. In particular, the separation of fast and slow waves and the semi-implicit interface coupling is validated and the efficiency gain compared to a fully explicit scheme is evaluated.

Validation of Domain Coupling
First, the correct coupling of the fully implicit and the IMEX domain is validated for onedimensional settings. Acoustic wave interaction The first setup is taken from [15] and describes a single acoustic pulse which interacts with a phase boundary of water and air. With this setup, inconsistencies in the discretization, especially in the coupling of the different parts of the domain, would show up as they can have a large influence on amplitudes and propagation speeds of the occurring waves. The velocity of the phase boundary itself is small, such that no discrete interface movement takes place.
At We discretize the domain with 175 elements with N = 5 , couple them with the standard TV flux and choose the 3rd order IMEX-ARK3 scheme from [39] for the temporal discretization. As a comparison, the explicit scheme according to [37] with the explicit 3rd order low storage Runge-Kutta scheme from [69] is chosen.
According to [15], the ratio between the amplitude of the transmitted and the reflected acoustic wave is Δp trans. water ∕Δp refl. air = 2.001 . We perform calculations with CFL = 1 with respect to the acoustic eigenvalue and a simulation with the timestep increased by a factor of 100. The results of the simulations are depicted in Fig. 4.
If CFL = 1 is chosen, a comparison with the explicit scheme is possible. One can hardly distinguish the two solutions (upper part of Fig. 4) and both have a similar ratio of the transmitted and reflected waves Δp trans.
water ∕Δp refl. air = 2.016 , which matches the analytical solution very well. Thus, our novel method reproduces the expected behavior and predicts the correct dynamics, if the associated waves are resolved sufficiently. For the increased timestep no calculations with the explicit scheme are possible and deviations from the afore observed dynamics are to be expected for the IMEX method (lower part of Fig. 4). One can see that due to the large timestep the dissipation and dispersion of the waves are increased. Nevertheless, the ratio of the waves still matches the analytical solution with Δp trans. water ∕Δp refl. air = 2.007 with high accuracy.
Accelerated droplet The next setup, taken from [3], describes the acceleration of a liquid due to a pressure difference in the surrounding gas. Differently to the previous setup, the phase interface moves and thus allows further validation of the IMEX LSGFM.
On the domain = [0, 1 m] with Dirichlet boundary conditions, the initialization is given by The gaseous phase ( x 1 < 0.4 m or x 1 > 0. 6 m ) is modeled as an ideal gas with gas = 1.4 and the liquid phase as a stiffened gas, using liquid = 7.15 and p ∞ = 3 308 Pa . The discretization is done with 100 elements with N = 3 , TV flux and the 3rd order Runge-Kutta schemes from [39] and [69]. For the explicit scheme, the timestep is restricted with traveling waves inside of the droplet, the IMEX scheme predicts an almost linear pressure profile. This is to be expected as the IMEX scheme does not temporally resolve those waves, instead the linear pressure distribution of an incompressible material [3] is approximated.

Validation of Geometrical Coupling
Further validation of the novel scheme is provided by the simulation of the free droplet oscillation. This setup is particularly well suited for the validation of the consistent coupling of the different building blocks visualized in Fig. 2, as the whole dynamics is driven by the forces at the interface. Hence, a consistent coupling of interface properties, level set advection and the bulk phases is crucial. The test setup describes an ellipsoidal shaped droplet at rest in a periodic domain = [−2 m, 2 m] 3 . Due to surface tension, the droplet starts to oscillate. The initial data of the level set field is given by and the remaining initial data are summarized in Fig. 6, with denoting the surface tension coefficient.
The choice of the material parameters ensures that the ratios of the speeds of sounds are similar to the ones given in [57] for water and air. For this setup the maximum occurring Mach numbers are Ma ≈ 5 × 10 −3 and the analytical oscillation period for the lth mode can be calculated with the formula by Lamb [44] Φ 0 (x) = x 2 1.25 2 + with the radius of the droplet in equilibrium r equi . For the given initial data, the oscillation period of the second mode is T = 5 s. A grid convergence study is performed to validate the novel scheme. For that purpose, grids with 15 3 , 20 3 and 30 3 elements with N = 3 and TV or TVHLL-LM flux are used. The temporal discretization is done with the 4th order scheme IMEX-ARK4 from [39]. Instead of calculating the capillary timestep [14], the maximum timestep is prescribed with Δt max = 10 −3 s for the coarse meshes and Δt max = 7 × 10 −4 s for the finest mesh. Figure 7 shows the results in kinetic energy E kin = ∫ 1 2 ‖u‖ 2 2 d (left) and the mass loss (right). One can note that for a finer resolution the numerical viscosity is decreased and hence the damping of the oscillation is reduced, resulting in a smaller oscillation period. Considering the mass loss, also a convergent behavior can be observed, i.e., the mass loss decreases for a finer resolution. The figure shows that choosing the TVHLL-LM Riemann solver instead of the TV Riemann solver significantly improves the solution. While the mass loss is almost the same, much less numerical dissipation is introduced, resulting in a more accurate solution.

3
A quantitative analysis of the convergence study is provided in Table 1. One can see that the numerical solution approaches the analytical solution of T = 5 s and the mass loss decreases under grid refinement.
Summing up, the presented studies validate the novel scheme. Investigations on the efficiency are provided in the next section.

Influence of the Approximation of the Jacobi-Vector Product
Choosing implicit time discretization opens the field of optimization of the (non-) linear solver parameters. One of the DOFs of the implicit method is the definition of the finite difference, approximating the Jacobian-vector product of the GMRES method, see Eq. (10). Here, one can choose a first or a second order finite difference and the step size of the finite difference (Eq. (11)). In literature, several options how to choose the step size can be found in Eq. (12) with n DOF denoting the total number of DOFs. For the present scheme, the simple choice (1) scale ∶= p 0 air , denoting the background pressure in the gaseous domain, is proposed. The choice of the scaling factor is a trade-off between truncation error of the finite difference and round-off errors due to limited machine accuracy.
A comparison between the different options is provided by the analysis of the required iterations for a certain amount of timesteps. For that purpose, the setup from the previous Sect. 4.2 with 20 3 elements, N = 3 and TV flux is used. The amount of iterations are counted to advance the solution from t = 2 s for 50 timesteps with IMEX-ARK4 and Δt = 10 −3 s each. The relative convergence tolerance of Newton's method is set to Newton = 10 −4 and the relative convergence tolerance of the linear solver GMRES is varied. Two sets of simulations are executed: one with the parameters summarized in Fig. 6 and one with an increased background pressure of p 0 air = 10 5 Pa , p 0 water = 1.004 567 × 10 5 Pa and p ∞ = 5.7 × 10 −8 Pa , decreasing the occurring Mach numbers by a factor of 10. This allows an evaluation of the influence of the Mach number on the iterations. The results of the simulations are shown in Fig. 8.
scale ∶= One can see that the choice of scale strongly influences the required amount of iterations. For the larger Mach number almost no differences between the different options are visible, except for the choices (4) scale and (6) scale , which require more iterations. In those two cases, the second order finite difference decreases the number of iterations. For the smaller Mach number, it is only possible to obtain a solution with (1),(2), (5) scale . With the remaining options, convergence of the linear solver was not obtained with a maximum of 500 iterations per Newton step. From Fig. 8 one can conclude that with (1) scale or (2) scale a second order finite difference is not necessary for both considered setups. As the calculation of (1) scale does not require any further computational effort, it is chosen for the following computations.

Comparison with Fully Explicit Scheme
Finally, a comparison with a fully explicit scheme regarding the required computational time is presented. With this comparison, the influence of the block-Jacobi preconditioner and the chosen timestep size on the computational resources is evaluated.
Again, the setup from Sec. 4.2 is taken with Newton = 10 −3 , GMRES = 10 −1 and the CFL number is varied from CFL = 8 × 10 −4 to CFL = 8 × 10 −2 . Note that the CFL number is calculated with the maximum convective eigenvalue ̂ max = |u ⋅ n| max . The explicit reference scheme uses a 4th order low storage Runge-Kutta scheme from [6] for the time discretization with CFL = 0.9 with respect to the acoustic eigenvalue. A comparison of the computational results is provided in Fig. 9.
It is clearly visible that the simulation results with the different CFL numbers and the explicit schemes are almost not distinguishable. Due to the similar solution quality those setups can be used to evaluate the influence of the timestep size on the required GMRES iterations and computational time.
For that purpose, five simulations from t = 1.99 s to t = 2.04 s are executed for each configuration on one node of the HPE "Hawk" system with 128 cores and the mean value is taken for the comparison. This choice leads to approximately 4 000 DOFs per core, what has turned out to be a good compromise between processor local work, communication and required memory. First, the iterations of the linear solver are considered. Figure 10 (left) shows that the preconditioner has only a small influence on the convergence properties of the linear solver for small timesteps. For larger timesteps the block-Jacobi preconditioner allows to reduce the required amount of iterations significantly. While the required iterations are almost constant without preconditioner, the iterations are reduced for larger timesteps if the preconditioner is chosen. The amount of evaluations of the spatial operator for the explicit scheme is marked with blue. One can see, that with the IMEX scheme with large timesteps up to a factor of ≈ 2 less operator evaluations than with the explicit scheme are necessary.
Considering the required wall-clock time as a measure of efficiency shows that large speed-ups can be obtained using the novel IMEX scheme. Figure 10 (right) shows that a speed-up up to a factor of approximately 25 can be obtained compared to the explicit scheme. This large gain in wall-clock time is mainly due to the reduced amount of expensive operations that have to be executed in each timestep: velocity extrapolation, the evaluation of the interface geometry (i.e., Φ and Φ ) are done once per timestep and the level set reinitialization is done only each second timestep. As those operations have a significant  influence on the overall computational time, large timesteps are a key feature to reduce the computational costs. Additionally, one can see that also here the use of the block-Jacobi preconditioner is beneficial for large timesteps. Although for small timesteps the computational time is increased due to the costs of building and applying the preconditioner, this is counterbalanced for large timesteps by the fewer iterations.
Concluding, the investigations in this section have shown that the novel IMEX LSGFM allows for an efficient and accurate simulation of low Mach number droplet dynamic phenomena. Compared to a fully explicit scheme, large gains in computational time can be achieved if large timesteps are chosen.

Application of the IMEX LSGFM
In this section, the novel scheme is applied to two more typical applications of low Mach number droplet dynamics to illustrate its capabilities.

Drop Impingement onto a Deep Pool
A typical field of low Mach number droplet dynamics is the impingement of droplets onto pools of various depth. In [58] a two dimensional numerical study on the influence of the Weber number, pool depth and impingement angle has been carried out. We utilize one specific setup of this study to illustrate and validate the capabilities of the IMEX LSGFM.
The impingement Weber number, defined as with radius r, density and velocity of the droplet of the simulation setup is We = 217 . Gravity with the gravitation constant g = 9.81 m ⋅ s −2 is taken into account. The initial data and the boundary conditions of the computational domain are summarized in Fig. 11. The domain = [−20r, 20r] × [−5r, 15r] is discretized with 400 × 200 elements with N = 2 , which are coupled with the TVHLL-LM flux. The pool surface is located at x = 0 m and has a depth of 5r and the droplet center is initialized 2r above the pool surface. With the parameters summarized in Fig. 11, an initial Mach number of approximately Ma = ‖u 0 ‖ 2drop ∕c air = 1.7 × 10 −2 can be calculated. The characteristic momentum ( u) * ∶= ‖ u‖ 2 ∕‖( 0 u 0 ) water ‖ 2 and the characteristic timescale * ∶= 2r∕‖u 0 ‖ 2 are used to relate the occurring quantities to the initial conditions. The temporal discretization is done with the IMEX-ARK3 [39]. For stability reasons, the CFL number had to be reduced to CFL = 0.1 for this example. Results of the simulation are shown in Fig. 12.
The figure shows that when the droplet hits the pool surface a jet is formed; in three dimensions named crown. Following, the crown height and the crater depth increase and a satellite droplet is ejected. This has also been observed in [58]. At t = 13.5 * = 0.99 ms the crater depth, the crown height and the angle of the crown are measured. A comparison with the results from [58], provided in Table 2, shows a good agreement. Deviations, especially in the crater depth and crown height are most likely due to the neglected viscous effects. Moreover, deviations caused by different numerical resolutions are possible as the number of DOFs is not reported in [58].

Binary Droplet Collision
A comparison with experimental data is possible for a three dimensional setup, where viscosity plays a minor role. This is the case for the binary head on collisions in [2] with a collision Weber number of We = 23 . For this Weber number one can expect the coalescence of the droplets followed by a separation without ejecting further satellite droplets [2], see [63] for incompressible simulations of this setup. The initial data are summarized in Fig. 13. The domain = [0, 1 m] 3 is discretized with 30 3 elements with N = 2 and describes one eighth of the total domain. At the boundaries, symmetry boundary conditions are prescribed and the TVHLL-LM flux is chosen as numerical flux function. With the parameters summarized in Fig. 13, an initial Mach number of Ma = ‖u 0 ‖ 2water ∕c air = 5.7 × 10 −2 can be calculated. The temporal discretization is done with the IMEX-ARK3 scheme [39] and CFL = 0.9 . The results of the simulation are reported in Fig. 14, together with the experimental results from [2]. One can see that the experimental results are reproduced with high accuracy: first the droplets merge and form a torus; due to the high surface tension they then elongate and finally separate again. All different type of shapes are correctly captured with the IMEX LSGFM.

Conclusion and Outlook
In this work it has been shown how an explicit LSGFM can be extended to deal with low Mach number flows. To increase the prohibitive small timestep of the explicit method, an IMEX flux splitting technique has been applied. It has been shown that by separating the domain into a fully implicit narrow band around the phase interface and into an IMEX flux splitting region, the timestep can be increased significantly. The separating of the domain became necessary as the interface flux and the narrow band region cannot be treated partially implicit and partially explicit. A mixed discontinuous Galerkin/finite volume method has been used for the spatial discretization and IMEX Runge-Kutta schemes with a Jacobian-free Newton-GMRES method have been used for the temporal discretization.
The novel IMEX level set ghost fluid method has been validated and it has been shown that a large gain in computational time can be achieved compared to a fully explicit method. Applications to typical droplet dynamic problems illustrate the broad applicability of the developed scheme.  with the interface density, velocity and pressure with and A is given by Eq. (6). In the fully implicit region, the numerical flux is simply the sum of both parts TVHLL-LM Riemann solver For small Mach numbers, the low dissipation TVHLL-LM solver can be used. A low dissipation numerical flux for low Mach numbers is obtained by modifying the HLL flux from [65] with the low Mach number fix according to [9]. For the explicit system the original flux of the TV flux is used, see Eq. (13). The implicit numerical flux is given by The approximate wave speeds S L and S R , which are required for the HLL scheme can be calculated with with the wave speed estimates