Kinetic simulations of compressible non-ideal fluids: from supercritical to phase-change and exotic behavior

In this study, we present a thorough investigation of a compressible kinetic model for non-ideal fluids [DOI:10.1103/PhysRevE.102.020103]. The model imposes the local thermodynamic pressure through appropriate rescaling of the particle's velocities, which accounts for both long- and short-range effects and hence full thermodynamic consistency. The model is fully Galilean invariant and treats mass, momentum, and energy as local conservation laws. After detailed analysis and derivation of the hydrodynamic limit, the model's accuracy and robustness is assessed for various benchmark simulations ranging from Joule-Thompson effect, phase-change and high-speed flows. We show that our model can operate in the entire phase diagram, including super- as well as sub-critical regimes and inherently captures phase-change phenomena.


Introduction
The lattice Boltzmann method (LBM) is a kinetic-based approach for the simulation of hydrodynamic phenomena with applications ranging from turbulence [1,2] to micro [3,4] and multiphase flows [5][6][7][8]. The fully discretized kinetic equations evolve particle distribution functions (populations) f i (x, t), associated with a set of discrete velocities c i , according to a simple stream-and-collide algorithm, which recovers the Navier-Stokes equations in the hydrodynamic limit [9].
Despite significant success of the LBM in isothermal incompressible flow applications [10][11][12][13], the simulation of compressible flows with the LBM is still an active field of research. Furthermore, the majority of studies concern ideal gases, which is inherent to the LBM. However, the ideal-gas assumption is no longer valid in many scientific and engineering applications and real-gas effects have to be taken into account. This includes phenomena such as rarefaction shock waves [14][15][16][17], the acoustic emission instability [18,19], inversion line (change of sign of the Joule-Thomson coefficient), phase transition, surface tension and super-critical flows. Hence, developing a compressible scheme with a non-ideal equation of state (EOS) remains an uncharted territory within the LBM and is subject to this paper.
For incompressible flows, two main approaches have been proposed so far to model non-ideal gas flows within the LBM context: pressure-based methods [20,21] and forcing methods [22][23][24][25]. In pressure-based methods, the equilibrium populations are altered directly and constructed such that the full non-ideal pressure tensor, including the non-ideal EOS and the Korteweg stress, are recovered. However, as shown in [26] these methods lack thermodynamic consistency. On the other hand, forcing methods account for the deviation from the ideal-gas pressure by an appropriate (non-local) force term, which is introduced in the kinetic equations. Promising results have been obtained for various of applications, ranging from droplet collisions at relatively large density ratios [27,28] to droplet impact on textured [29] and flexible surfaces [30]. Furthermore, these methods have been extended to thermal multiphase flows including phase-change [31,32]. Note however that these models are still assuming incompressible flow and thus do not recover the correct energy equation of the two-phase system. The temperature field is then typically modeled as a passive scalar field under various assumptions, which leads to a tailored form of the energy equation with additional correction terms to account for the latent heat [33].
To mitigate these shortcomings, we recently proposed a novel method for non-ideal compressible fluid dynamics [34] based on adaptive discrete velocities in accordance to local flow conditions. In contrast to the aforementioned schemes, the model features full Galilean-invariance and is thermodynamicaly consistent. As a consequence, the full energy-equation of a non-ideal fluid is recovered, which means that no additional phase-change model is required. This enables us to capture a large range of flow regimes, which we aim to explore in this paper. In particular, we will assess the model's performance and range of validity for benchmarks ranging from super-critical flows, throttling, phase change to shock-stability.
The paper is structured as follows: Section 2 provides an in-depth analysis of the model, starting by a presentation of the discrete kinetic equations in Sec. 2.1. Finally, the Chapman-Enskog analysis and the derivation of the hydrodynamic limit of the model is carried out in Sec. 2.2. Numerical benchmarks including the simulation of the Joule-Thomson effect, phase-change and high-speed flows are presented in Sec. 3. Finally, conclusions are drawn in Sec. 4.

