Conservation of energy in the direct numerical simulation of interface-resolved multiphase flows

In the context of deformable bubbles, surface tension produces a dynamic exchange between kinetic and surface elastic energy. This exchange of energy is relevant to bubble dynamics, like bubble induced turbulence or drag reduction. Unfortunately, the underlying physical mechanism is not exactly explained by the state-of-the-art numerical methods. In particular, the numerical violation of energy conservation results in an uncontrolled evolution of the system and yields well-known numerical pathologies. To remedy these problems, we tackle two of the most problematic terms in the numerical formulation: convection and surface tension. We identify the key mathematical identities that imply both physical conservation and numerical stability, present a semi-discretization of the problem that is fully energy preserving, and assess their stability in terms of the discrete energy contributions. Numerical experiments showcase the stability of the method and its energy evolution for stagnant and oscillating inviscid bubbles. Results show robust as well as bounded dynamics of the system, representing the expected physical mechanism.


Introduction
Turbulent bubbly flows appear often in engineering applications, like bubble columns or bubble-laden pipes, and lead to many interesting physical phenomena that need further understanding. Two of the most practical applications of bubbly flows are drag reduction and chemical selectivity. Since its observation in the 1970s, drag reduction has attracted significant research enthusiasm, mostly motivated by naval applications. Promising reductions in drag being as high as 80% (Madavan et al., 1984), it could dramatically reduce ship's fuel consumption and emissions. While a comprehensive understanding of the physical mechanisms behind it is still missing (Lohse, 2018;Mathai et al., 2020), recent studies (Spandan et al., 2018) suggest that deformation of the bubbles is responsible for the significant drag reduction, in particular in the bubble wake. Another remarkable flow feature is the intensification of local micro-processes in the chemical industry (Schlüter et al., 2021), for which a better understanding of local transport phenomena in the bubble wake is also required.
Due to the delicate physical setup of the bubbly flows and the small scale phenomena involved in such processes, numerical simulations are an essential tool for the development of these technologies. However, the simulation of multiphase flows is prone to several numerical pathologies. These pathologies arise both from the difference in physical properties and the inclusion of surface tension, and manifest themselves by the appearance of parasitic currents, the miscalculation of bubble deformation, and ultimately the numerical instabilities. These issues not only jeopardize the viability of the simulation, but even when able to converge, they may misrepresent important flow physics, such as the exchange of surface and turbulent kinetic energy between carrier and dispersed phases (Elghobashi, 2019) or the contribution of surface tensions to vorticity in deformed bubbles (Hasslberger et al., 2018). The remediation of instabilities involves mostly the inclusion of stabilization methods, which adds numerical diffusion to stabilize the system. However, this approach may ultimately compromise the physical reliability of the simulation (Trautner et al., 2021). Vol. 5, No. 4, 2023, 333-343 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-022-0148-4 N. Valle,R. Verstappen 334 To tackle these numerical problems, a method for the direct numerical simulation (DNS) of multiphase flows is presented. We believe that, since the physics of the problem are stable, numerical methods should be stable, too, by mimicking the same physical mechanisms of the original system without adding any extra non-physical terms. By identifying the key mathematical identities which are responsible for the conservation of the physical variables of interest in the discrete setting, the method mimics most of them to ensure the conservation of basic physical quantities.
It turns out that most standard numerical methods for multiphase flows do not preserve conserved quantities at the discrete level. Understanding energy as a weighted norm of the variables of the system, energy preserving discretization can also be seen as a method producing inherently bounded quantities, which directly translates to stability.
Departing from the conservative level set (CLS) method (Olsson and Kreiss, 2005) for interface capturing, novel numerical techniques are introduced to enforce conservation of primary and secondary physical properties. The transport of variable density flows is approached with a consistent mass and momentum transfer scheme (Rudman, 1998), which also results in the conservation of kinetic energy (Mirjalili and Mani, 2021). Surface tension is approached with the novel energy preserving scheme (Valle et al., 2020), which substantially removes parasitic currents.
By adopting setting a fully energy-preserving discretization in this work, we assess its performance in several canonical cases as a first step towards the full DNS of bubble induced turbulence. The method capabilities are stressed in the simulation of a canonical flow configuration. Focus is put on the stability and the conservation of both primary and secondary physical quantities.

