{'fraction_non_alphanumeric': 0.036632747456059205, 'fraction_numerical': 0.012617946345975949, 'mean_word_length': 4.601243523316062, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A fully conservative sharp-interface method for compressible mulitphase flows with phase change

A fully conservative sharp-interface method is developed for multiphase flows with phase change. The coupling between two phases is implemented via introducing the interfacial fluxes, which are obtained by solving a general Riemann problem with phase change. A novel four-wave model is proposed to obtain an approximate Riemann solution, which simplifies the eight-dimensional roo-finding procedure in the exact solver to a sole iteration of the mass flux. Unlike in the previous research, the jump conditions of all waves are imposed strictly in the present approximate Riemann solver so that conservation is guaranteed. Different choices of the fluid states used in the phase change model are compared, and we have shown that the adjacent states of phase interface should be used to ensure numerical consistency. To the authors' knowledge, it has not been reported before in the open literature. With good agreements, various numerical examples are considered to validate the present method by comparing the results against the exact solutions or the previous simulations.


Introduction
The multiphase flows with phase change are ubiquitous in scientific and industrial applications such as cavitation in the naval industry [1], bubble condensation in electronic cooling [2], the combustion of fuel in rocket engines [3], etc. Over the last decades, many numerical methods have been developed to investigate the underlying mechanisms of these flow phenomena, which can be generally divided into two classes, i.e. diffuse interface methods and sharp interface methods. In the diffuse interface methods, the material interface is smeared out over a finite number of grid cells where different fluids are mixed [4]. For the mixed fluids at the interface region, a mixture model or an artificial equation of state (EOS) needs to be developed to obtain a thermodynamically consistent description [5], which is difficult as the mixture states may be even unphysical. A typical example is the homogeneous equilibrium method (HEM) proposed by Menikoff and Plohr [6], in which the mixed fluids were assumed to be a mixture of liquid and vapor in thermodynamic phase equilibrium. This approach was further extended to the van der Waals EOS by Müller and Voß [7] via using the Maxwell construction. However, the phase equilibrium assumption suffers a kink at the coexistence curves of two phases [8]. The seven-equation model [9][10][11] is a full non-equilibrium model, in which instantaneous relaxation procedures are imposed for pressure and velocity. By including temperature and chemical potential relaxation effects, Zein et al. [12] modified this model to take into account phase transition. Nevertheless, the solutions of the seven-equation model usually exhibit a complex wave structure [13]. In addition, in the diffuse interface methods, the interface profile tends to be thickened over time due to numerical diffusion. Although several interface sharpening schemes have been developed to suppress the interface diffusion [14][15][16], it remains a big challenge to conduct long-time simulations with the diffuse interface methods [17].
In the sharp-interface methods, the interface is modeled as a discontinuity within the flow filed [18][19][20], at which appropriate jump conditions need to be imposed to ensure the physical consistency of the overall model. To track the interface with a non-smearing representation, different approaches have been developed, e.g., the front-tracking method [21,22], where the interface is represented by massless markers, the arbitrary Lagrangian Eulerian (ALE) method [23,24], where the interface coincides with mesh lines, and the level-set method [25], where the interface is defined as the zero level-set of a signed distance function. When phase change is considered, although there have been many efforts in the context of incompressible multiphase flows (see Refs. [26][27][28][29][30], to name but a few), only a few numerical schemes for compressible regimes are reported in the open literature and all of them are based on the level-set method. In contrast to other approaches, the level-set method is simple to implement and can handle large interface deformations and topological changes in an automatical way [31], which make it more suitable for modelling phase change in compressible multiphase flows. In Ref. [32], the level-set method was combined with a ghost fluid method (GFM) [33] to investigate the compressible bubble growth with phase change. The states of ghost cells, i.e., the cells defined in the vicinity of the interface where the ghost fluid and real fluid co-exist, were determined by the jump conditions accounting for phase transfer. However, as shown by Liu et al. [34], the ghost cells should be populated according to the solution of an interfacial Riemann solver so that waves can be correctly transmitted and reflected at the interface. Following this strategy, Houim and Kuo [35] extended the GFM to reacting flows with phase change by solving a modified interfacial Riemann problem. The resulting velocity and pressure jumps due to phase change were incorporated in the conservation equations of the gas-liquid interface.
This work was further modified by Das and Udaykumar [36] via simplifying the implementation of the interfacial jump conditions. The numerical instabilities caused by the rotation of the deviatoric stress tensor were avoided by rotating the velocity field. Since the velocity and pressure jumps are calculated using the left and right states of the interfacial Riemann problem, one common issue in Refs. [35,36] is that the conservation equations are not truly fulfilled at the interface. By introducing an additional wave representing a phase interface, which was observed in several experiments [37][38][39], Fechter et al. proposed an exact Riemann solver [40] and an approximate solver [41] for a general two-phase Riemann problem including phase change and surface tension. Note that, in Refs. [40,41], the coupling between two phases at the interface is also imposed through the GFM.
However, for all the GFM-based methods, an unavoidable problem is the lack of conservation properties, which is crucial for compressible multiphase flows with strong interface interactions. For compressible multiphase flows without phase change, a fully conservative sharp-interface method was pro-posed by Hu et al. [42]. This method was extended by Lauer et al. [43] to simulate oscillating bubbles with phase change. They assumed that phase change is much slower than the hydrodynamic interface interaction, which may not hold for other cases. A more general model taking account of evaporation is independently developed by Paula et al. [44]. Since the heat flux at the phase boundary is neglected in Ref. [44], this model will lead to inconsistent stationary solutions, see Ref. [40] for a detailed derivation.
In this paper, we aim to develop a fully conservative sharp-interface method based on [42] for compressible multiphase flows with phase change.
Unlike in Refs. [43,44], we obtain interfacial fluxes through solving the general Riemann problem proposed by Fechter et al. [40] to ensure thermodynamic consistency. The exact Riemann solver proposed in Ref. [40] suffers from poor efficiency due to the resulting eight-dimensional root-finding problem. To address this issue, a novel approximate Riemann solver is developed in this paper. Unlike the approximate Riemann solver of fechter et al. [41], the present solver employs a four-wave model so that the energy coupling in the vapor phase is achieved. Since the specific form of an EOS is not used in this Riemann solver, it can be used to simulate real fluids. In addition, the choices of the liquid and vapor states used in the phase change model is investigated in detail, which is not discussed in the previous researches. As we will show in Section. 4.1, the intermediate states in the Riemann solution fan should be used to calculate the mass flux to ensure numerical consistency. The remainder of this paper is organized as follows. In Section 2, we briefly review the conservative sharp-interface framework for multiphase flows without phase change. Subsequently in Section 3 we introduce the novel approximate Riemann solver for the general Riemann problem including phase change and surface tension effects. A number of numerical tests are carried out to validate the present method in Section 4, followed by the concluding remarks in Section 5.

Governing equations
For inviscid compressible flows, the governing equations are where U = (ρ, ρu, ρv, ρw, ρE) T is the vector of conserved variables, in which ρ, u, v, w and E denote the density, the three Cartesian velocity components and the total energy with relation E = e + 1 2 (u 2 + v 2 + w 2 ), respectively. Note that e stands for the internal energy per unit mass. The inviscid flux tensor F reads and p is the pressure. To close this system, an EOS is needed to describe the thermodynamic properties of the materials. Thanks to the sharp-interface treatment employed in this paper, different EOS can be used in the bulk region of each phase so that any unphysical mixing is avoided. The specific thermodynamic description is given in the following section.

Two-phase thermodynamics
With an evolving interface Γ(t), the flow domain Ω is divided into two sub- domains Ω liq and Ω vap , which represent the regions occupied by the liquid phase and the vapor phase, respectively. Assuming that the fluid in each region is in a local thermodynamic equilibrium, its state can be determined by any two independent thermodynamic variables. In this paper, the fluids water and n-dodecane are considered, which are described by the stiffenedgas EOS [44,12] and the Helmholtz-energy-based EOS [45], respectively.
For the stiffened-gas EOS, the relations between different thermodynamic variables are given by where T is the temperature, degenerate to the ideal-gas EOS [44].
In the Helmholtz-energy-based EOS, a non-analytical formulation of the residual Helmholtz energy is employed, whose parameters are obtained by fitting the underlying EOS to the experimental data. Although it is very inefficient, by using this kind of EOS, the prediction error for fluid properties is usually less than %1 [46]. Moreover, the terms used in this kind of EOS varies for different materials so that we cannot simulate different fluids by simply changing the parameters. Here n-dodecane is considered to validate the present method for real-fluid simulations. The specific formulation is not given here for simplicity and the interested readers are referred to Ref. [45].
As phase change is considered, only the fluid states below the critical point [40] are considered such that the liquid phase and the vapor phase coexist. We also require that the so-called saturation curves exist, which can be characterized by where G is the Gibbs free energy, σ the surface tension coefficient, and κ the interface curvature. As in Ref. [41], we restrict ourselves to admissible EOS regions, where, with S being the entropy, the fundamental derivative for sound speed ∂p ∂ρ S is positive. By doing so, the split and composite waves [8,40] are excluded from the two phase Riemann problem considered in Section. 3, whose treatment is beyond the scope of the present study.

Conservative discretization
The main principle of the conservative sharp-interface method [42] is to modify the spatial discretization in the Cartesian cells cut by the interface, i.e., the so-called cut cells. Take a two-dimensional Cartesian grid with spacing ∆x and ∆y as an example. By applying the first-order forward Euler method to Eq. (1) leads to the following discretized equations, where ∆t is the time step and α i,j U i,j is the conservative variables in the cut cell, with α i,j being the volume fraction of the corresponding phase. As illustrated in Fig. 1, A and ∆Γ represent the cell-face apertures and the interface segment of the cut cell, respectively. The reconstructed cell-face numerical inviscid fluxF is obtained by a high-order shock-capturing scheme such as the WENO [47] scheme. The interfacial fluxX, which is crucial for modelling phase change, is obtained by solving a two-phase Riemann problem along the normal direction n (see green arrow in Fig. 1). Let W l i,j (ρ, V, p, e) and W v i,j denote the primitive variables of the liquid phase and the vapor phase in the cut cell (i, j), respectively, where V = u · n is the normal velocity. Then the two-phase Riemann problem can be formulated as in which the superscript (g) indicates the ghost states obtained by using the extending algorithm [48].

Interface description
The interface is tracked by a level set function φ(x, t), which represents the signed distance from the interface to each cell center. The time evolution of φ(x, t) [25] is governed by ∂φ ∂t +û · ∇φ = 0, whereû represents the interface velocity obtained by solving Eq. (6). The normal direction and the curvature can be calculated by respectively.

Approximate Riemann solution with phase change
3.  Figure 2: The wave structures for the general Riemann problem with evaporation and condensation.

Condensation
As mentioned above, a one-dimensional two-phase Riemann problem is solved along the normal direction to obtain the interfacial flux. When phase change is considered, an additional wave is introduced in the Riemann solution fan to model the phase interface. Assuming a left liquid state and a right vapor state, we plot the Riemann wave patterns for evaporation and condensation in Fig. 2 with a (x, t)-diagram. For the other three waves, standard jump conditions [49] can be imposed directly. To ensure thermodynamic consistency [40], the jump conditions at the phase interface are given where [a] = a * vap − a * liq is the jump operator, Q lat the latent heat set as 2240000 J/kg for water and 249410 J/kg for n-dodecane, σ the surface tension coefficient, S p the velocity of the phase interface and j the mass flux with Here the superscript * denotes the states in the intermediate region, the socalled star region, and the subscripts liq and vap are chosen according to the direction of phase change, which is formulated as As in Refs. [40,41], the phase interface is assumed to be subsonic and non- in which c is the sound speed. With the phase interface, the resulting system is underdetermined and an additional phase change model is needed to decide the mass flux j.

The four-wave approximate Riemann solver
Based on the work of Fechter et al. [41], a novel approximate Riemann solver is developed by reconstructing a contact wave in the Riemann wave fan. As illustrated in Fig. 3, the solution consists of five regions with constant states, which are separated by four waves including two shock waves, a phase interface, and a contact wave. In the approximate solution, the wave speeds of the shock waves are estimated a priori. In this paper, we use the simple estimates suggested by Davis [50], which read The jump conditions at the shock waves are given by where k = L, R. For the contact wave, the standard Rankine-Hugoniot equations can be simplified to in which S c denote the speed of the contact wave, k = L or k = R represents the condensation case or the evaporation case, respectively. Up to now, with the jump conditions of the four waves Eqs. (9,14,15) and the relation of the mass flux Eq. (10), we obtain a simplified system with thirteen independent equations. However, there are fifteen known variables in total: Note that the pressure in the star region does not agree with the thermodynamic pressure calculated by using the EOS, which is dropped due to the HLL approximation [51].
To solve the current system, two additional equations need to be added.
For the first one, in this paper, we make an assumption that The second equation is a phase change model for calculating the mass flux, Following Refs. [35,36], we employ the Hertz-Knudsen relation (HKR) [52], where p sat is the saturation pressure, R v is the specific gas constant, λ evap and λ cond are the evaporation and condensation coefficients, respectively. In this paper, R v is set to 461.52 J/(kg · K) for water and 48.81 J/(kg · K) for n-dodecane. As in Ref. [53], the saturation pressure of water is calculated by while the saturation pressure of n-dodecane is obtained via the open-source library CoolProp [54]. Special attention needs to be paid to the calculation of the temperature. As mentioned above, in the star region, the density, the pressure, and the internal energy do not agree with the EOS. In consequence, we have three ways to calculate the temperature: , and T * = T * (p * , e * ). The last one is used in the present method to eliminate the approximate error of density introduced by Eq. (16).
By using all equations except the phase change model Eq. (17), the unknown variables W * M , W * R , S p , and S c can be expressed in terms of j, which are given in detail in Appendix A. The one-dimensional reduced system can be expressed by where j is the only unknown variable on the right side. In this paper, the adjacent states of the phase interface W * liq and W * vap are used to calculate the mass flux (see Eq. (10)), while the initial states W L and W R are used in Refs. [35,36]. In fact, as shown in Eq. (20), we can obtain the approximate Riemann solution non-iteratively if j is evaluated with W L and W R . However, as we will show in Section. 4.1, this simple strategy will lead to an inconsistent numerical solution, which is not discussed in the previous researches. With the relations W * liq = W * liq (j) and W * vap = W * vap (j), the target equation Eq. (17) is rewritten as j = f m (j) so that it can be solved with a standard root-finding algorithm. In this paper, we employ the steffense method [55], which is formulated as where the initial value j 0 is estimated by using W L and W R . This iteration procedure will be terminated when |j k+1 − j k | < 10 −6 and 3-4 iterations are typically enough. Nevertheless, since the Euler equation is used and the heat conduction is neglected, the temperature jump at the contact wave can be arbitrary large. As reported in Ref. [40], being the potential source of numerical instability, unphysical large temperature jumps will be predicted at the contact wave in some complex problems. To ensure numerical stability, we add an additional constraint for the temperature difference at the contact wave. When it is larger than 50 K, the iteration will be aborted with the mass flux fixed at j = j 0 .
After obtaining the approximate Riemann solution, the interface velocity is set as S p . Following Ref. [17], we extend the interface velocity to the cells in the vicinity of the interface along the normal direction. In addition, the interfacial fluxes can be computed bŷ where + and − are applied for the vapor phase and the liquid phase, respectively. The conservation property of the present method is confirmed by the relation as follows, Note that the momentum conservation and the energy conservation are not satisfied, and the reason is twofold. First, the pressure jump caused by surface tension leads to the violation of the momentum conservation and the energy conservation. However, as the present method is based on the conservative sharp-interface framework, it maintains the so called zero-order consistency [56], i.e., the total force acting on the drop with constant curvature is zero.
Second, to model the phase interface in a thermodynamic consistent way with the heat conduction neglected [41], the latent heat is introduced in the energy jump condition as a source term. The phase interface will liberate thermal energy in the evaporation process while the opposite situation happens in the condensation process, which are known as the heat of evaporation and condensation, respectively. When phase change and surface tension are not considered, the interfacial fluxes Eq. (22) will recover the original scheme of Hu et al. [42].

Numerical examples
In this section, various numerical examples are presented to validate the present method. The inviscid term is discretized by a fifth-order accurate WENO scheme [47], which is also used to solve the level-set advection equa- and the final time is t = 0.7 ms. As mentioned above, the mass flux can be evaluated through the adjacent states of the phase interface (W * liq and W * vap ) or the initial states (W L and W R ). The results obtained with the two strategies are given in Fig. 4, which are labeled by "Star" and "LR", respectively. Note that there is a common exact solution for the two situations, which is achieved via setting the model coefficients in Eq. (18) , the numerical results agree well with the exact solution. In contrast, the results obtained with j = f m (W L , W R ) will converge to an incorrect solution. This issue is not even noticed in the previous researches, and it will also occur in the condensation case. For simplicity, the evaporation case is chosen to elucidate this issue.   Let S denotes the unknown states in the star region. Then the exact Riemann solution can be expressed as where the superscript 0 represents the initial conditions given in Eq. (24).
Therefore, in the numerical simulation, the initial states of the interfacial Riemann problem are W 0 L and W 0 R at the beginning and then change to W * e L and W * e M within a few time steps. If they are used to evaluate the mass flux, an inconsistent numerical solution will arise due to f m (W 0 . Conversely, the mass flux obtained by using the adjacent states of the phase interface is consistent during the whole simulation since the relation holds. In conclusion, no matter which phase change model is adopted, the mass flux should be calculated through the adjacent states of the phase interface to ensure numerical consistency.

Riemann problem of n-dodecane with condensation
We consider the Riemann problem of n-dodecane with condensation here.
With λ evap = 0.6 and λ cond = 1.0, the initial conditions are given by

Riemann problem of water with condensation
Here we consider the Riemann problem of water with condensation. The The computational domain is a unit square with symmetry boundary conditions. Note that only 1/4 of the droplet is considered here, whose center coincides with the origin. With λ evap = 1.0 and λ cond = 0.9, the twodimensional simulation is carried on a MR grid where the effective resolution where the superscripts 0 and n represent the initial condition and the states at the n-th time step. Fig. 10(a) plots the relative mass of each phase versus simulation time, in which the mass transfer from the liquid phase to the vapor phase can be observed. The conservation of the present method is confirmed by the total relative mass equaling to one during the whole simulation, see

Circular bubble with condensation
The two-dimensional condensation of an n-dodecane bubble is computed here. As in Section  grid points, see Fig. 9. The multi-dimensional implementation of the present method reproduces well the results of the one-dimensional approach assuming radially symmetry [49]. Note that the two-dimensional results are obtained with a coarse resolution of only 256 × 256. The mass transfer from the vapor phase to the liquid phase can be observed from Fig. 10(b), which also confirms that the present method has no conservation error.

Oscillating droplet
As in Refs. [40,41], the oscillating droplet case is considered to validate the present method for problems dominated by surface tension. Initially, with the semi-axes a = 0.6 mm and b = 0.4 mm, an elliptical droplet is  placed at the center of a 5 mm × 5 mm square domain with zero kinetic energy. Due to surface tension, this droplet will oscillate periodically along with the translation between the potential energy and the kinetic energy. The symmetry property of the present case is utilized by computing only 1/4 of the physical domain. We set the surface tension coefficient as σ = 20 Nm −1 , and the initial conditions as A for details. Fig. 11(a)-(c) depict the shapes of the droplet at different time instants. Qualitatively, no obvious differences are observed between the interfaces obtained with and phase change, as the phase change effects are small according to the initial conditions. For a quantitative analysis, the lengths of the semi-axes are plotted over time in Fig. 11(d). With the phase change effects, the oscillation amplitude is changed by the additional forces acting on the surface while the oscillation frequency remains the same.
According to the Rayleigh formula extended to two-phase flows [58,59], the analytical oscillation period can be calculated by where o is the oscillation mode and R is the radius of the liquid drop at the equilibrium state. In this case, we have o = 2 and R = 0.4899 mm so that the theoretical oscillation period is T = 0.204 ms. As shown in Fig. 11(d),

Two-dimensional shock and water droplet interaction
Here we study the interaction between a shock wave of Mach number 1. and the initial conditions are

Three-dimensional shock and water droplet interaction
This case is a three-dimensional extension of the case in Section. 4.4.2.
All the initial conditions are the same except the interface geometry. Here we consider a three-dimensional spherical water droplet, which is initialized

Concluding remarks
In this paper, we have developed a fully conservative sharp-interface method for compressible multiphase flows with phase change. In the previous researches [30,35,36,40,41], the coupling between the liquid phase and the vapor phase is usually implemented through the GFM, which is inherently non-conservative. Based on the work of Hu et al. [42], the present method hire the interfacial fluxes to ensure strict conservation, which are obtained by solving a general two-phase Riemann problem with phase change.
In compare with the methods in Refs. [43,44], the present method makes no assumption for the time-scale feature of the phase change process and can guarantee the thermodynamic consistency. A novel approximate Riemann solver employing a four-wave model is developed to solve the Riemann problem efficiently. The original eight-dimensional root-finding procedure in the exact Riemann solver [40] is simplified as only an iteration of the mass flux. Unlike the approximate Riemann solver of Fechter et al. [41], where the energy coupling in the vapor phase is dropped, the present approximate solver impose the jump conditions of all waves strictly. Different strategies for choosing the liquid and the vapor states used in phase change model are discussed in this paper. We have shown that, no matter which model is con- the two-phase interface will be considered in the future work.
For this limit case, the interfacial fluxes are calculated by using W * M and W * R .