Kinetic equations
Our thermokinetic model of non-ideal fluids [34] is based on the so-called Particles-on-Demand (PonD) method [35], which constructs the particle's velocities relative to the reference frame (gauge) λ = {T, u}, where T is the local temperature and u is the local velocity. While the former leads to thermodynamic consistency, the latter guarantees Galilean invariance. In addition, in the local reference frame, the local equilibrium becomes exact and solely dependent on the density. This is in contrast to classical LBM where one typically resorts to a truncated polynomial. The populations can be transformed between different reference frames by requiring the moments to be independent of the reference frame [35]. In [34], we generalized this concept to encompass the thermodynamics of non-ideal fluids by defining the new set of discrete velocities as Define reference frame based on the values from the previous time-step: Define discrete velocities: Semi-Lagrangian advection: Compute filed variables: Compute new refernce frame: where e = e(ρ, T) is the internal energy of a non-ideal fluid. Subsequently, the pressure can be evaluated using the EOS of choice using the updated values for density and temperature. Finally, we can define the corrector gauge as

Convergence of reference frame? Collision
This predictor-corrector step is iterated until the convergence of the gauge is achieved, where the limit gauge is the co-moving reference frame. Once the co-moving reference frame is determined, the advections (2) and (3) are completed. Finally, the collision step in the co-moving reference frame follows. For the f populations the local equilibrium takes the exact form f where W i are lattice weights, which are known for any velocity set. The equilibrium for the g population is derived using Grad's approximation [38][39][40] for the new discrete velocities (1). Thus, the equilibrium populations are constructed from the moments where the explicit relations for the equilibrium moments are given by Here, we shall define a second-order polynomial based on the general discrete velocities v iα = √ θc iα + u α , where it can easily be observed that the zeroth-order moment is already recovered. To satisfy the higher order moments, we must solve for M α and M αβ . Substituting (18) in (14) gives Considering the explicit relations of the moments given by Eqs. (15)- (17) and taking into account the scaling θ = p/ρT L one obtains Finally, upon substitution in Eq. (18), g eq i is obtained as, where D is the space dimension. As mentioned earlier, the source term G i in Eq. (5) is responsible for correction of the energy equation. Here, we only attempt to derive the formulation as we will discuss the details later. To derive an expression for the correction term in the co-moving RF, we follow the same method employed to express g eq i . The pertinent moments of the correction term are, Using relations (19) and (20), one can write, which leads to the following solution, Eventually, we can write the polynomial form of the correction term as, where M 0 is the correction term in the energy equation.

Energy conservation in a non-ideal gas
In what follows, we briefly comment on the necessity of using the two-population approach as a result of modeling a generic, non-ideal gas. The recovery of the correct energy equation from the kinetic equations is crucial for modeling non-ideal fluids and in particular, two-phase systems with large variations of density and temperature as encountered in supersonic flows with shocks. To that end, two-population approaches have long been used in the literature to split the conservation laws, where one population accounts for mass and momentum, whereas a separate population carries the total or internal energy only (see, e.g., [41] in the ideal-gas framework or [32] for multiphase flows). We remind that the internal energy of a non-ideal fluid is now a function of both density as well as the temperature, where C v = (∂e/∂T) v is the specific heat at constant volume and v = 1/ρ is the specific volume. One can evaluate that the the second term on the RHS of Eq. (28) vanishes for an ideal-gas. For the sake of presentation, we shall at first neglect the interface energy and only consider E = u 2 /2 + e, where E is the total energy, e = e(s, v) is the local internal energy per unit of mass, s is the entropy and temperature is defined by T = (∂e/∂s) v . The pressure tensor in our approach takes the form P αβ = ρu α u β + pδ αβ , where p is the thermodynamic pressure corresponding to the chosen EOS. Accordingly, the trace of this tensor represents the total energy 2ρE = ρu 2 + 2ρe with In the case of an ideal-gas, this reduces to e = DT, which, as expected, is only a function of temperature. Thus, equation (29) is fully valid for ideal gases with a fixed specific heat C v = D/2. However, the internal energy of a real gas is a function of both density and temperature. In this paper, we use the classical van der Waals (vdW) EOS p = ρRT/(1 − bρ) − aρ 2 to model non-ideal behavior of real-gases but others can be used analogously. The constants are set to a = 2/49, b = 2/21 and R = 1, where a is the long-range attraction parameter, b represents the excluded-volume effect and R is the specific gas constant.
Considering the vdW EOS, we have which suggests that e vdw = F(T) − aρ, where F(T) is an arbitrary function of temperature. In other words, the internal energy of a vdW fluid is the sum of a density-dependent function and temperature-dependent function. The same analogy can be applied to all EOS's with a linear dependence on the temperature and density-dependent tail part; such as Carnahan-Starling EOS. However, evaluating the internal energy of a vdW fluid using equation (29) (for D = 2) yields, which in inconsistent with (31) since the first terms depends on both, the density and the temperature. Hence, the pressure tensor alone cannot represent the total energy in general. This is only true for ideal gas. Consequently a second set of population is needed to represent the energy.