Formulation
The bubble interface is captured by means of a classical CLS method (Olsson and Kreiss, 2005), which introduces a smoothed volume fraction θ. Such a volume fraction is considered as the convolution of the distance function: where d is the signed distance function used in a classical level set and ε determines the thickness of the regularized interface. Such a thickness is a function of the mesh spacing (Olsson and Kreiss, 2005). The limit of 0 ε  is the Heaviside step function. An exemplification of the regularization process can be seen in Fig. 1.
The volume fraction gradient is used to compute the surface area, the interface normal 1 n θ θ -=|  |  , and the curvature κ n = ⋅ (Olsson and Kreiss, 2005). Assuming an incompressible flow: where u  stands for the velocity, the transport of θ in conservative form reads as Adopting a one-fluid model, both liquid and gas phases are modeled as a single fluid whose physical properties are calculated from θ . Local density ρ and viscosity μ are computed from the marker function θ and the properties of each phase, i ρ and i μ , as Owing to the linear density variation with θ in Eq. (4) and the conservative nature of Eq. (3), the conservation of mass is implicit. The dimensionless form of the inviscid Navier-Stokes equations introduces the Weber number b We Ru g γ =  , where R and b u  are the bubble's radius and velocity, respectively, g is gravity, and γ is the surface tension; and finally reads Here, P stands for the pressure, and S is the strain tensor. Note that the stress discontinuity at the phase interface, resulting from surface tension phenomena, is introduced in the formulation by means of the continuum surface force (CSF) (Brackbill et al., 1992). This results in a fully conservative discretization of mass and momentum, while energy should be preserved owing to the inviscid flow.

Energy-preserving discretization
The formulation described above is brought into a computational setting by mimicking each of their differential operators in a staggered, second-order finite volume framework.
Discrete fields are arranged as vectors, and discrete differential operators adopt the form of matrices. Collocated quantities are located at the cell centers (e.g., c θ , c p ) and staggered ones at the face centers (e.g., f u , f ρ , f n ). Both collocated and staggered arrangements define vector spaces that can be equipped with an inner product ( ) ⋅,⋅ . This applies to both continuum (( , ) d ) f g f g V = ⋅ ò and discrete T (( , ) ) M = f g f g fields. When needed, suitable interpolators can be constructed from cells to faces (Π) and vice versa (Π ) * . These may include high-order interpolations.
The rest of the operators are: G, the cell-to-face discrete gradient; D, the face-to-cell discrete divergence; f L , the staggered diffusion, which already includes viscosity; while C c ( ) ⋅ and C f ( ) ⋅ are the collocated and staggered convective terms, which depend on a flow field. For a detailed discussion on the construction of such operators, see Valle et al. (2020). In this framework, the semi-discretized governing equations (2), (3) and (5) read This formulation is prone to numerical instabilities. Two of the most problematic terms are momentum convection f f ( ( )) C ρ u u (discussed in Section 2.2.1) and surface tension G f c ( ) k θ (discussed in Section 2.2.2).

Convective term
The instability from the discretization of the convective term typically arises from an inconsistent mass and momentum transport. This has been previously reported in the context of high density ratios (Rudman, 1998). The same results have been presented from an energy-preserving approach (Veldman, 2019;Mirjalili and Mani, 2021).
As a summary, we present the energy-preserving discretization of the convective term reported by Mirjalili and Mani (2021), where proof of energy-conservation and stability can be found.
Disregarding viscosity and surface tension from Eq. (5), the conservation of kinetic energy ( , ) ρu u   can be expressed in terms of conserved variables ρ and ρu  as From the continuum point of view, the idea is that the contribution of mass transport to kinetic energy needs to be compensated with the contribution of momentum transport to kinetic energy.
At the discrete level, the semi-discretized form of Eq. (9) requires bringing density from cells to faces: where  represents element-wise multiplication. From a numerical perspective, the discrete energy Equation (10) then states that the evolution of k E is bounded and the values of f u will remain finite, and thus stable, throughout the simulation.
For the contribution of mass transfer to Eq. (10), we combine Eqs. (4) and (7) to obtain Eq. (11): which implies that the mass transport equation uses the same advection scheme as the marker function. This defines the density at the face f ρ in terms of the approximation of c θ at the face employed by c f ( ) C u in its transport.
The key in consistent mass and momentum transport lies in matching the (staggered) mass flow implied by c f ( ) C u with the definition of the (staggered) momentum. Note that both mass flux and momentum actually have the same form: ρu  .
For the contribution of momentum to Eq. (10), we introduce Eq. (8) to obtain which represents the contribution of the convective term to kinetic energy. This term is known to be neglected for incompressible and constant density flows owing to the use of symmetry-preserving schemes (Verstappen and Veldman, 2003). When density is variable, the use of the consistent mass and momentum transport described above, in conjunction with the well-known symmetry-preserving scheme (Verstappen and Veldman, 2003), allows to reformulate f f ( ) C ρu as (Mirjalili and Mani, 2021;Valle et al., 2021): (Mirjalili and Mani, 2021), which is well-known to produce a null contribution to energy (Verstappen and Veldman, 2003), and so this term vanishes.
Finally, substituting Eqs. (11) and (12) into Eq. (10) results in a null contribution to kinetic energy at the discrete level, since the non-null contributions, due to density variations, cancel each other as far as the mass and momentum fluxes are consistent.

Surface tension
For surface tension, we adopt the formulation presented in Valle et al. (2020), which balances kinetic energy with surface (potential) energy variations. Surface energy arises from the tensional state at which the interface surface is subject to, which is due to surface tension (Isachenko et al., 1974) and is thus proportional to the surface area of the interface, which we denote by Γ .
The key idea in Valle et al. (2020) is to satisfy the first variation of area (Frankel, 2011): at the discrete level. Such a geometric identity enforces the match between kinetic and surface (potential) energy transfers. By disregarding viscosity in Eq. (5), the evolution of kinetic and surface energies read as (Valle et al., 2020): At the discrete level, Eq. (15) requires to define the discrete potential energy: From a numerical point of view, p E can be seen as a norm of c θ . Thus, the discrete form of Eq. (15) can be seen as stating that the combined evolution of both f n and c θ remains bounded, and thus the simulation is stable.
Finally, Eq. (17) is enforced by constructing the appropriate cell-to-face curvature interpolator  . This is required since curvature is naturally computed at cell centers, c f D = k n (Olsson and Kreiss, 2005), whereas in a staggered formulation curvature is required at the faces. In the context of a CLS-CSF method, we obtain a discrete counterpart of the first variation of area as (Valle et al., 2020): Introducing Eq. (11), and exploiting the duality between G and D of the symmetry-preserving formulation (Verstappen and Velman, 2003), we can rearrange as n . Finally, after taking the adjoint of the left hand side, we obtain From where we can define  in terms of the cell-to-face interpolation of c θ used for the construction of the fluxes in f c (Valle et al., 2020). In summary, any upwinding used for interpolation of c θ corresponds with a reciprocal downwinding for the interpolation of c k (Valle et al., 2020).

Numerical method
In the transport of the marker function, we adopt a classic high resolution Superbee flux limiter for the spatial discretization, while temporal integration is performed with a Runge-Kutta integrator.
We have omitted the recompression step (Olsson and Kreiss, 2005;Olsson et al., 2007) of the original CLS formulation, which includes an additional integration in pseudo-time τ : This is entirely artificial compression aimed at maintaining a sharpness of thickness proportional to ε , which is otherwise spread as a result of the transport scheme (Olsson and Kreiss, 2005). However, the change of θ in pseudo-time does also change the amount of surface area, which in turn introduces an uncontrolled amount of energy to the system (Valle et al., 2020). We acknowledge that this does compromise the sharpness of the interface, particularly in the presence of high shear flows, but we have given precedence to energy consistency over interface sharpness. Further discussion about this particular issue is addressed in Section 4.
Pressure-velocity decoupling is achieved by a fractional step method (Chorin, 1968). Contrary to Mirjalili and Mani (2021), in this work we do impose incompressibility of the velocity field rather than of the momentum field. It results in a variable-coefficient pressure Poisson equation, which reduces the set of algebraic solvers suitable for it. In this work, we employ a conjugate gradient method with a Jacobi preconditioner, which aims at mitigating the stiff condition number, particularly for regions away from the interface.

Results
In order to test the robustness of the presented energypreserving discretizations, the following canonical cases are presented for validation.
1 We = was assumed for simplicity.

Stagnant bubble
A bubble of diameter 0 15 R = . was placed at the center of a fully periodic box [1 1 1]´, which was represented by a uniform cartesian mesh of [64 64 64]´ cells. Viscosity and gravity were set to 0. The exact analytical solution predicts a fully stagnant velocity field. For a bubble of the density 1 ρ , surrounded by a liquid of the density 0 ρ , linear perturbation theory predicts the oscillation period on an ellipsoid as (Lamb, 1945): ( 1)( 2)