Correction of the energy equation
We must comment that the expression of heat flux recovered from the Chapman-Enskog analysis (see Section 2.2) without the correction term in the g population is found as q CE α = −µ∂ α h where µ is the shear viscosity and h = e + p/ρ is the specific enthalpy. At the limit of an ideal-gas, this is equivalent to the Fourier law q ig α = −k ig ∂ α T where k ig = µC ig p and hence the Prandtl number is fixed to Pr = µC ig p /k ig = 1 due to the single relaxation time BGK collision model. However, considering the enthalpy of a real-gas as a function of pressure and temperature, we have, While one could eliminate the pressure part of the enthalpy by the correction term and only retain the temperature dependent part, it must be noted that the thermal expansion coefficient at the critical point diverges, β → ∞ and so does the specific heat C p → ∞. Hence to recover the well-known Fourier law, the post-collision of the g population is augmented by the correction term and k is the conductivity to be set independently.

Surface tension
In order to describe two-phase flows in the sub-critical part of the phase-diagram, the collision step for the f -populations (4) is augmented with a source (forcing) term S λ i , where δu = F/ρδt is the change of the local flow velocity due to the force F α = ∂ β K αβ , where is the Kortwewg stress [8], ∆ = ∇ 2 is the laplacian operator and κ is the surface tension coefficient. The first term on the R.H.S of equation (33) denotes the transformation of equilibrium populations residing at the reference frame "u + δu" to the reference frame "u" which could be realized as the Exact Difference Method (EDM) [42] adapted to the comoving RF. Having included the source term S i , the actual fluid velocity is now shifted In the presence of the interface, the local equilibrium (22) is extended to account for the forcing F. In order to do that, the same analogy used in the absence of the force term is employed here with the only difference that the velocity terms in the pertinent moments (15)- (17) are replaced by the modified velocities, With the new settings considered, the solution to relations (19) and (20) gives, which can be written in a compact form using the expressionû whereĤ =Ê + p/ρ = h +û 2 /2. Finally, the extended equilibrium takes the form, where ρÊ = ρE +û 2 /2 is based on the actual velocity of the flow. Note that in the absence of the force, the equilibrium (40) simplifies to Eq. (22). Consequently, the corresponding work of the added force is taken into account in the energy equation by modifying the correction term (27) Finally, with the above modifications, the hydrodynamic equations for a two-phase system are recovered in its correct form. The evolution equations (69)-(71) together with the stress tensor (72) remain intact however, all the velocity terms u are replaced by the actual velocityû. Furthermore, the standard form of the total-energy conservation for a two-phase system [43] is recovered, where ρÊ = ρÊ + κ 2 |∇ρ| 2 accounts for the excess energy of the interface, as well. It is essential to mention that the van der Waals formulation of a real-gas, features negative values of pressure for a range of temperatures in the subcritical region (T r < 0.84375) where a constant base-pressure must be added in order to have a meaningful evaluation of discrete velocities (1). To resolve this issue, we redefine the pressure as p = p vdw +p which is positive. This will contribute to the internal energy according to relation (28), where the pressure-dependent part is reevaluated as, Finally, the internal energy of a vdW fluid with base-pressurep and constant specific heat C v is derived, However, the enthalpy of such a fluid remains intact since, where the effect of the base-pressure is cancelled out in evaluation of the enthalpy.

Excluding the forcing term
Here we aim at recovering the macroscopic Navier-Stokes equations from the dynamics of kinetic equations (2)(3)(4)(5). To this end, the pertinent equilibrium moments of f and g populations are required, which are computed as follows: where [uδ] αβγ = u α δ βγ + u β δ αγ + u γ δ αβ and H is the total enthalpy. First, we introduce the following expansions: Applying the Taylor expansion up to second order and separating the orders of results in: i }.
The local conservation of density, momentum and energy imply Applying conditions (57) and (58) on equation (55), we derive the following first order equations, where D (1) α is the first order total-derivative. Subsequently, we can derive a similar equation for pressure considering that p = p(ρ, T). This yields where ς = ∂p ∂ρ s is the speed of sound given by, The second order relations are obtained by applying the conditions (57) and (58) to equation (56), Equations (59) and (64) Finally, summing up the contributions of density, momentum and temperature at the and 2 orders and taking into account the correction to the energy equation (27), we get the hydrodynamic limit of the model, which reads where D t = ∂ t + u α ∂ α is the material derivative, q α = −k∂ α T is the heat flux and the nonequilibrium stress tensor reads, The shear and bulk viscosity are, respectively and ς = (∂p/∂ρ) s is the speed of sound. As expected, the bulk viscosity vanishes in the limit of ideal monatomic gas. Furthermore, one can observe that only the excluded volume part of the pressure T(∂p/∂T) v contributes to the temperature equation (71), as expected.