R s ρ sρ T γs s s
which was used to present the time evolution results. This case is known to exhibit spurious currents, which arise form errors incurred into the calculation of curvature κ (Magnini et al., 2016). Most state-of-the-art methods typically diverge in the absence of viscosity after a few iterations.

Unit density ratio
As a first test, density ratio was set to 1 in order to remove the effects of the convective term and thus isolate the impact of the surface tension force into the dynamics of the system.
Results show the inception of a non-zero velocity field at the beginning of the simulation, which nonetheless vanishes as the interface is accommodated to the grid resolution. This is remarkable since no diffusion is present in the system. Instead, the only mechanism taming the spurious currents in the system is the exchange between kinetic and potential energy (Valle et al., 2020). As it can be seen in Fig. 3, the total energy is effectively preserved.
This initial velocity field arises from the mismatch between the initialization of the bubble profile and the numerical equilibrium one, which results from the minimization of surface energy. While sphere is the shape that already minimizes the surface potential energy, this does not necessarily hold at the discrete level. From a numerical perspective, the initialization of the bubble profile may inadvertently introduce additional elastic energy in the bubble, which provides the impulse to the initial velocity field. Nonetheless, owing to the periodic domain and the conservative discretization, the system is quickly brought to a new equilibrium state by accommodating the bubble interface to the new minimal energy. This results in the vanishing of the initial velocity field and the corresponding adjustment of the bubble's interface, which can be seen from the kinetic energy evolution in Fig. 3.
Interestingly enough, the currents produced in the vicinity of the bubble in the geometric adjustment described above appear in a ring region around the bubble, which expands and fades away as the system evolves in time in a "shock"-like fashion.

Non-unit density ratio
A more challenging test involves assuming a density ratio different from 1. Since we are interested in the rising by buoyancy of lighter bubbles, we will focus on density ratios 1 0 1 ρ ρ  / . The method results in a solid bubble representation, which retains its shape throughout the simulation, even for longer simulation periods.
As in the unit density ratio case, spurious currents appear at the beginning of the simulation and fade away after a short period. Due to the high density ratios, these velocities achieve higher values than in the unit density case, since the numerical imbalance is introduced in the momentum equation. As it can be seen from the right halves of Fig. 4, this results in a higher magnitude of the velocity field in the region inside the bubble. Nonetheless, they remain bounded within reasonable limits owing to the energy-preserving formulation.
Contrary to the unit density ratio case, the structure of the currents does not show anymore the ring-like shape, but rather exhibits a more complex pattern which retains connectivity with the bubble interface. This can be explained by the smooth transition zone around the bubble customary of the CLS, which produces a smooth density gradient. Such a density gradient dampens the "shock"-like evolution of the spurious currents that was observed in the unit density ratio case. Nonetheless, the structures still fade away as the system evolves, even though its dynamics are slower, as it can be seen from Fig. 4.
High density ratios push the numerical solutions to the limit, since small perturbation in momentum may produce huge accelerations in lighter fluid due to the 1 ρ / factor in the pressure Poisson equation. Nonetheless, results in Fig. 5 show a robust method and a bounded energy evolution. Even when spurious currents are an order of magnitude higher compared with the unit density ratio, the system is remarkably stable considering that the density ratio is three orders of magnitude smaller. The different density ratios introduce a new source of error in this configuration. Due to the variable density, errors of the pressure Poisson equation are unevenly scaled, and are actually amplified at the lighter phase.

Oscillating bubble
Now, the original bubble is slightly deformed to turn into the shape of an ellipsoid. For small perturbations, the exact analytical solution predicts an oscillatory ellipsoid that evolves from prolate to oblate and vice versa (Lamb, 1945). The ellipsoid is stretched by 20% on the y axis. Using Eq. (21), we can estimate the oscillation period for small deformations. Due to the dynamic equilibrium solution, this test introduces the complexity of the dynamic interface capturing. This is the reason why the mesh has been refined to [96 96 96]´ cells.