Including the forcing term
As mentioned before, the force terms in the kinetic equations represent the interface dynamics. First, we recast the post-collision state of the f population in the following form [44], Here we expand the forcing termŜ i in addition to the expansions (50,52,53). Similarly, we get the following relations at the orders of 0 , 1 , 2 respectively, It is important here to assess the solvability conditions imposed by the local conservations. Considering the moment-invariant property of the transfer matrix between the two gauges λ = {p/ρ, u} andλ = {p/ρ,û}, one can easily compute, This implies that, According to the definition of S i , the following moments can be computed: Similarly, the first order equations of density and momentum are derived by applying the solvability conditions (83) and (84) on equations (79) and (80), whereD (1) α . At this point it is necessary to mention that since there is a force added to the momentum equation (in this case the divergence of the Korteweg stress), it should also be considered in the energy equation as well. Hence as mentioned before, M 0 in Eq. (27) is modified to, The equilibrium moments are modified as, whereÊ = e +û 2 /2 andĤ =Ê + p/ρ. With the changes mentioned so far, the first-order equation of temperature is derived asD αûα .
Finally, in a similar manner as the case without the force, the macroscopic equations are recovered by collecting the equations of density, momentum and temperature at each order, It should be noted that the error terms associated with the forcing are not shown here. For instance, as reported in the literature [44][45][46] one can show that the error term in the momentum equation appears as ∇ · (δt 2 FF/4ρ). The total energy of the fluid is formulated byÊ = e (T, v) +û 2 /2 + E λ where E λ = κ|∇ρ| 2 /2 is the non-local part corresponding to the excess energy of the interface. The evolution equation for the specific internal energy e(T, v) can be obtained by considering equations (28,95,97) ρD t e = −p∂ αûα −τ αβ ∂ αûβ − ∂ α q α .
From the momentum equation (96) we get, and the evolution of the excess energy can be computed using the continuity equation, by summing up the contribution of each three part we get the full conservation equation for the total energy, ∂ t ρÊ + ∂ α ρÊû α + pû α +τ αβûβ + K αβûβ + κρ∂ βûβ ∂ α ρ + q α = 0.

Results and discussion
In this section, we show validity of our model in a broad range of challenging problems, which are chosen to probe the correct thermodynamics as well as Galilean invariance: • As a first test of basic thermodynamic consistency for non-ideal fluids, we simulate the inversion line of a vdW fluid, which is one of the classic thermodynamical concepts of non-ideal fluids. To capture this phenomenon it is crucial that the model recovers the correct energy equation and can operate in a wide range of pressures and temperatures in the super-critical part of the phase diagram. • Phase-change is the next fundamental process that is tested with our model. It is important to remind that since the full energy equation is recovered by our kinetic equations, phase-chase emerges naturally in the proposed scheme and no additional ad-hoc phase-change model is required. In addition, we probe fast dynamics with temperatures near the critical point, where phase-change happens in short time-scales. • As a final test case we probe both thermodynamic consistency as well as Galilean invariance in supersonic flows. In particular, we study the behavior of a perturbed shock-front in both an ideal-gas as well as a vdW fluids at Mach number Ma=3. In agreement with theory our model shows to capture all regimes, including the exotic behaviors of a real fluid.

Inversion line
When a fluid passes through a throttling device, the value of the enthalpy remains constant in the absence of work and heat. During this so-called throttling process, the pressure of the fluid drops and the behavior of the temperature is characterized by the well-known Joule-Thomson coefficient µ = (∂T/∂P) h . Depending on the sign and value of this coefficient, the temperature may increase, decrease or remain the constant through the process. It is trivial to talk about the ideal-gas assumption here since the coefficient is zero by definition and the temperature of an ideal-gas will thus not change. On the other hand, for real gases, we need to distinguish between three different regions in the T − P diagram, corresponding to the different signs of the Joule-Thomson coefficient. Hence, let us start by defining the inversion line as the locus of points where µ = 0. Hence, crossing the inversion line will lead to a change of sign of µ. For the vdW EOS, one can derive the expression for the inversion line as where the subscript "r" indicates that the quantities are reduced with respect to their values at the critical point. The critical values of pressure, temperature and density for a vdW fluid are P cr = a/27b 2 , T cr = 8a/27Rb and ρ cr = 1/3b, respectively. In addition to the reduced variables, it is useful to define the non-dimensional enthalpy as,ĥ where δ = R/C v is a constant. A closer assessment of (103) reveals that the point with the maximum pressure on the inversion line has the following quantities P r , ρ r , T r ,ĥ = (9.0, 1.0, 3.0, 11.25).
To test that our model captures these phenomena also numerically, we measure the value of Joule-Thomson coefficient at different points in the T − P diagram. This will then allow us to infer the inversion line, which can be compared to its analytical solution. We do this in two steps: in a first simulation the flow is subjected to a positive acceleration under fixed density; hence the pressure drops and the quantity (∂P/∂T) ρ is measured. In a second simulation, the isothermal speed of sound (∂P/∂ρ) T is computed by introducing a very small perturbation in the pressure field and measuring the velocity of the subsequent shock front. The Joule-Thomson coefficient may now be computed by, where C p is the specific heat at constant pressure. Finally, we use a simple Euler scheme to construct the isenthalpic lines with, ∆T ≈ µ∆P.
The simulations are conducted for three different enthalpies,ĥ = 5,ĥ = 11.25 andĥ = 15. Figure 2 shows the measured values of the dimensionless Joule-Thomson coefficient at different reduced pressures up to the far supercritical value P r = 15. The comparison between the van der Waals theory and the simulation is excellent and thus validates our scheme.