Unit density ratio
For a unitary density ratio, the convective term should not introduce new difficulties associated with the evolution of the variable density flow, which allows showcasing the role of the surface tension. Results in Fig. 6 show the first oscillation of the ellipsoid, from where good qualitative interface representation can be observed.
Looking at the energy evolution in Fig. 7, we observe a bounded, and thus stable, time evolution of the system's energy. Even after an initial energy peak, which we link with the adjustment of the interface representation described in the case of the stagnant bubble, we see how the system can well capture the oscillatory behavior of the system.
However, an artificial damping of the oscillatory solution is observed. This can be explained by several mechanisms: First, while the space discretization technique is conservative, as it can be seen from Fig. 7(a), the time integration scheme is not. A second mechanism may lie in the degradation of the quality of the interface representation, which may broaden the interface artificially and eventually misrepresent the interface location. This is expected, since the energy-preserving formulation employed here does not include any additional  recompression steps, as they result in an artificial increase of the mechanical energy (Valle et al., 2020). As a result, an artificial thickening of the interface is expected, which brings the system to a new stagnant equilibrium in the long run. Finally, the lack of linear momentum conservation, which is so far elusive, may produce inaccurate transfers of  Conservation of energy in the direct numerical simulation of interface-resolved multiphase flows 341 potential and kinetic energies. Even when those transfers are enforced to be consistent, preserving the semi-discretized mechanical energy, the actual intensity of those transfers may not be accurate, resulting in an inaccurate share of potential and kinetic energies.

Non-unit density ratio
Now, stressing even more the previous case, we set a density ratio of 3 1 0 10 ρ ρ -/ = . As in the case of the stagnant bubble, high density ratios amplify the numerical errors coming from the momentum equation which will go into the velocity field. This results in higher velocity magnitudes in the interior of the bubble, as can be seen by comparing Figs. 8(b) and 8(e), which show nearly symmetric geometric states. However, velocities are way higher inside the bubble. Once again, results show a bounded evolution of the system's energy, even for high density ratios. Even when the system is energetically bounded, and thus stable, we observe a frequency degradation on the oscillating behavior, which may be due to the mechanisms discussed above for the stagnant case, and intensified by the high density ratio at the light phase.

Conclusions
A full energy-preserving method has been introduced for the calculation of several canonical flow configurations. This is achieved by identifying the corresponding mathematical identities of the governing equations at the continuum level and mimicking them at the discrete one. This results in the following: For the discretization of the convective form, the method imposes the same discretization for mass flux and momentum transport, since they are inherently the same quantity ρu  . For the surface tension term, the method links the cell-to-face interpolation of the curvature with the cell-to-face interpolator used for the transport of the marker function. This work highlights how energy-preserving methods not only mimic the features of the actual physical system, but also achieve stability by the same mechanisms of the physical setup. While stability is guaranteed by enforcing the transfers between kinetic and elastic energies to match each other exactly, there is no guarantee that the quantity of such a transfer is correct. We conclude that this requires two main conditions: an accurate interface representation and a conservative surface tension discretization. The accuracy of the interface representation depends on the transport scheme for the marker function. Equation (14), and its discrete counterpart Eq. (19), are linear in velocity (Valle et al., 2020), and so should the transport scheme be linear in velocity as well. This is satisfied by flux-limited schemes, but not by other more sophisticated transport schemes, such as the original CLS (Olsson and Kreiss, 2005). In physical terms, this corresponds with the fact that any deformation of the interface that is not due to a velocity field is inherently non-physical. Consequently, unless special care is taken to produce a null energy contribution in this artificial step, the conservation of energy is seriously compromised. This latter result suggests the development of either energetically neutral recompression steps or new physics-compatible high resolution schemes, potentially including recompressive and anti-diffusive fluxes.
The conservative discretization of surface tension is not guaranteed by using the CSF (Popinet, 2018), a property which closed surfaces should satisfy exactly (Blackmore and Ting, 1985). The reason behind this lack of trivial conservation lies in the embedding of the interface into the background Eulerian mesh, which results in the lack of an explicit and closed representation of the interface surface. This implies that while the transfers between the kinetic and surface energies are balanced owing to the energy-preserving discretization, such a transfer may not be physically realistic.
Finally, since in the development of the energy-preserving discretization presented in Section 2 we have worked with the semi-discretized (i.e., in space only) form of the equations, conservation errors introduced by discrete time integration have not been considered. These could be removed by the use of symplectic time integrators, but since those are necessarily implicit, its computational cost is prohibitive. For this reason, as a future task, we may consider the inclusion of pseudo-symplectic Runge-Kutta methods for the temporal integration of the system of equations (Capuano et al., 2017), which are aimed at reducing errors in conservative quantities introduced by discrete time integration.