Phase change: one dimensional Stefan problem
In this section, we validate our model for phase-change problems, starting from the one dimensional Stefan problem, where a liquid-vapor system is subjected to a heated wall on the vapor side. The heat transfer from the wall leads to evaporation of liquid and the interface is moving away from the wall. The analytical solution for the liquid-vapor interface location with time is given by where α v is the diffusivity of the vapor and β is the solution to [47] where St = C pv ∆T/h f g is the Stefan number, C pv is the specific heat capacity of the vapor phase, ∆T is the temperature difference between the wall and the saturation temperature and h f g denotes the latent heat of evaporation. Simulations have been carried out for three different Stefan numbers at fixed diffusivity. The choice of the Stefan number is directly related to the velocity of the interface: and hence the Mach number of the flow. Note that since our model is not restricted to low-speed flows, we can accurately capture a wide range of Stefan numbers. Figure ( given temperature difference ∆T as we approach the critical point due to vanishing of h f g and diverging of C p at the critical point.

Phase change: nucleate boiling
Due to its importance in engineering and real life applications, boiling and its various regimes have always been the focus of many studies, both numerically and theoretically [43,48,49]. To validate our model, we choose a two dimensional setup where a liquid is in direct contact with a wall with high temperature at the middle of the wall. The schematic of the setup is shown in Fig. 4a. The non-uniform temperature of the wall triggers a two dimensional flow and the nucleus starts to appear and rise under the gravitational field. The nucleus continues to rise and grow until the necking is achieved where the bubble is detaching from the nucleus. Once the first bubble is detached from the nucleus and released into the liquid, the nucleus continues to grow and releases a second bubble. This is a periodic process of bubble release, which is a function of surface tension, density ratio and the gravity. An empirical correlation for the bubble release frequency was found experimentally by Zuber [50] and reads: where d is the departure diameter and itself proportional to g −0.5 [51]. Hence, the bubble release period is proportional to g −0.75 . We consider domain of 121 × 601 points with time step δt = 0.3, conductivity k = 0.6, specific heat C v = 3, viscosity ν = 0.005, surface tension coefficient κ = 0.0234 and gravity g = 0.0001 in lattice units. The Jacob number is defined as where C pl is the specific heat of the liquid phase. The wall temperature is set to T w = 1.5T sat , where T sat is the saturation temperature and the initial temperature of the liquid is T sat = 0.9T cr , which fixes the latent heat of evaporation. This choice of parameters leads to the Jacob number Ja = 2.21. Figure 4 b illustrates a sequence of the bubble interface from the early times of the first nucleus development until the first bubble is released into the liquid. The bubble release period was measure for different values of gravity and the results are presented in in Fig. (5). The comparison shows that our numerical simulations agree well with the empirical correlation.

Phase change: Single-mode film boiling
As a final phase-change validation, we conduct simulations of film boiling, where a heated horizontal surface is covered by a thin layer of vapor. The liquid rests on top of the vapor and both phases are initially saturated. Phase-change then takes place at the liquid-vapor interface, where the heat is transported from the hot wall with temperature T w , which is set to be above its saturation temperature T sat . The governing non-dimensional numbers are the Jacob, Prandtl and the Grashof number [52], defined for the vapor phase. l s is the non-dimensional capillary length defined as, and t * = t/ l s /g is the dimensionless time. The well-known experimental correlation as proposed in Klimenko [53] has the following form in the laminar flow regime (Gr ≤ 4.03 × 10 5 ), where Nu is the Nusselt number. In our simulation, we consider domain of 129 × 257 points with δt = 0.3, k = 0.6, C v = 3, ν = 0.005, κ = 0.0234 and g = 0.0001 in lattice units. The non-dimensional numbers are Ja = 0.064, Pr = 0.094 and Gr = 2482.58. Based on the value of the Jacob number, the Nusselt number as computed from the correlation (115) amounts to Nu k = 2.6085861. Initially, the liquid-vapor interface is perturbed with the function where W is the width of the domain. The space-averaged Nusselt number is computed throughout the simulation using where the gradient of the temperature is computed at the wall using finite differences. The evolution of the liquid-vapor interface is shown at three different times in Fig. 6. The first bubble is released at t * ≈ 15, which is then followed by a periodic release of bubbles. The space-averaged Nusselt number is computed during the simulations until the first bubble is released. The results are presented in Fig. 7, where we have compared the time-averaged Nusselt number with the correlation (115). We can confirm the time-averaged Nusselt number is in very close agreement with the experimental correlation while according to [53], there is an acceptable ±25% margin on the correlation.

On the stability of the shock waves
The stability of planar shock waves subject to small perturbations have long been investigated since the pioneering work of D'yakov [54] and later modifications of Kontorovich [55], which were the first attempts to study the conditions under which a planar shock with corrugations on its surface would become unstable. The key parameter in the analysis of shock instabilities is the so called "D'yakov parameter" [18], defined as evaluated at the post-shock (downstream) state, where j 2 = (p 1 − p 0 )/(v 0 − v 1 ) is the square of the mass flux across the shock front, v is the specific volume, p is the thermodynamic pressure and the subscript H denotes that the derivative is taken along the "Hugoniot" curve. Furthermore, the subscripts "0" and "1" refer to the upstream (pre-shock) and downstream (post-shock) states, respectively. It has been shown that the necessary condition for stability of a shock wave is [19,56]: where 0 < M 1 < 1 is the downstream Mach number, which is measured in the reference frame that is moving with the shock. Under this condition, linear perturbations imposed on the shock front will asymptotically decay in time as t −3/2 [56]. According to the theory, a planar shock wave is unconditionally stable when propagating through an ideal-gas medium [57]. This can be easily evaluated, where h for an ideal gas EOS yields h D,ig = −1/M 2 0 , which always falls within the stability range (120). On the other hand, for non-ideal fluids, these stability conditions (120) can be violated, which leads to amplification of the perturbations until the structure of the flow filed is altered [56]. It has been shown that the violation of the upper limit of the stability condition (120) corresponds to the splitting of the shock-front into two counter-propagating waves [58], while the violation of the lower limit is associated with the splitting of the shock front into two waves, travelling in the same direction [59]. Extensive theoretical investigations have been carried out to study the dynamics of the isolated planar shock waves propagating in an inviscid fluid medium. Namely, Bates [19,60] has derived analytical expressions for the amplitude of the ripples on the shock front, for initial sinusoidal perturbations.
According to Ref. [60], two families of solutions emerge depending on the sign of the following non-dimensional parameter: where are non-dimensional parameters as a function of downstream conditions and η = ρ 1 /ρ 0 is the compression ratio through the shock. For the case Λ > 0, the solution is given by: where, δx is the amplitude of the ripple, τ = Ukt/η is the non-dimensional time, U is the speed of the shock front in the laboratory reference of frame, k is the wave number of the initial ripple, J 1 is the first-kind Bessel function and the parameters a and b are defined as: and are real numbers. It is interesting to mention that the ideal-gas EOS belongs to this class of solutions. Using the asymptotic approximation J 1 (x) ∼ √ 2(πx) −1/2 cos (x − 3π/4) as x → ∞ and considering Eq. (123), one can confirm that the amplitude of the ripple in a fluid medium with Λ > 0 (such as an ideal gas) will decay in time with the negative power law τ −3/2 in the long-time limit [60]. However, the situation can be different for fluids with an EoS that accept a negative Λ. Finally, for the case Λ < 0, the solution to the initial value problem is [60]: where, σ = −ib is a real number. The presence of the exponential function implies a stronger damping compared to Eq. (123). However, the long time asymptotic is still a function of τ −3/2 in both cases [60].
These theoretical considerations give us the opportunity to test and validate our numerical model also in the high-speed regime for the exotic shock-wave behavior of non-ideal gases. Our simulations consist of a long channel with periodic boundary conditions in the vertical direction. In all cases, the conductivity is set to zero and the viscosity is chosen to take the lowest possible value as long as the simulations are stable. In order to capture the shocks and avoid oscillations at the shock front, a third-order WENO scheme based on a 4-point stencil has been used in the reconstruction process instead of the third-order In all cases, the shock front is initially perturbed with a single-mode sinusoidal function. The ratio of the amplitude to the wave length of the perturbation is 10%. The first two cases fall into the category of stable shocks, where the perturbations on the shock front are expected to decay in time. However, the last case is an example of shock-splitting, which will be discussed below. We now limit our focus on the first two cases. At time t = 0, the shock starts to propagate while it oscillates as it moves further towards the low-pressure side. We then measure the oscillation amplitude and compare it to the analytical expressions. The initial shape of the shock front and its evolution in time is shown in Fig. 8 for the first case, where the oscillation of the front and the damping effect is apparent. Figure 9 shows the bird-eye view of this simulation.
According to the sign of Λ, the magnitude of the ripple for case (1) and case (2) were compared with analytical solutions (123) and (126), respectively. The results are presented in Fig. (10) and are in good agreement with the theory. It is apparent that our model captured the two distinct damping effects accurately, with more pronounced damping for the non-ideal fluid, as expected. We now consider the exotic case (3). As mentioned earlier, non-ideal fluids can show exotic behaviors under certain conditions. Regarding case (3), due to its large specific heat value, the Hugoniot curve passes through regions, where the relation h D < −1 is satisfied. As argued in [14], this, together with the fact that the Hugoniot curve has more than two intersection points with the Rayleigh line (see Fig. 11), can cause the shock front to split into two traveling waves. Figure 12 presents our simulation for this case, where it is visible that the initial perturbation on the shock front has become unstable leading to splitting of the shock. The resulting waves, travel in the same direction with different speeds, as expected. This validates our model also for high-speed flows and shock-waves in real-gas media.

Conclusion
In this paper, we have presented a thorough study of our recently proposed model for compressible non-ideal flows. The model features full Galilean-invariance and the full energy equation is recovered for a non-ideal fluid, accounting also for two-phase systems and the presence of interfaces. It has been shown that the model is able to handle flows which are far into supercritical states. The effect of the inversion line on the T − P diagram was correctly captured for a vdW fluid in a wide range of reduced-pressure; from p r = 3 to p r = 15. In addition, owing to the full energy conservation, the latent heat is already included in the model. This was shown on two different phase-change benchmarks: The one-dimensional Stefan problem and Boiling. As one of the advantages of the model, we were able to choose relatively large Stefan and Jacob numbers, which are scarce in the literature. Finally, the stability of an initially perturbed shock front in ideal gas and vdW fluid mediums with the Mach number Ma ≈ 3 were studied and compared to theoretical predictions. It was observed that the damping effect is much stronger in a vdW fluid as predicted by the inviscid theory. Beside from the fact that all of these simulations were implemented by taking only 9 discrete velocities, the results show that the real-gas effects have been properly captured by It is visible that the Hugoniot curve has more than two intersection points with the Rayleigh line. Also, as the volume decreases, there are regions where h D < −1.
the proposed model.