A positivity-preserving high-order weighted compact nonlinear scheme for compressible gas-liquid flows
Man Long Wong, Jordan B. Angel, Michael F. Barad, Cetin C. Kiris
AA positivity-preserving high-order weighted compactnonlinear scheme for compressible gas-liquid flows
Man Long Wong a,1, ∗ , Jordan B. Angel b,1 , Michael F. Barad b , Cetin C. Kiris b a Science and Technology Corporation, Moffett Field, CA 94035, United States b NASA Ames Research Center, Moffett Field, CA 94035, United States
Abstract
We present a robust, highly accurate, and efficient positivity- and boundedness-preservingdiffuse interface method for the simulations of compressible gas-liquid two-phase flowswith the five-equation model by Allaire et al. [1] using high-order finite difference weightedcompact nonlinear scheme (WCNS) in the explicit form. The equation of states of gasand liquid are given by the ideal gas and stiffened gas laws respectively. Under a mildassumption on the relative magnitude between the ratios of specific heats of the gasand liquid, we can construct limiting procedures for the fifth order incremental-stencilWCNS (WCNS-IS) with the first order Harten–Lax–van Leer contact (HLLC) flux suchthat positive partial densities and squared speed of sound can be ensured in the solutions,together with bounded volume fractions and mass fractions. The limiting procedures arediscretely conservative for all conservative equations in the five-equation model and canalso be easily extended for any other conservative finite difference or finite volume scheme.Numerical tests with liquid water and air are reported to demonstrate the robustness andhigh accuracy of the WCNS-IS with the positivity- and boundedness-preserving limiterseven under extreme conditions.
Keywords: positivity-preserving, boundedness-preserving, weighted essentiallynon-oscillatory (WENO), diffuse interface method, multi-phase, shock-capturing
1. Introduction
In the numerical computations of compressible flows, simulations cannot proceedwhen negative density or squared speed of sound appears because the system of equa-tions becomes ill-posed. This problem is more pronounced for extreme applications suchas those in astrophysics where strong shocks, rarefactions, or blast waves may exist inthe simulations. While one can have successful simulations with some robust (or evenpositivity-preserving) first or second order shock-capturing schemes [11, 27, 3, 28, 15],these schemes are numerically very dissipative and are inefficient for scale-resolving sim-ulations, such as large eddy simulations (LESs) or direct numerical simulations (DNSs). ∗ Corresponding author. These authors contributed equally to the development of positivity- and boundedness-preservinglimiters.
Preprint submitted to Journal of Computational Physics February 11, 2021 a r X i v : . [ phy s i c s . c o m p - ph ] F e b ver the decades, many high-order accurate shock-capturing schemes with more localizednumerical dissipation and higher resolution were developed for scale-resolving simulations[42, 22, 10, 4, 25, 19, 24, 12, 49, 44]. While these high-order shock-capturing schemes havea certain degree of robustness for problems involving shocks and other kinds of discon-tinuities, there is still no guarantee of having successful simulations for severe problemsusing those schemes. Replacing the negative density or squared sound speed with pos-itive ones is not conservative and may trigger other numerical issues such as spuriousoscillations.In recent years, many positivity-preserving limiters [53, 54, 55, 20] have been de-veloped for high-order shock-capturing schemes. These limiters can preserve positiv-ity of density and squared speed of sound for compressible flows. Motivated by thepositivity-preserving technique in Perthame and Shu [33], Zhang and Shu [53] devel-oped positivity-preserving high-order discontinuous Galerkin (DG) schemes for Eulerequations. The method was later extended to Euler equations with source terms [54].Positivity-preserving limiters specifically designed for high-order finite difference schemesdidn’t appear until the works by Zhang and Shu [55] and Hu et al. [20] where extremeproblems could be successfully simulated with the finite difference weighted essentiallynon-oscillatory (WENO) schemes using positivity-preserving limiters. These methods areconservative as the limiters are applied to the fluxes directly in the conservation form.All of the positivity-preserving limiters aforementioned were designed for compressiblesingle-phase flows. The appearance of non-physical states also happens in compressiblemulti-phase simulations. In general it is more likely for non-physical states to appear dueto higher density gradients across material interfaces in the related applications such assupersonic combustion, cavitation erosion, break-up of high-speed liquid jets, water-basedacoustic suppression systems, etc. In addition to having negative density and squaredspeed of sound, the solutions are also considered non-physical if mass fractions or volumefractions are not bounded between zero and one. In order to address the numerical issues,some boundedness-preserving diffuse interface methods [39, 21] in the Eulerian frame-work have been proposed. Shen et al. [39] adopted the maximum-principle-satisfyinglimiter by Zhang and Shu [52] for a space-time conservation element and solution element(CE/SE) scheme. The scheme can ensure the boundedness of volume fractions in thefive-equation model by Allaire et al. [1] for multi-phase flows. Another thermodynamics-consistent boundedness-preserving scheme by Jain et al. [21] for the same flow modelwas developed with the use of interface-regularization terms. Although both methodscan preserve boundedness of volume fractions, non-physical states can still appear in themulti-phase simulations since partial densities or squared sound speed can still becomenegative. A positivity-preserving high-order method by Cheng and Shu [7] was proposedfor multi-phase simulations in the Lagrangian framework. Compared to Eulerian diffuseinterface methods, Lagrangian methods can be more accurate at material interfaces sincethe computational mesh moves with the fluids. On the other hand, diffuse interface meth-ods in the Eulerian framework is more attractive for flows involving large deformationsas the degree of deformations is limited by mesh distortions in Lagrangian methods [35].Motivated by the need for the simulations of water sound suppression systems in rocketlaunch environments that involve interactions between strong shocks and complex air-water interfaces, we propose a high-order positivity-preserving diffuse interface methodin the Eulerian framework targeting gas-liquid two-phase flows with large deformations,where the gas and liquid are described by the ideal gas and stiffened gas equation of2tates respectively. Unlike the previous works [17, 18] that are based on the homoge-neous relaxation model, the flows in this work are described by the five-equation modelby Allaire et al. [1]. Under a necessary but generally valid assumption that the ratioof specific heat of the gas is smaller than that of the liquid, the numerical method canensure physically admissible states with positive partial density of each phase, positivesquared sound speed, and bounded volume fractions and mass fractions.The high-order shock-capturing scheme used in this work is based on the explicit finitedifference formulation of weighted compact nonlinear schemes (WCNSs) [10, 31, 51, 29, 9,30, 49] and the nonlinear weighting technique of the incremental-stencil WENO (WENO-IS) scheme [48]. The use of WCNS as a diffuse interface method for the five-equationmodel [1] has already been demonstrated by Wong and Lele [49] but it is only appliedfor single-phase flows with mixture of ideal gases. The WENO-IS scheme was originallydesigned as a finite volume scheme by Wang et al. [48] for compressible multi-phase flowswith shocks and material interfaces using the same five-equation model. Although therobustness of the scheme was demonstrated in that paper, the finite volume approach iscomputationally more expensive than the finite difference WCNS for multi-dimensionalproblems, while the orders of accuracy are similar [36]. The WCNS with the WENO-ISnonlinear weighting technique, WCNS-IS, presented in this work is more efficient andhas similar robustness in minimizing spurious oscillations in simulations.In this work, we first show the convexity of the physically admissible set of solu-tion states, under the mild assumption on the relative magnitude between the ratios ofspecific heats of the ideal gas and the liquid. We then prove the positivity-preservingand boundedness-preserving properties of the first order Harten–Lax–van Leer contact(HLLC) flux for gas-liquid flows with our choice of the advection velocity of the ma-terial interface using the convexity of the admissible set. Based on the positivity- andboundedness-preserving properties of the first order HLLC flux, we propose a limiter toblend a flux from any Cartesian conservative high-order shock-capturing schemes with theHLLC flux. The flux limiter together with a limiter for WENO interpolation can ensurethe positivity-preserving and boundedness-preserving properties of the overall scheme.The fifth order accurate WCNS-IS formulation presented in this work is also provedmathematically and shown numerically to be high-order accurate in smooth advectionproblems, while robust because of the use of lower order interpolation near discontinuitiessuch as shocks or material interfaces. We have demonstrated that the WCNS-IS schemewith the positivity- and boundedness-preserving limiters, PP-WCNS-IS, can successfullysimulate very intense one-dimensional (1D) and two-dimensional (2D) air-water problemssuch as Mach 10 shock-water column interaction and Mach 100 water jet problems. Theresults also show that smaller errors are produced at shocks and material interfaces, andfine-scale flow features such as vortices are better captured with the high-order schemecompared to the first order HLLC scheme due to more localized numerical dissipation andhigher resolution of the former method. All of the numerical tests highlight the robust-ness of the overall positivity- and boundedness-preserving finite difference scheme, anddemonstrate the method as a highly accurate diffuse interface method for scale-resolvingcompressible gas-liquid simulations. 3 . Governing equations The five-equation single-velocity, single-pressure model proposed by Allaire et al. [1]for compressible two-phase flows is considered in this work. The flow model has thefollowing form: ∂ t ( α ρ ) + ∇ · ( α ρ u ) = 0 , (1) ∂ t ( α ρ ) + ∇ · ( α ρ u ) = 0 , (2) ∂ t ( ρ u ) + ∇ · ( ρ u ⊗ u ) + ∇ p = 0 , (3) ∂ t E + ∇ · [( E + p ) u ] = 0 , (4) ∂ t α + u · ∇ α = 0 , (5)where u is the velocity vector and p is the mixture pressure. u = u and u = ( u v ) T for1D and 2D cases respectively (similar extension for the three-dimensional case). α k and ρ k are respectively the volume fraction and phasic density of phase k , where k = 1 , α k ρ k = ρY k is called partial density of phase k , where ρ = α ρ + α ρ is the mixturedensity and Y k is mass fraction of phase k . E = ρ ( e + | u | /
2) is the mixture total energyper unit volume, where e is the mixture specific internal energy. Also, α + α = 1 and ρe = α ρ e + α ρ e , where e k is the phasic specific internal energy of phase k . Thesystem is closed by the mechanical equilibrium and equation of state of each phase . Thestiffened gas equation of state is chosen in this work because of its popularity for gasesand liquids. The equation of state was first proposed by Harlow and Amsden [16]. Thestiffened gas equation of state of each phase is given by: p k γ k − γ k p ∞ k γ k − ρ k e k , (6)where γ k and p ∞ k are fitting parameters for each of the fluids. γ k is the ratio of specificheats that is greater than one and p ∞ k is non-negative. The stiffened gas equation of stateis reduced to the ideal gas equation of state if p ∞ k = 0. With the mechanical equilibriumassumption, i.e. p = p = p , we can obtain the mixture equation of state by multiplyingequation (6) by α k for each phase and then summing over all the phases: pγ − γ p ∞ γ − ρe, (7)where γ and p ∞ are properties of the mixture. They can be defined by the followingrelations: 1 γ − α γ − α γ − , (8) γ p ∞ γ − α γ p ∞ γ − α γ p ∞ γ − . (9) In this work, each phase consists of one species.
4e define the conservative variable vector as W = ( ρ α ρ α ρ u E α ) T . In thiswork, we also define the set of admissible states as, G = W = α ρ α ρ ρ u Eα (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ α ≤ , α ρ ≥ , α ρ ≥ , ρc > . (10)This requires boundedness of the volume fractions, positivity of the partial densities, andpositivity of the squared speed of sound c . The positive squared speed of sound impliesthat the system of equations remains hyperbolic with real wave speeds. Note that thepositive partial densities also mean that all mass fractions Y k are bounded between zeroand one. The flow model has a mixture speed of sound c given by [1]: ρc = γ ( p + p ∞ )= γ ( γ −
1) ( ρe − p ∞ ) = γ ( γ − (cid:18) E − | ρ u | ρ − p ∞ (cid:19) . (11)Both γ and p ∞ only depend on α in W .It is obvious that α k and α k ρ k are both concave functions of the conserved variables W . Since all γ k are greater than one, γ > ≤ α ≤
1. Therefore, an equivalentphysically admissible set is, G = W = α ρ α ρ ρ u Eα (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ α ≤ , α ρ ≥ , α ρ ≥ , ρe − p ∞ > . (12) Lemma 1. If p ∞ = 0 ( p ∞ = 0 ), γ ≥ γ ( γ ≥ γ ) and ≤ α ≤ , the function ρe − p ∞ is a concave function of the conserved variables W .Proof. The non-zero eigenvalues of the Hessian matrix of the function are: (cid:26) − α ρ + α ρ , − | ρ u | + ( α ρ + α ρ ) ( α ρ + α ρ ) , − ∂ p ∞ ∂α (cid:27) . (13)The first non-zero eigenvalue does not exist for the 1D case. If the last non-zero eigenvalueis non-positive, then the lemma is proved since all eigenvalues are non-positive. Here, weconsider the p ∞ = 0 case and the other case is symmetric. For p ∞ = 0 we have, ∂ p ∞ ∂α = 2 γ p ∞ γ − γ − − γ − (cid:16) α γ − + α γ − + 1 (cid:17) α (cid:16) γ − − γ − (cid:17)(cid:16) α γ − + α γ − + 1 (cid:17) − , (14)and it is clear that if γ ≥ γ , the last eigenvalue is less than or equal to zero. Notethat if γ > γ , the function cannot be concave. In this work, we assume that one of thephases is an ideal gas with p ∞ k = 0. Strictly speaking, α is not a conservative variable. heorem 1. The set G is a convex set.Proof. This is an immediate consequence of Lemma 1 and Jensen’s inequality.The proof of convexity of G relies on the assumption that the ratio of specific heats ofthe liquid ( p ∞ k >
0) is larger than that of the ideal gas ( p ∞ k = 0). However, this is a verymild assumption based on most gas-liquid test problems found in previous literature. Theideal gases considered in the previous works are usually either monatomic gases or air,where the ratios of specific heats are around 1.67 and 1.4 respectively. When the stiffenedgas equation of state is chosen for the liquid in tests, the most popular liquid used iswater, where the ratio of specific heats mostly ranges between 4–7 [34, 43, 8, 48, 32]. Arare but justified choice of γ = 1 .
932 is utilized in [5] but that ratio of specific heats isstill larger than those of monatomic gases or air. Other liquids commonly found in theliterature include ethanol [32] with γ = 2 . γ = 8 .
2. Both liquidsalso have larger specific heat ratios than the ideal gases considered.In the following sections, we make the assumption that ratio of specific heats of liquidis larger than that of ideal gas and formulate our numerical method by taking advantageof the convexity of G such that the solution update is a convex combination of statesalready in G , thus is also in G .
3. First order positivity- and boundedness-preserving scheme with HLLCRiemann solver
The flux given by the HLLC Riemann solver is illustrated in this section. For sim-plicity, a 2D case with domain [ x a , x b ] × [ y a , y b ] is considered with the following equationin compact form: ∂ t W + ∂ x F x ( W ) + ∂ y F y ( W ) + Σ ( W , ∇ W ) = , (15)where W = α ρ α ρ ρuρvEα , F x ( W ) = α ρ uα ρ uρu + pρvu ( E + p ) u , F y ( W ) = α ρ vα ρ vρuvρv + p ( E + p ) v , Σ ( W , ∇ W ) = u · ∇ α . (16)If the domain is discretized uniformly into a Cartesian grid with N x × N y grid points,we have the domain covered by cells I i,j = (cid:2) x i − / , x i +1 / (cid:3) × (cid:2) y j − / , y j +1 / (cid:3) for 1 ≤ i ≤ N x , 1 ≤ j ≤ N y , where the grid midpoints are given by: x i + = x a + i ∆ x, y j + = y a + j ∆ y (17)6nd ∆ x = x b − x a N x , ∆ y = y b − y a N y . (18)To obtain the numerical scheme with an exact or approximate Riemann solver, thenumerical discretizations in different directions are treated independently. The flux in y direction has similar formulation. Therefore, only the numerical discretization in x direction is discussed in details in this section. The discretization in x direction is conducted by considering a generalized Riemannproblem with a planar discontinuity initially at each grid midpoint in the x direction, asshown in figure 1. The generalized Riemann problem is reduced to a quasi-1D problemdue to assumed homogeneity in other directions. The reduced quasi-1D generalizedRiemann problem at midpoint x = x i +1 / between grid cells located at ( x i , y j ) and( x i +1 , y j ) is hence formulated as: ∂ t W + ∂ x F ( W ) + Σ ( W , ∂ x W ) = 0 , W ( x, t = 0) = (cid:40) W L , x i + < , W R , x i + ≥ . (19)where superscript “ x ” in F x and index “ j ” are dropped for convenience. For first orderaccurate spatial approximation, W L = W i and W R = W i +1 . Higher order spatialapproximation can be obtained from high-order interpolation to construct a high-orderscheme which will be discussed in another section. The first order solution of Riemannproblem given by equation (19) is self-similar: W ( x, t ) = R ( x/t, W i , W i +1 ).The exact solution to the generalized Riemann problem is computationally expen-sive and challenging to obtain. Approximate Riemann solvers can be used to provideapproximate solutions in a less expensive way in Godunov-type schemes. Assuming theapproximate waves generated at the two midpoints x i ± / do not interact under suit-able Courant–Friedrichs–Lewy (CFL) condition as shown in figure 2, the approximatenumerical solution at grid cell i is the cell-averaged value: W n +1 i = 1∆ x (cid:90) x i x i − R (cid:18) x − x i − ∆ t , W ni − , W ni (cid:19) dx + 1∆ x (cid:90) x i + 12 x i R (cid:18) x − x i + ∆ t , W ni , W ni +1 (cid:19) dx. (20)Specifically since the last equation is an advection equation, the approximate solutionof the volume fraction is given by: α n +1 i = (cid:20) α ni − u + ,i − ∆ t ∆ x (cid:0) α ni − α ni − (cid:1)(cid:21) + (cid:20) α ni − u − ,i + ∆ t ∆ x (cid:0) α ni +1 − α ni (cid:1)(cid:21) (21) α n +1 i = α ni − u + ,i − ∆ t ∆ x (cid:0) α ni − α ni − (cid:1) − u − ,i + ∆ t ∆ x (cid:0) α ni +1 − α ni (cid:1) , (22)7here u + ,i − / = max { , u ∗ ,i − / } and u − ,i +1 / = min { , u ∗ ,i +1 / } . u ∗ ,i ± / are theapproximate material wave speeds at edges x i ± / for the advection equation. x i − x i + x i + x i x i +1 t n t n +1 xt Figure 1: Illustration of a scheme with an approximate Riemann solver. Waves are generated fromdiscontinuities initially located at grid midpoints. x i − x i x i + t n t n +1 xt Figure 2: Approximate waves generated at the left and right edges of a cell I i = (cid:2) x i − / , x i +1 / (cid:3) . Thered solid lines and green dotted lines represent the locations of the approximate acoustic waves andcontact waves respectively. The HLLC discretization in a particular direction for a multi-dimensional problemcan be approximated by the solutions of a quasi-1D generalized Riemann problem inthat direction. The approximate solutions, W HLLC , of a quasi-1D generalized Riemannproblem in x direction with an initial planar discontinuity, is illustrated in figure 3. Theapproximate solution to the generalized Riemann problem contains three discontinuities:one contact wave and two acoustic waves. The speed of the contact wave is denotedby s ∗ while the smallest and largest acoustic wave speeds are represented by s L and s R respectively. The material wave speed for the advection equation is chosen as the contactwave speed u ∗ = s ∗ .The HLLC approximate solutions in the four different regions separated by the threediscontinuities are given by: W HLLC = W L , if s L > , W ∗ ,L , if s L ≤ < s ∗ , W ∗ ,R , if s ∗ ≤ ≤ s R , W R , if s R < , (23)8here L and R are the left and right states respectively at a midpoint. With K = L or R , the star state for the five-equation model for a 2D problem is given by: W ∗ ,K = χ ∗ ,K ( α ρ ) K χ ∗ ,K ( α ρ ) K χ ∗ ,K ρ K s ∗ χ ∗ ,K ρ K v K χ ∗ ,K (cid:104) E K + ( s ∗ − u K ) (cid:16) ρ K s ∗ + p K s K − u K (cid:17)(cid:105) α ,K (24) χ ∗ K is defined as: χ ∗ K = s K − u K s K − s ∗ . (25)We use the wave speeds suggested by Einfeldt et al. [11]: s L = min (¯ u − ¯ c, u L − c L ) , s R = max (¯ u + ¯ c, u R + c R ) , (26)where ¯ u and ¯ c are the arithmetic averages from the left and right states. For instance, ¯ c is the average of c L and c R . Following Batten et al. [3], the wave speed in the star regionis given by: s ∗ = p R − p L + ρ L u L ( s L − u L ) − ρ R u R ( s R − u R ) ρ L ( s L − u L ) − ρ R ( s R − u R ) . (27) x i + x L x R W L W R W ∗ ,L W ∗ ,R s L s R s ∗ t n t n +1 xt Figure 3: HLLC wave diagram at interface x i +1 / . The red solid lines represent the locations of theapproximate acoustic waves and the green dotted line represent the location of the contact wave. We now introduce a flux-source form that is convenient for the derivation of the flux-based numerical discretization for the non-conservative system of equations and also theextension for high-order methods. The equation given by (15) can be rewritten as: ∂ t W + ∂ x G x ( W ) + ∂ y G y ( W ) = S ( W , ∇ W ) , (28)9here G x ( W ) = α ρ uα ρ uρu + pρvu ( E + p ) uα u , G y ( W ) = α ρ vα ρ vρuvρv + p ( E + p ) vα v , S ( W , ∇ W ) = α ∇ · u . (29)The relation between fluxes F x and G x and that between fluxes F y and G y are givenby: G x = F x + (0 0 0 0 0 f xα ) T , (30) G y = F y + (0 0 0 0 0 f yα ) T , (31)where f xα = α u and f yα = α v .The fully discretized form of equation (28) with first order accurate forward Eulertime integration is given by: W n +1 i,j − W ni,j ∆ t + ˆ G xi + ,j − ˆ G xi − ,j ∆ x + ˆ G yi,j + − ˆ G yi,j − ∆ y = ˆ S i,j , (32)where ˆ S i,j = α n ,i,j (cid:18) ˆ u i + 12 ,j − ˆ u i − ,j ∆ x + ˆ v i,j + 12 − ˆ v i,j − ∆ y (cid:19) . (33)Only numerical approximation in x direction is discussed in the following part as dis-cretization in y direction is similar. The first order accurate approximation of ˆ G xi ± ,j with the HLLC solutions is given by:ˆ G xi − ,j = G x, HLLC i − ,j ( W i − ,j , W i,j )= F x, HLLC i − ,j ( W i − ,j , W i,j ) + (0 0 0 0 0 f x, HLLC α,i − ,j ( W i − ,j , W i,j )) T , (34)ˆ G xi + ,j = G x, HLLC i + ,j ( W i,j , W i +1 ,j )= F x, HLLC i + ,j ( W i,j , W i +1 ,j ) + (0 0 0 0 0 f x, HLLC α,i + ,j ( W i,j , W i +1 ,j )) T , (35)where F x, HLLC i − ,j ( W i − ,j , W i,j ) = F x, HLLC (cid:0) R HLLC (0 , W i − ,j , W i,j ) (cid:1) , (36) F x, HLLC i + ,j ( W i,j , W i +1 ,j ) = F x, HLLC (cid:0) R HLLC (0 , W i,j , W i +1 ,j ) (cid:1) . (37)10he conservative HLLC fluxes F x, HLLC i ± ,j can be obtained with the divergence theorem [3]: F x, HLLC (cid:0) R HLLC (0 , W L , W R ) (cid:1) =1 + sign( s ∗ )2 ( F L + s − ( W ∗ ,L − W L )) + 1 − sign( s ∗ )2 ( F R + s + ( W ∗ ,R − W R )) , (38)where s − = min (0 , s L ) , s + = max (0 , s R ) . (39)Note that the last component of F x, HLLC for the advection equation is zero. The dis-cretization of the advection equation is contributed by ˆ u i ± ,j and f x, HLLC α,i ± ,j which aregiven by the first order accurate approximations as:ˆ u i − ,j = u ∗ ,i − ,j = u HLLC ∗ ( W i − ,j , W i,j ) = s ∗ ( W i − ,j , W i,j ) , (40)ˆ u i + ,j = u ∗ ,i + ,j = u HLLC ∗ ( W i,j , W i +1 ,j ) = s ∗ ( W i,j , W i +1 ,j ) , (41)and f x, HLLC α,i − ,j ( W i − ,j , W i,j ) = 1 + sign( s ∗ ,i − ,j )2 ( α ,i − ,j s ∗ ,i − ,j )+ 1 − sign( s ∗ ,i − ,j )2 ( α ,i,j s ∗ ,i − ,j ) , (42) f x, HLLC α,i + ,j ( W i,j , W i +1 ,j ) = 1 + sign( s ∗ ,i + ,j )2 ( α ,i,j s ∗ ,i + ,j )+ 1 − sign( s ∗ ,i + ,j )2 ( α ,i +1 ,j s ∗ ,i + ,j ) . (43)The expressions given above form the first order accurate solution of volume fractiongiven by equation (22).Finally, the non-conservative flux ˆ G x, ± i ∓ ,j (similarly for ˆ G y, ± i,j ∓ ) is introduced:ˆ G x, ± i ∓ ,j = ˆ G xi ∓ ,j − α ,i,j (0 0 0 0 0 ˆ u i ∓ ,j ) T , (44)where equation (32) can be simplified to: W n +1 i,j − W ni,j ∆ t + ˆ G x, − i + ,j − ˆ G x, + i − ,j ∆ x + ˆ G y, − i,j + − ˆ G y, + i,j − ∆ y = . (45)Note that the components of ˆ G x, ± i ∓ ,j (or ˆ G x, HLLC , ± i ∓ ,j more precisely) for all conservativeequations, except the last advection equation, are conservative numerical fluxes and areequivalent to the corresponding components of F x, HLLC i ∓ ,j . For a quasi-1D problem, it is shown in equation (20) that the solution update isthe convex averaging of the exact or approximate solutions to the generalized Riemann11roblem. Therefore, the HLLC Riemann solver gives physically admissible solution if allstates generated are physically admissible using Jensen’s inequality for integral equations.Here, the left star state is considered and the right star state can be proved to bephysically admissible by symmetry.With the definition of s ∗ given by equation (27), it can be shown that s L < s ∗ [3].Also, since s L = min (¯ u − ¯ c, u L − c L ), s L < u L . As a result, the partial densities in thestar state are positive: ( α ρ ) ∗ ,L = s L − u L s L − s ∗ ( α ρ ) L ≥ , (46)( α ρ ) ∗ ,L = s L − u L s L − s ∗ ( α ρ ) L ≥ . (47)Since ρ ∗ ,L = ( α ρ ) ∗ ,L + ( α ρ ) ∗ ,L , ρ ∗ ,L = s L − u L s L − s ∗ ρ L ≥ . (48)The positivity of mixture and partial densities implies all mass fractions are boundedbetween zero and one. As for the volume fraction, since α , ∗ ,L = α ,L ,0 ≤ α , ∗ ,L ≤ . (49)The only remaining requirement is ( ρe ) ∗ ,L > p ∞ L for positive squared speed of sound.From the definition of ( ρe ) ∗ ,L :( ρe ) ∗ ,L = E ∗ ,L − | ρ u ∗ ,L | ρ ∗ ,L = E ∗ ,L −
12 ( ρu ) ∗ ,L + ( ρv ) ∗ ,L ρ ∗ ,L (50)( ρe ) ∗ ,L = s L − u L s L − s ∗ E L + s L − u L s L − s ∗ ( s ∗ − u L ) (cid:18) ρ L s ∗ + p L s L − u L (cid:19) − s L − u L s L − s ∗ ρ L (cid:0) s ∗ + v L (cid:1) . (51)Therefore, we require the following inequality: s L − u L s L − s ∗ E L + s L − u L s L − s ∗ ( s ∗ − u L ) (cid:18) ρ L s ∗ + p L s L − u L (cid:19) − s L − u L s L − s ∗ ρ L (cid:0) s ∗ + v L (cid:1) > p ∞ L (52) E L + ( s ∗ − u L ) (cid:18) ρ L s ∗ + p L s L − u L (cid:19) − ρ L (cid:0) s ∗ + v L (cid:1) > p ∞ L ( s L − u L − s ∗ + u L )( s L − u L ) (53)( ρe ) L + 12 ρ L (cid:0) u L + v L (cid:1) + ( s ∗ − u L ) (cid:18) ρ L s ∗ + p L s L − u L (cid:19) − ρ L (cid:0) s ∗ + v L (cid:1) > p ∞ L ( s L − u L − s ∗ + u L )( s L − u L ) (54)12 ρ L ( s ∗ − u L ) − p L + p ∞ L u L − s L ( s ∗ − u L ) + ( ρ L e L − p ∞ L ) > . (55)12et β = s ∗ − u L , then the inequality above is a quadratic function of β . We can showthat this quadratic has no real roots by ensuring the discriminant is negative. That is, (cid:18) p L + p ∞ L u L − s L (cid:19) − ρ L ( ρ L e L − p ∞ L ) < . (56)This implies that we require: s L < u L − p L + p ∞ L (cid:112) ρ L ( ρ L e L − p ∞ L ) . (57)It should be noted that: p L + p ∞ L (cid:112) ρ L ( ρ L e L − p ∞ L ) = (cid:115) ( γ L − p L + p ∞ L )2 ρ L < (cid:115) γ L ( p L + p ∞ L ) ρ L = c L . Since s L = min (¯ u − ¯ c, u L − c L ), the constraint on s L given by equation (57) is alreadysatisfied and we have proved the solutions given by the left star state have positive partialdensities and squared speed of sound. Besides, the volume fractions are bounded. Thus,the HLLC Riemann solver is positivity- and boundedness-preserving.Equation (45) can be re-written as: W n +1 i,j = σ x (cid:104) W ni,j + λ x (cid:16) ˆ G x, + i − ,j − ˆ G x, − i + ,j (cid:17)(cid:105) + σ y (cid:104) W ni,j + λ y (cid:16) ˆ G y, + i,j − − ˆ G y, − i,j + (cid:17)(cid:105) , (58)where λ x = ∆ t/ (∆ xσ x ) and λ y = ∆ t/ (∆ yσ y ). σ x and σ y are partitions of the contri-bution in the x and y directions respectively where σ x + σ y = 1. They can be definedas [20]: σ x = τ x τ x + τ y , σ y = τ y τ x + τ y , τ x = ( | u | + c ) max ∆ x , τ y = ( | v | + c ) max ∆ y , (59)such that 0 < σ x < < σ y <
1. If the time-step size ∆ t is given by a chosen CFLnumber, CF L , with the following equation:∆ t = CF Lτ x + τ y , (60)one has the relations for the equivalent 1D time step sizes in different directions, ∆ t x and ∆ t y : λ x = CF L ( | u | + c ) max = ∆ t x ∆ x , (61) λ y = CF L ( | v | + c ) max = ∆ t y ∆ y . (62)13e can define: W xi,j = W ni,j + λ x (cid:16) ˆ G x, + i − ,j − ˆ G x, − i + ,j (cid:17) , (63) W yi,j = W ni,j + λ y (cid:16) ˆ G y, + i,j − − ˆ G y, − i,j + (cid:17) . (64)Since ˆ G x, ± i ∓ ,j are obtained from the quasi-1D HLLC solutions in the x direction, theapproximate waves from the edges at x i ± do not interact if the CFL condition, CF L ≤ .
5, is satisfied (same for ˆ G y, ± i,j ∓ in the y direction). Thus, equation (20) is satisfied for W xi,j (similar for W yi,j in y direction). The first part of the RHS of equation (20) is thesolution in the half cell from x i − to x i at t + ∆ t x and the second part is the solutionin another half cell. Using the finite volume approach on the left and right half cells, wewill get:1∆ x (cid:90) x i x i − R (cid:18) x − x i − ∆ t x , W ni − ,j , W ni,j (cid:19) dx = 12 (cid:104) W ni,j + 2 λ x (cid:16) ˆ G x, + i − ,j − F xi,j (cid:17)(cid:105) , (65)1∆ x (cid:90) x i + 12 x i R (cid:18) x − x i + ∆ t x , W ni,j , W ni +1 ,j (cid:19) dx = 12 (cid:104) W ni,j − λ x (cid:16) ˆ G x, − i + ,j − F xi,j (cid:17)(cid:105) . (66)If we define: W x, − i,j = 2∆ x (cid:90) x i x i − R (cid:18) x − x i − ∆ t , W ni − ,j , W ni,j (cid:19) dx, (67) W x, + i,j = 2∆ x (cid:90) x i + 12 x i R (cid:18) x − x i + ∆ t , W ni,j , W ni +1 ,j (cid:19) dx, (68)we will get: W xi,j = 12 W x, − i,j + 12 W x, + i,j . (69) W x, ± i,j are at physically admissible states since they are convex averaging of the approxi-mate HLLC solutions. This also means W xi,j is also physically admissible. Equation (58)becomes: W n +1 i,j = σ x W xi,j + σ y W yi,j . (70)Therefore, W n +1 i,j is also physically admissible since it is a convex combination of W xi,j and W yi,j .
4. Incremental-stencil WCNS
In this section, a high-order finite difference scheme for discretizing equation (28) isintroduced. The high-order scheme belongs to the family of weighted compact nonlinearschemes (WCNSs) which is a variant of the WENO schemes for discontinuity-capturing.It was first proposed by Deng and Zhang [10] in which compact (spatially implicit) fi-nite difference schemes are combined with WENO interpolation. Since then, WCNSsare extended to higher order of accuracy [31, 51]. In principle, WCNSs can be used14ith both explicit or compact finite difference schemes. Nonomura and Fujii [29] sug-gested that explicit finite difference schemes are more efficient and later also proposeda family of robust explicit midpoint-and-node-to-node finite difference schemes [30]. Inthis section, an explicit WCNS with the explicit hybrid cell-midpoint and cell-node finitedifference scheme [9, 50] and nonlinear interpolation adapted from incremental-stencilreconstruction [48] for finite volume WENO scheme is presented.The semi-discretized finite difference form of equation (28) is given by: ∂ W ∂t (cid:12)(cid:12)(cid:12)(cid:12) i,j + (cid:100) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j + (cid:100) ∂ G y ∂y (cid:12)(cid:12)(cid:12)(cid:12) i,j = ˆ S i,j , (71)where ˆ S i,j is given by: ˆ S i,j = α ,i,j (cid:18) (cid:99) ∂u∂x (cid:12)(cid:12)(cid:12) i,j + (cid:99) ∂v∂y (cid:12)(cid:12)(cid:12) i,j (cid:19) . (72)A high-order discretization is considered consistent (and conservative for any conservativeequation) if: (cid:100) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = ˆ G xi + ,j − ˆ G xi − ,j ∆ x , (cid:100) ∂ G y ∂y (cid:12)(cid:12)(cid:12)(cid:12) i,j = ˆ G yi,j + − ˆ G yi,j − ∆ y , (73) (cid:99) ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = ˆ u i + ,j − ˆ u i − ,j ∆ x , (cid:99) ∂v∂y (cid:12)(cid:12)(cid:12)(cid:12) i,j = ˆ v i,j + − ˆ v i,j − ∆ y . (74)Therefore, equation (71) becomes: ∂ W ∂t (cid:12)(cid:12)(cid:12)(cid:12) i,j + ˆ G xi + ,j − ˆ G xi − ,j ∆ x + ˆ G yi,j + − ˆ G yi,j − ∆ y = ˆ S i,j , (75)Equation (75) looks as same as equation (32) if it is further discretized in time withforward Euler method but ˆ G xi ± ,j , ˆ G yi,j ± and ˆ S i,j (composed of ˆ u i ± ,j and ˆ v i,j ± ) arein high-order accurate approximations and are given by the WCNS introduced in thissection. It should be noted that equation (75) is a conservative discretization for allconservative equations, except the last advection equation that is given by: ∂α ∂t (cid:12)(cid:12)(cid:12)(cid:12) i,j + ˆ f xα,i + ,j − ˆ f xα,i − ,j ∆ x + ˆ f yα,i,j + − ˆ f yα,i,j − ∆ y = α ,i,j (cid:18) ˆ u i + ,j − ˆ u i − ,j ∆ x + ˆ v i,j + − ˆ v i,j − ∆ y (cid:19) . (76)15 .1. Explicit hybrid cell-midpoint and cell-node scheme The sixth order accurate explicit scheme from the hybrid cell-midpoint and cell-node compact scheme (HCS) [9] family is used for the approximation of the first orderderivatives ∂ G x /∂x | i,j and ∂ G y /∂y | i,j . The sixth order explicit HCS formulation isgiven by: ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j ≈ (cid:100) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = 1∆ x (cid:20) ψ (cid:16) ˜ G xi + ,j − ˜ G xi − ,j (cid:17) − ψ − (cid:0) G xi +1 ,j − G xi − ,j (cid:1) + 35 ψ − (cid:0) G xi +2 ,j − G xi − ,j (cid:1) − ψ − (cid:0) G xi +3 ,j − G xi − ,j (cid:1)(cid:21) . (77)If we replace ˜ G xi ± ,j with the exact fluxes, (cid:100) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j − (cid:18) ψ − (cid:19) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j ∆ x − (cid:18) ψ − (cid:19) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j ∆ x + O (cid:0) ∆ x (cid:1) . (78)If ψ = 256 / (cid:100) ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j − ∂ G x ∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j ∆ x + O (cid:0) ∆ x (cid:1) . (79) ψ = 256 /
175 is adopted in this work.Any central explicit or compact finite difference scheme can be rewritten into theflux-difference forms given by equation (73) and it is derived in [44]. Following thatwork, the HCS given by equation (77) has implied reconstructed flux ˆ G xi +1 / ,j given by:ˆ G xi + ,j = ψ ˜ G xi + ,j − (cid:18) ψ − (cid:19) (cid:0) G xi,j + G xi +1 ,j (cid:1) + (cid:18) ψ − (cid:19) (cid:0) G xi − ,j + G xi +2 ,j (cid:1) − (cid:18) ψ − (cid:19) (cid:0) G xi − ,j + G xi +3 ,j (cid:1) . (80)Note that the equation above is also used for reconstructing the flux ˆ f xα,i + ,j in theadvection equation.High-order finite difference approximations of the velocity components are also re-quired for ˆ S i,j . Following the idea of [49], the numerical derivatives of the velocitycomponents are also given by the same finite difference scheme as the flux derivatives: ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j ≈ (cid:99) ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = 1∆ x (cid:20) ψ (cid:16) ˜ u i + ,j − ˜ u i − ,j (cid:17) − ψ − u i +1 ,j − u i − ,j )+ 35 ψ − u i +2 ,j − u i − ,j ) − ψ − u i +3 ,j − u i − ,j ) (cid:21) . (81)16he implied reconstructed velocity component ˆ u i +1 / ,j is given by:ˆ u i + ,j = ψ ˜ u i + ,j − (cid:18) ψ − (cid:19) ( u i,j + u i +1 ,j ) + (cid:18) ψ − (cid:19) ( u i − ,j + u i +2 ,j ) − (cid:18) ψ − (cid:19) ( u i − ,j + u i +3 ,j ) . (82)The discretizations for flux and velocity component derivatives in y direction aresimilar. High-order accurate approximations are required for ˜ G xi +1 / ,j and ˜ u i +1 / ,j toform high-order discretization for equation (75), which are discussed in the followingsections. Finally, by using equation (44), equation (71) can be rewritten as: ∂ W ∂t (cid:12)(cid:12)(cid:12)(cid:12) i + ˆ G x, − i + ,j − ˆ G x, + i − ,j ∆ x + ˆ G y, − i,j + − ˆ G y, + i,j − ∆ y = . (83) In WCNSs, the fluxes at the midpoints are obtained with aid of explicit nonlinearinterpolations, which can also be interpreted as nonlinear filtering processes to avoidspurious oscillations near shocks and other discontinuities. For simplicity, the implemen-tation details of a WCNS is explained with a 1D version of equation (83) in this and thefollowing sub-section. Thus, superscript “ x ” in G x is dropped for convenience.For the 1D five-equation model, the algorithm to obtain the high-order fluxes ˆ G ± i +1 / with the WCNS approach is given below:1. Convert all W i in the stencils of left-biased and right-biased nonlinear WENOinterpolations to primitive variable vectors V i .2. Perform characteristic decomposition by transforming all V i in the stencils of inter-polations to characteristic variable vectors U i with the projection matrix R − i +1 / : U i = R − i +1 / V i .3. Compute ˜ U L and ˜ U R at each midpoint with U i using left-biased and right-biasednonlinear WENO interpolations respectively.4. Transform ˜ U L and ˜ U R to ˜ V L and ˜ V R with the projection matrix R i +1 / : ˜ V L = R i +1 / ˜ U L , and ˜ V R = R i +1 / ˜ U R .5. Convert ˜ V L and ˜ V R to ˜ W L and ˜ W R .6. Compute the high-order flux and velocity at each midpoint with the Riemannsolver. If the HLLC Riemann solver is used: ˜ G i + = G HLLC (cid:16) ˜ W L , ˜ W R (cid:17) and˜ u i + = u HLLC ∗ (cid:16) ˜ W L , ˜ W R (cid:17) .7. Compute the flux and velocity at the nodes: G i = G ( W i ) and u i = u ( W i ).8. Reconstruct the flux ˆ G i + and velocity ˆ u i + at the midpoints using the flux dif-ferencing approach (equations (80) and (82) in this work).9. Compute ˆ G ± i +1 / using equation (44):ˆ G − i + = ˆ G i + − α ,i (0 0 0 0 ˆ u i + ) T , ˆ G + i + = ˆ G i + − α ,i +1 (0 0 0 0 ˆ u i + ) T .
17n this work, only the incremental-stencil WENO interpolation of left-biased midpointvalues is presented. The interpolation of right-biased midpoint values is similar due tosymmetry and can be obtained by flipping the stencils and corresponding coefficients.The projection matrices for transformation between primitive variables and characteristicvariables can be found in Appendix A. The projection matrices R i +1 / are computed atmidpoints x i +1 / with the arithmetic averages of partial densities, mixture density andspeed of sound at x i and x i +1 . The finite volume WENO scheme with the incremental-stencil reconstruction (WENO-IS) was proposed by Wang et al. [48]. The finite volume WENO-IS is robust for com-pressible multi-phase problems with shocks and is also accurate for those problems dueto high-order WENO reconstruction with the use of HLLC Riemann solver, which iswell-known for its accuracy in capturing material interfaces. However, in general a finitevolume WENO scheme is more expensive compared with finite difference WCNS andWENO schemes with similar orders of accuracy for multi-dimensional simulations. Thisis due to the fact that a finite volume scheme requires multi-dimensional reconstruc-tions to obtain point values at many Gaussian points on the cell boundaries from cellaverages when one desires third or high order of accuracy [41, 46, 8]. Generally, finitedifference WENO schemes or WCNSs in explicit forms are four times cheaper in 2Dand nine times cheaper in 3D compared to finite volume WENO schemes since multi-dimensional reconstructions are not required for the former schemes [36]. While the costsof WENO reconstruction in finite difference WENO schemes and WENO interpolation infinite difference WCNS methods are similar, a finite difference WENO scheme can onlybe used with flux-vector splitting methods, such as Lax–Friedrichs flux splitting, whenhigh-order of accuracy is desired. The use of flux-difference splitting methods such asRiemann solvers in a finite difference WENO scheme for multi-dimensional simulationsdegenerates to only second order but a WCNS can still maintain high order of accuracywhen used with a Riemann solver. Therefore, in this work we propose a WCNS witha WENO interpolation adapted from the robust incremental-stencil WENO reconstruc-tion such that the HLLC Riemann solver can be applied for upwinding while the overallscheme is still high-order accurate and efficient.The incremental-stencil (IS) interpolation approximates the midpoint values by non-linear combination of linearly interpolated values from four different sub-stencils, S – S (shown in figure 4). The interpolated values at the midpoints ˜ u j + from the four differentsub-stencils are given by: S : ˜ u i + = 12 ( u i + u i +1 ) , (84) S : ˜ u i + = 12 ( − u i − + 3 u i ) , (85) S : ˜ u i + = 18 (3 u i + 6 u i +1 − u i +2 ) , (86) S : ˜ u i + = 18 (3 u i − − u i − + 15 u i ) . (87)The interpolated values from stencils S and S are second order accurate and thosefrom stencils S and S are third order accurate. The variable u can either be fluxes,18onservative variables, primitive variables or variables that are projected to the charac-teristic fields. In this work, the primitive variables projected to the characteristic fieldsare employed in the interpolation process. i + i − i − i + i + i − i − i i + 1 i + 2 i + 3 S S S S S S S is in grey since it is not considered in the actual interpolation. A fifth order 5-point linear interpolation from S can be obtained from linear combi-nation of the lower order interpolations: S : ˜ u i + = (cid:88) k =0 d k ˜ u ki + , (88)where the linear weights are given by: d = 1532 , d = 532 , d = 516 , d = 116 . (89)The expanded form of the linear interpolation is given by: S : ˜ u i + = 1128 (3 u i − − u i − + 90 u i + 60 u i +1 − u i +2 ) . (90)The 5-point linear interpolation scheme may generate spurious oscillations due toGibbs phenomenon or is even unstable near shocks or discontinuities. To improve ro-bustness, a nonlinear WENO interpolation is suggested in this work by replacing thelinear weights with the incremental-stencil (IS) nonlinear weights [48]. The nonlinearweights have the following form: ˜ u i + = (cid:88) k =0 ω k ˜ u ki + . (91)19he IS nonlinear weights are given by [48]: ω k = η k (cid:80) s =0 η s , (92) η k = d k (cid:16) τ β k + (cid:15) · τ β + (cid:15) (cid:17) , if k < ,d k (cid:16) τ β k + (cid:15) (cid:17) , otherwise . (93)The smoothness indicators are defined by: β k = (cid:82) x i + 12 x i − ∆ x (cid:0) ∂∂x ˜ u k ( x ) (cid:1) dx, k = 0 , , (cid:80) l =1 (cid:82) x i + 12 x i − ∆ x l − (cid:16) ∂ l ∂x l ˜ u k ( x ) (cid:17) dx, otherwise (including k = 01) , (94)where ˜ u k ( x ) are the Lagrange interpolating polynomials from stencils S k . The nonlinearweights require smoothness indicator computed with the Lagrange interpolation poly-nomial from stencil S besides those from S – S . The integrated forms of smoothnessindicators are given by: β = ( u i − u i +1 ) , (95) β = ( u i − − u i ) , (96) β = 1312 ( u i − − u i + u i +1 ) + 14 ( u i − − u i +1 ) , (97) β = 1312 ( u i − u i +1 + u i +2 ) + 14 (3 u i − u i +1 + u i +2 ) , (98) β = 1312 ( u i − − u i − + u i ) + 14 ( u i − − u i − + 3 u i ) . (99)The reference smoothness indicator, τ , is defined by [12]: τ = (cid:88) l =3 (cid:90) x i + 12 x i − ∆ x l − (cid:18) ∂ l ∂x l ˜ u ( x ) (cid:19) dx (100)where ˜ u ( x ) is the Lagrange interpolating polynomial from stencil S . The integratedform of the reference smoothness indicator is given by: τ = 1312 ( u i +2 − u i +1 + 6 u i − u i − + u i − ) + 14 ( u i +2 − u i +1 + 2 u i − − u i − ) . (101)The robustness of incremental-stencil WENO interpolation presented above is as highas the WENO incremental-stencil reconstruction, as both of them can choose one ofthe 2-point stencils when there are closely located discontinuities, as explained in [48].In Appendix B, it is proved that the IS nonlinear interpolation with the HCS finitedifferencing is fifth order accurate for a 1D scalar hyperbolic conservation law with perfectupwinding, provided not at critical points. This WCNS is termed WCNS-IS in this work.20 . Positivity- and boundedness-preserving limiting procedures Similar to the finite volume WENO-IS, the finite difference WCNS-IS is generallyrobust in capturing shocks and material interfaces in multi-phase flow simulations. How-ever, it cannot be guaranteed that these high-order schemes are free from numericalfailures due to negative squared speed of sound, partial densities (hence also out-of-bounds mass fractions), and out-of-bounds volume fractions. In this section, positivity-and boundedness-preserving limiting procedures are introduced to improve the robust-ness of the WCNS-IS. The procedures are conservative for the corresponding equationsin the system that are conservative, i.e. except the volume fraction advection equation.In the algorithm of WCNS, there are two stages where positivity and boundednesscan be violated. The first stage is the WENO interpolation step for the left-biasedand right-biased interpolated conservative variable vectors ( ˜ W L and ˜ W R ), where theinterpolated conservative variables may not be physically admissible. Another stageis the flux reconstruction step of ˆ G ± i +1 / using the Riemann solver and the high-orderHCS finite difference scheme. The positivity- and boundedness-preserving interpolationlimiter and flux limiter are introduced in this section to respectively deal with the twoissues mentioned. The incremental-stencil WENO interpolation is robust but it is still not positivity-preserving for partial densities and squared speed of sound. It is also not boundedness-preserving for volume fractions. However, the first order interpolation with the leftand right node values are in the admissible state set and are positivity-preserving andboundedness-preserving. Therefore, the high-order WENO interpolation can be limitedwith a convex combination of itself and the first order interpolation. For simplicity, thissub-section only discusses the limiting procedures for the left-biased WENO interpolationin the x direction for the governing equations. Hence, the subscripts “ L ” and “ j ” areomitted.The first stage of the interpolation limiting procedures is to obtain a limited conserva-tive variable vector with positive partial densities at each midpoint, ˜ W ∗ i + . As describedin algorithm 1, ˜ W ∗ i + is first initialized as the WENO interpolated conservative variablevector. It is then limited for positive partial densities through repeated convex combi-nation of ˜ W ∗ i + with the first order interpolated conservative variable vector for eachphase using a user-defined small tolerance (cid:15) α k ρ k . In the next stage, limiting procedureon the volume fractions can be applied similarly through another set of successive convexcombinations given by algorithm 2 with another user-defined small threshold (cid:15) α k . Afterthis stage, the limited vector ˜ W ∗∗ i + should have all partial densities (including mixturedensity) positive and all volume fractions bounded between (cid:15) α k and 1 − (cid:15) α k . This alsomeans that all mass fractions are bounded between 0 and 1. The remaining quantity tolimit is the squared speed of sound.We can define a helper variable ˜ c where ˜ c = ( ρe − p ∞ ) /ρ such that squared soundspeed and ˜ c are related through c = γ ( γ −
1) ˜ c . As the final stage of the limitingprocedures, ˜ W ∗∗∗ i + with positive ρ ˜ c can be obtained using algorithm 3. This is carriedout through convex combination of ˜ W ∗∗ i + and the first order interpolated conservative21ariable vector with tolerance (cid:15) ρ ˜ c by utilizing the convexity of the admissible set provedwith the Jensen’s inequality. The squared speed of sound of the final limited conservativevariable vector, ˜ W ∗∗∗ i + , is also positive as c = γ ( γ − ρ ˜ c ) /ρ . Note that γ is larger than1 since volume fractions of ˜ W ∗∗ i + are already bounded between 0 and 1. The positivityand boundedness limiting procedures for the right-biased WENO interpolation in the x direction can be performed with W i +1 instead of W i in a similar way. We havechosen the tolerances for the limiting procedures as (cid:15) α k ρ k = 1 . − (cid:15) α k = 1 . − (cid:15) ρ ˜ c = 1 . −
8. Note that in practice, user should make sure the tolerances are chosento be smaller than the minimum values of the corresponding initial fields.Due to numerical round-off, the ˜ W ∗∗∗ i + may still have negative partial densities, nega-tive squared speed of sound, or out-of-bounds volume fractions. Therefore, a hard switchis suggested by setting ˜ W ∗∗∗ i + = W i when α k ρ k (cid:16) ˜ W ∗∗∗ i + (cid:17) < (cid:15) HS α k ρ k , α k (cid:16) ˜ W ∗∗∗ i + (cid:17) < (cid:15) HS α k , or ρ ˜ c (cid:16) ˜ W ∗∗∗ i + (cid:17) < (cid:15) HS ρ ˜ c . The tolerances of the hard switch are chosen as (cid:15) HS α k ρ k = 1 . − (cid:15) HS α k = 1 . −
11, and (cid:15) HS ρ ˜ c = 1 . − Algorithm 1:
Left-biased interpolation limiting procedure for positive partialdensities.Set ˜ W ∗ i + = ˜ W i + ; for k = 1 , dofor all midpoints doif α k ρ k ( W i ) < (cid:15) α k ρ k then θ i + = 0; else if α k ρ k (cid:16) ˜ W ∗ i + (cid:17) < (cid:15) α k ρ k then Solve θ i + from the formula: (cid:16) − θ i + (cid:17) α k ρ k ( W i ) + θ i + α k ρ k (cid:16) ˜ W ∗ i + (cid:17) = (cid:15) α k ρ k ; else θ i + = 1; end Perform convex combination:˜ W ∗ i + = (cid:16) − θ i + (cid:17) W i + θ i + ˜ W ∗ i + ; endend The high-order flux reconstruction step of ˆ G ± i +1 / using a Riemann solver and theHCS finite difference scheme is not positivity- and boundedness-preserving in generaland may cause numerical failures. Therefore, a flux limiter is critical to make sure that22 lgorithm 2: Left-biased interpolation limiting procedure for bounded volumefractions.Set ˜ W ∗∗ i + = ˜ W ∗ i + ; for k = 1 , dofor all midpoints doif α k ( W i ) < (cid:15) α k then θ i + = 0; else if α k (cid:16) ˜ W ∗∗ i + (cid:17) < (cid:15) α k then Solve θ i + from the formula: (cid:16) − θ i + (cid:17) α k ( W i ) + θ i + α k (cid:16) ˜ W ∗∗ i + (cid:17) = (cid:15) α k ; else θ i + = 1; end Perform convex combination:˜ W ∗∗ i + = (cid:16) − θ i + (cid:17) W i + θ i + ˜ W ∗∗ i + ; endendAlgorithm 3: Left-biased interpolation limiting procedure for positive ρ ˜ c . for all midpoints doif ρ ˜ c ( W i ) < (cid:15) ρ ˜ c then θ i + = 0; else if ρ ˜ c (cid:16) ˜ W ∗∗ i + (cid:17) < (cid:15) ρ ˜ c then Solve θ i + from the formula: (cid:16) − θ i + (cid:17) ρ ˜ c ( W i ) + θ i + ρ ˜ c (cid:16) ˜ W ∗∗ i + (cid:17) = (cid:15) ρ ˜ c ; else θ i + = 1; end Perform convex combination (applying Jensen’s inequality):˜ W ∗∗∗ i + = (cid:16) − θ i + (cid:17) W i + θ i + ˜ W ∗∗ i + ; end W n +1 i,j = σ x (cid:18) W x, − i,j + 12 W x, + i,j (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) W xi,j + σ y (cid:18) W y, − i,j + 12 W y, + i,j (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) W yi,j , (102)where W x, ∓ i,j = W ni,j ± λ x (cid:16) ˆ G x, ± i ∓ ,j − F xi,j (cid:17) , (103) W y, ∓ i,j = W ni,j ± λ y (cid:16) ˆ G y, ± i,j ∓ − F yi,j (cid:17) . (104)Note that ˆ G x, ± i ∓ ,j and ˆ G y, ± i,j ∓ here are the high-order reconstructed fluxes in contrastto the first order HLLC fluxes, ˆ G x, HLLC , ± i ∓ ,j and ˆ G y, HLLC , ± i,j ∓ . Since both σ x and σ y arepositive and σ x + σ y = 1, W n +1 i,j is a convex combination of W x, ± i,j and W y, ± i,j . If all fourconservative variable vectors are in the physically admissible set, W n +1 i,j is also in thephysically admissible set. For simplicity, only positivity- and boundedness-preservingflux limiting in x direction for W x, ± i,j is discussed here and hence “ x ” superscript and “ j ”subscript are dropped in the following part.If the first order flux from the approximate Riemann solver is positively preservingsuch that all intermediate states given by the approximate Riemann solutions are in theadmissible state set, such as the HLLC Riemann solver presented in this work, we can firstconstruct the positivity flux limiting procedure for partial densities through the convexcombination of the high-order solution W ± i and the first order solution W HLLC , ± i . Thisfirst stage of the flux limiting procedures to obtain the limited flux ˆ G ∗ , ± i + is detailed inalgorithm 4 with tolerance (cid:15) α k ρ k . After this stage, both solutions W ∗ , ± i time-advancedwith the limited fluxes have all partial densities (including mixture density) positive,where W ∗ , ∓ i = W ni ± λ (cid:16) ˆ G ∗ , ± i ∓ − F i (cid:17) . (105)Note that in the last step of the algorithm, the intention is to hybridize ( ˆ G ∗ , ± i + − F i ) with( ˆ G HLLC , ± i + − F i ) but the F i on both sides of the equation cancel each other. Also, theflux limiting process is conservative for all equations except the last advection equationof volume fraction.In the next stage, we can apply the boundedness flux limiting for volume fractionssimilarly to obtain the limited flux ˆ G ∗∗ , ± i + , which is given in algorithm 5 with thresh-old (cid:15) α k . After this step, the limited solutions W ∗∗ , ± i should have all volume fractionsbounded in addition to positive partial densities (including mixture density), where W ∗∗ , ∓ i = W ni ± λ (cid:16) ˆ G ∗∗ , ± i ∓ − F i (cid:17) . (106)Finally, the positivity flux limiting for squared sound speed is conducted throughthe helper variable ρ ˜ c similarly with the previous sub-section. This step is described24n algorithm 6 with tolerance (cid:15) ρ ˜ c and makes use of the fact that the admissible set ofthe conservative variable vector is convex to obtain the final limited flux ˆ G ∗∗∗ , ± i + . Thesquared speeds of sound c of W ∗∗∗ , ± i , where W ∗∗∗ , ∓ i = W ni ± λ (cid:16) ˆ G ∗∗∗ , ± i ∓ − F i (cid:17) , (107)are limited to be positive since c = γ ( γ − ρ ˜ c ) /ρ .Note that the tolerance values ( (cid:15) α k ρ k , (cid:15) α k , and (cid:15) ρ ˜ c ) are as same as those used inthe positivity- and boundedness-preserving interpolation limiter. Similar to the inter-polation limiter, a hard switch is used by setting ˆ G ∗∗∗ , ± i + = ˆ G HLLC , ± i + if either W ∗∗∗ , + i or W ∗∗∗ , − i +1 are at states that are not bounded by smaller tolerances ( (cid:15) HS α k ρ k = 1 . − (cid:15) HS α k = 1 . −
11, and (cid:15) HS ρ ˜ c = 1 . − The extension of positivity- and boundedness-preserving limiters for high-order fluxwith the high-order strong stability preserving Runge–Kutta (SSPRK) time steppingmethods [40, 14, 13] is trivial since SSPRK time stepping schemes are convex combina-tions of Euler forward steps. However, the upper limit of CFL number is still constrainedby 0.5.
6. Test problems
1D and 2D test problems are conducted with the first order HLLC scheme and thehigh-order WCNS-IS with the positivity- and boundedness-preserving limiters. The twoschemes are termed HLLC and PP-WCNS-IS respectively in this work. All tests involveliquid water and air as an ideal gas. The properties of the fluids are given in table 1. Notethat the fluids satisfy the requirement from the positivity- and boundedness-preservinglimiters as the ratio of specific heats of water is larger than that of air. The three-stagethird order SSPRK scheme (TVDRK3) [40] is used for time stepping for both schemes.The CFL number is chosen to be 0.5 unless constant time step size is used. Whenconstant time step size is used, the corresponding CFL number is always less than 0.5until the end of simulations.
To verify the formal order of accuracy of each scheme, advection of volume fractiondisturbance in a 2D periodic domain [ − ,
1) m × [ − ,
1) m is used as the test problemsimilar to that in [49]. The initial conditions are given by table 2 and the exact solutionsare given by table 3.Simulations using different schemes are conducted up to t = 0 . N x = N y = 8 to N x = N y = 256. All simulations are run with verysmall constant time steps in order to observe the spatial orders of accuracy of differentnumerical schemes. ∆ t/ ∆ x = 4 . e − − is used.25 lgorithm 4: Flux limiting procedure for positive partial densities.Set ˆ G ∗ , ± i + = ˆ G ± i + ( W ∗ , + i = W + i , W ∗ , − i +1 = W − i +1 ); for k = 1 , dofor all midpoints doif α k ρ k (cid:16) W HLLC , + i (cid:17) < (cid:15) α k ρ k then θ + i + = 0; else if α k ρ k (cid:0) W ∗ , + i (cid:1) < (cid:15) α k ρ k then Solve θ + i + from the formula: (cid:16) − θ + i + (cid:17) α k ρ k (cid:16) W HLLC , + i (cid:17) + θ + i + α k ρ k (cid:0) W ∗ , + i (cid:1) = (cid:15) α k ρ k ; else θ + i + = 1; endif α k ρ k (cid:16) W HLLC , − i +1 (cid:17) < (cid:15) α k ρ k then θ − i + = 0; else if α k ρ k (cid:0) W ∗ , − i +1 (cid:1) < (cid:15) α k ρ k then Solve θ − i + from the formula: (cid:16) − θ − i + (cid:17) α k ρ k (cid:16) W HLLC , − i +1 (cid:17) + θ − i + α k ρ k (cid:0) W ∗ , − i +1 (cid:1) = (cid:15) α k ρ k ; else θ − i + = 1; end Set θ i + = min (cid:16) θ + i + , θ − i + (cid:17) ;Perform convex combination:ˆ G ∗ , + i + = (cid:16) − θ i + (cid:17) ˆ G HLLC , + i + + θ i + ˆ G ∗ , + i + ;ˆ G ∗ , − i + = (cid:16) − θ i + (cid:17) ˆ G HLLC , − i + + θ i + ˆ G ∗ , − i + ; endend Table 4 shows the L errors and the computed rates of convergence of volume fraction α respectively by the two schemes at t = 0 . L error of volume fraction is26 lgorithm 5: Flux limiting procedure for bounded volume fractions.Set ˆ G ∗∗ , ± i + = ˆ G ∗ , ± i + ( W ∗∗ , + i = W ∗ , + i , W ∗∗ , − i +1 = W ∗ , − i +1 ); for k = 1 , dofor all midpoints doif α k (cid:16) W HLLC , + i (cid:17) < (cid:15) α k then θ + i + = 0; else if α k (cid:0) W ∗∗ , + i (cid:1) < (cid:15) α k then Solve θ + i + from the formula: (cid:16) − θ + i + (cid:17) α k (cid:16) W HLLC , + i (cid:17) + θ + i + α k (cid:0) W ∗∗ , + i (cid:1) = (cid:15) α k ; else θ + i + = 1; endif α k (cid:16) W HLLC , − i +1 (cid:17) < (cid:15) α k then θ − i + = 0; else if α k (cid:0) W ∗∗ , − i +1 (cid:1) < (cid:15) α k then Solve θ − i + from the formula: (cid:16) − θ − i + (cid:17) α k (cid:16) W HLLC , − i +1 (cid:17) + θ − i + α k (cid:0) W ∗∗ , − i +1 (cid:1) = (cid:15) α k ; else θ − i + = 1; end Set θ i + = min (cid:16) θ + i + , θ − i + (cid:17) ;Perform convex combination:ˆ G ∗∗ , + i + = (cid:16) − θ i + (cid:17) ˆ G HLLC , + i + + θ i + ˆ G ∗∗ , + i + ;ˆ G ∗∗ , − i + = (cid:16) − θ i + (cid:17) ˆ G HLLC , − i + + θ i + ˆ G ∗∗ , − i + ; endend computed as: L error = (cid:118)(cid:117)(cid:117)(cid:116) N − (cid:88) i =0 N − (cid:88) j =0 ∆ x ∆ y (cid:0) α i,j − α i,j (cid:1) / N − (cid:88) i =0 N − (cid:88) j =0 ∆ x ∆ y, (108)where α i,j is the exact solution of volume fraction at the corresponding grid point. It27 lgorithm 6: Flux limiting procedure for positive ρ ˜ c . for all midpoints doif ρ ˜ c (cid:16) W HLLC , + i (cid:17) < (cid:15) ρ ˜ c then θ + i + = 0; else if ρ ˜ c (cid:0) W ∗∗ , + i (cid:1) < (cid:15) ρ ˜ c then Solve θ + i + from the formula: (cid:16) − θ + i + (cid:17) ρ ˜ c (cid:16) W HLLC , + i (cid:17) + θ + i + ρ ˜ c (cid:0) W ∗∗ , + i (cid:1) = (cid:15) ρ ˜ c ; else θ + i + = 1; endif ρ ˜ c (cid:16) W HLLC , − i +1 (cid:17) < (cid:15) ρ ˜ c then θ − i + = 0; else if ρ ˜ c (cid:0) W ∗∗ , − i +1 (cid:1) < (cid:15) ρ ˜ c then Solve θ − i + from the formula: (cid:16) − θ − i + (cid:17) ρ ˜ c (cid:16) W HLLC , − i +1 (cid:17) + θ − i + ρ ˜ c (cid:0) W ∗∗ , − i +1 (cid:1) = (cid:15) ρ ˜ c ; else θ − i + = 1; end Set θ i + = min (cid:16) θ + i + , θ − i + (cid:17) ;Perform convex combination (applying Jensen’s inequality):ˆ G ∗∗∗ , + i + = (cid:16) − θ i + (cid:17) ˆ G HLLC , + i + + θ i + ˆ G ∗∗ , + i + ;ˆ G ∗∗∗ , − i + = (cid:16) − θ i + (cid:17) ˆ G HLLC , − i + + θ i + ˆ G ∗∗ , − i + ; end Fluid Phase number γ p ∞ (Pa)Liquid water 1 6 .
12 3 . .
40 0
Table 1: Properties of the fluids. can be seen from the table that all schemes achieve the expected rates of convergence.28 (kg m − ) ρ (kg m − ) u (m s − ) v (m s − ) p (Pa) α . .
25 sin [ π ( x + y )] Table 2: Initial conditions of 2D convergence problem. ρ (kg m − ) ρ (kg m − ) u (m s − ) v (m s − ) p (Pa) α . .
25 sin [ π ( x + y − t )] Table 3: Exact solutions of 2D convergence problem.
Number of HLLC PP-WCNS-ISgrid points error order error order8 . −
04 3 . − . −
04 0.98 1 . −
07 7.6932 . −
04 0.99 4 . −
09 4.9764 . −
05 1.00 1 . −
10 4.85128 . −
05 1.00 5 . −
12 4.96256 . −
05 1.00 1 . −
13 4.97
Table 4: L errors and orders of convergence of volume fraction for the 2D advection problem fromdifferent schemes at t = 0 . The next multi-phase problem is a 1D problem with the advection of two materialinterfaces. The settings of this problem are similar to those in [8, 49, 2]. The initialconditions are given by table 5. Periodic conditions are applied at both boundaries. Thespatial domain is x ∈ [0 ,
1) m and the final time is at t = 0 .
01 s. Simulations are evolvedwith constant time steps ∆ t = 1 . − x = 0 .
005 m. The two material interfaces have exactly advected one period at the endof the simulations. α ρ (kg m − ) α ρ (kg m − ) u (m s − ) p (Pa) α . ≤ x < .
75 1000 1 . − − . − . − . − Table 5: Initial conditions of 1D material interface advection problem.
The density fields obtained with the two schemes at the final simulation time are com-pared with the exact solution in figure 5. It can be seen that the high-order PP-WCNS-IS29an capture the material interfaces with much smaller numerical widths compared to thefirst order HLLC scheme and no spurious oscillations are observed at the two materialinterfaces for PP-WCNS-IS. Both the velocity and pressure fields are uniform and con-stant in this advection problem. In figure 6, it can be seen that the relative errors invelocity and pressure fields for both schemes are insignificantly small and the uniformand constant fields are maintained well over time. . . . . . . x (m) ρ ( k g m − ) (a) Global density profile .
60 0 .
65 0 .
70 0 .
75 0 .
80 0 .
85 0 . x (m) ρ ( k g m − ) (b) Local density profileFigure 5: Material interface advection problem at t = 0 .
01 s using different schemes. Black solid line:exact; red circles: HLLC; blue squares: PP-WCNS-IS. . . . . . . x (m) − . − . − . − . . . . . . r e l a t i v ee rr o r i n u × − (a) Velocity . . . . . . x (m) − − r e l a t i v ee rr o r i n p × − (b) PressureFigure 6: Relative errors for the material interface advection problem at t = 0 .
01 s using differentschemes. Red circles: HLLC; blue squares: PP-WCNS-IS.
This gas/liquid shock tube problem is taken from Chen and Liang [6] and Wang et al.[48]. The initial conditions are given by table 6. Extrapolations are applied at bothboundaries. The spatial domain is x ∈ [0 , .
5] m and the final time is at t = 3e − t = 1 . − ρ (kg m − ) α ρ (kg m − ) u (m s − ) p (Pa) α x < . . − . − . − x ≥ . . − . . − Table 6: Initial conditions of 1D gas/liquid Sod shock tube problem.
Figures 7 and 8 compare the numerical solutions from the two schemes with thereference solutions. In figure 7(a), it can be seen that PP-WCNS-IS can capture both thematerial interface (the left larger density jump) and the shock (the right smaller densityjump) with only a few grid points. On the other hand, the first order HLLC schemeis too dissipative that both the material interface and shock are smeared out severely,and hence the two discontinuities cannot be distinguished at this grid resolution. Infigures 7(b) and 8(a), it can be seen that PP-WCNS-IS can give more accurate solutionsin velocity and pressure fields compared to first order HLLC scheme. However, as seenin figure 8(b), PP-WCNS-IS produces a slightly larger undershoot and overshoot at theexpansion fan. . . . . . . x (m) ρ ( k g m − ) (a) Global density profile . . . . . . x (m) u ( m s − ) (b) Global velocity profileFigure 7: Gas/liquid Sod shock tube problem at t = 3e − This is a multi-material version modified from the well-known single-phase 1D planarSedov blast wave problem [37, 55, 20]. Initially there is a singularity of highly pressurizedair at the center of the domain filled with very low pressure water. Blast waves are createdat the original position of the singularity and propagate towards the domain boundaries.The initial conditions are given by table 7. Extrapolations are applied at both boundaries.The spatial domain is x ∈ [0 , .
0] m and the final time is at t = 1e − t = 2 . − . . . . . . x (m) . . . . . p ( P a ) × (a) Global pressure profile .
20 0 .
25 0 .
30 0 .
35 0 .
40 0 .
45 0 .
50 0 .
55 0 . x (m) − − p ( P a ) × (b) Local pressure profileFigure 8: Gas/liquid Sod shock tube problem at t = 3e − α ρ (kg m − ) α ρ (kg m − ) u (m s − ) p (Pa) α x < − . x or x > . x . − − . − . − . − . / ∆ x . − Table 7: Initial conditions of 1D planar multi-material Sedov blast wave problem. that both schemes can capture the blast waves without spurious oscillations. However,the velocity and pressure profiles computed with PP-WCNS-IS are much sharper at theshock fronts while the shocks captured with first order HLLC are severely smeared outdue to excess numerical dissipation. . . . . . . . . . x (m) − − ρ ( k g m − ) (a) Global density profile . . . . . . . . . x (m) − − − u ( m s − ) (b) Global velocity profileFigure 9: 1D multi-material Sedov problem at t = 1e − . . . . . . . . . x (m) − − p ( P a ) × (a) Global pressure profile . . . . . . . . . x (m) . . . . . . α (b) Global water volume fraction profileFigure 10: 1D multi-material Sedov problem at t = 1e − The case of a Mach 2.4 planar shock interacting with a water cylinder in the paperby Sembian et al. [38] is simulated. The purpose of this test case is to investigate thereliability of the flow model with the high-order diffuse interface method for simulatingtwo-phase flows with shocks. Figure 11 shows the schematic of the initial flow fieldand domain. The water cylinder is initially placed at location [4 cm , × α ρ (kg m − ) α ρ (kg m − ) u (m s − ) v (m s − ) p (Pa) α pre-shock air 1 . − . . − . − . . − . − . − . − Table 8: Initial conditions of 2D Mach 2.4 shock water cylinder interaction problem [38].
The grey schlieren images from the experiment [38] are shown in the left column offigure 12. At the instance when the incident shock interacts with the water column, ashock is reflected upstream since the acoustic impedance of water is higher than that ofair. The reflected shock interacts with the incident shock to generate a triple point wherethe reflected shock, incident shock and a Mach stem along with its slip line conincide.Meanwhile, there is also a shock transmitted into the water column. The transmittedshock travels faster than the shocks outside the water column and it gets reflected asan expansion wave when the transmitted shock reaches the downstream water-air inter-face. The reflected expansion wave focuses at a point due to the column’s downstreamconcave geometry, where negative presure is produced due to tensile stresses. The re-flected expansion wave forms a “horse-shoe” structure after focusing and is reflected33gain at the upstream water-air interface. The expansion wave continues to get reflectedinside the water column repeatedly. As the water column is a buff body, the surround-ing air separates in the adverse pressure gradient region on the water column surface.Therefore, recirculation regions are created and two counter-rotating vortices are formeddownstream of the flow. In the right column of figure 12, the density gradients computedwith the simulation results using the high-order PP-WCNS-IS scheme at different timesare displayed. Compared with the schlieren images from the experiment, it can be seenthat most of the wave features, such as the incident, reflected and transmitted shocks,and expansion waves are captured accurately. Also, the two counter-rotating vortices arereproduced in the simulation.
Post-shock air Pre-shock airWater11 . . . xy Figure 11: Schematic diagram of 2D Mach 2.4 shock water cylinder interaction problem by Sembianet al. [38].
This is the more extreme case of the previous test case. In this problem, a Mach 10shock wave in air interacts with a water cylinder with diameter of 8 L in a domain of[0 , L ] × [ − L, L ], where L = 1 mm is chosen. Figure 13 shows the schematic of theinitial flow field and domain. The water cylinder is initially placed at location [13 L, ×
768 mesh.This is a very extreme problem due to the high incident shock Mach number initially.The problem is simulated with the first order HLLC and the high-order PP-WCNS-IS.While the first order scheme has no numerical difficulty in this test case, numerical failureis experienced with WCNS-IS without the positivity- and boundedness-preserving lim-iters as the speed of sound becomes imaginary. This can happen at the strong incidentshock, or at the low pressure regions created by the expansion waves inside the watercolumn and at the counter-rotating vortices. The positivity- and boundedness-preservinglimiters are necessary for the WCNS-IS scheme in this problem. Figure 14 shows the com-parison of numerical schlieren defined as exp ( |∇ ρ | / |∇ ρ | max ) between the two schemes.34 ρ (kg m − ) α ρ (kg m − ) u (m s − ) v (m s − ) p (Pa) α pre-shock air 1 . − . . − . − . . . − . − . − . − Table 9: Initial conditions of 2D Mach 10 shock water cylinder interaction problem.
From the figure, we can see that the PP-WCNS-IS has much thinner interface thicknessover time compared with first order HLLC scheme and this is consistent with other testproblems. Vortical features produced by the hydrodynamic instability due to baroclinictorque are observed at the interface for PP-WCNS-IS as time evolves. However, the firstorder scheme is too dissipative to produce the roll-up of the interfaces at the chosen gridresolution. The comparison of speed of sound between the two schemes can be seen infigure 15. Finally, the volume fraction fields of both schemes are shown in figure 16. Itshould be noted that the volume fraction field is also verified to be always bounded by thecorresponding threshold chosen in the positivity- and boundedness-preserving limiters forPP-WCNS-IS.
This test case is a multi-phase version of the popular Mach 2000 jet problem firstproposed by Zhang and Shu [53]. In this problem, a Mach 100 water jet enters a domainfull of ambient air. The domain size is [0 , L ] × [ − . L, . L ], where L = 1 m is chosen.The initial conditions of the ambient air are given by table 10. Constant extrapolationis used at top, bottom and right boundaries. The left boundary is described by Dirichletboundary conditions given by table 11. The speed of the jet is 1 . − , which isaround Mach 100 with respect to the sound speed in the water jet. The computationsare performed on a 1024 ×
512 mesh. α ρ (kg m − ) α ρ (kg m − ) u (m s − ) v (m s − ) p (Pa) α . − . . − Table 10: Initial conditions of 2D Mach 100 water jet problem.
The comparison of speed of sound between the two schemes is shown in figure 17.Despite the large jump in sound speed across the bow shock ahead of the high speedwater jet, none of the schemes fail due to the positivity-preserving properties of bothschemes for sound speed. The main difference between the two schemes is at the waterjet front where the interface at the water jet front produced by first order HLLC schemeis heavily smeared out while that of PP-WCNS-IS is reasonably captured with only a fewgrid points. There are also some small but obvious numerical artificts at the bow shock inthe solutions computed with the first order scheme. In figure 18, the numerical schilren35 ρ (kg m − ) α ρ (kg m − ) u (m s − ) v (m s − ) p (Pa) α y ∈ [ − . L, . L ] 1000 1 . − . . − . − . − . . − Table 11: Left boundary conditions of 2D Mach 100 water jet problem. between the two schemes are compared. Since the interface water jet front is seriouslydiffused, the shape of the water-air interface at the water front cannot be visualizedat all. Figure 19 compares the volume fraction field of water at different times. Thevolume fraction field is verified to be bounded in PP-WCNS-IS. Similar to other fields,volume fraction interface at the water jet front is very diffused for the first order schemecompared to the high-order PP-WCNS-IS.
7. Concluding remarks
In this work, limiting procedures were proposed on a high-order finite differencescheme that can preserve the positivity of partial density of each phase and squared speedof sound and also the boundedness of the volume fractions in 1D and multi-dimensionalgas-liquid two-phase problems under a mild assumption on the material properties of thegas and liquid. The procedures consist of two stages which limit the WENO interpola-tion and flux reconstruction respectively in the high-order WCNS-IS algorithm. Discreteconservation of solutions for the conservation equations in the five-equation model is stillmaintained even with the limiting. The overall positivity- and boundedness-preservingscheme, PP-WCNS-IS, was tested with different severe problems, under suitable CFLconditions. Comparison between the results of the first order HLLC scheme and PP-WCNS-IS showed the low dissipation and high resolution properties of the latter schemewhile its robustness is also ensured. The positivity- and boundedness-preserving limitingprocedures can also be potentially used with any conservative finite difference and finitevolume schemes for gas-liquid two-phase flows. Future work includes generalization ofthe positivity- and boundedness-preserving limiters to more general equation of statesand multi-phase flows with more than two species. There will also be future investiga-tions on the use of the diffuse interface method for simulations of space vehicle launcheswith water-based sound suppression systems [47, 26].
8. Acknowledgments
This work was partially supported by the NASA Exploration Ground Systems (EGS)program and the NASA Engineering and Safety Center (NESC). Computer time hasbeen provided by the NASA Advanced Supercomputing (NAS) facility at NASA AmesResearch Center. We also gratefully acknowledge Dr. Bruce T. Vu, Dr. Jeffrey A.Housman and Dr. Oliver M. Browne for valuable discussions.36 ppendix A. Characteristic decomposition
The choice of variables for WENO reconstruction and interpolation is very criti-cal to avoid spurious oscillations across discontinuities, especially across the materialinterfaces. It was shown in [23] that if conservative variables are chosen for WENO re-construction, spurious oscillations will appear at material interfaces. Primitive variableswere suggested [23, 8, 49] for reconstruction and interpolation in order to maintain pres-sure and velocity equilibria across interfaces. Furthermore, WENO reconstruction andinterpolation of characteristic variables projected from primitive variables can avoid theinteraction of discontinuities in different characteristic fields. To illustrate how the prim-itive variables are converted into characteristic variables, we follow previous works [8, 49]by first rewriting the 2D governing equations in the quasi-linear primitive form: ∂ V ∂t + A ( V ) ∂ V ∂x + B ( V ) ∂ V ∂y = 0 , (A.1)where V is the vector of primitive variables. V and matrix A are given by: V = α ρ α ρ uvpα , A = u α ρ u α ρ u ρ
00 0 0 u ρc u
00 0 0 0 0 u , B = v α ρ v α ρ v v ρ
00 0 0 ρc v
00 0 0 0 0 v . (A.2)The eigenvectors of the matrices A and B have to be determined first in order to trans-form primitive variables to characteristic variables. The eigenvalue decompositions ofthe matrices are given by: A = R A Λ A R − A , B = R B Λ B R − B , (A.3)37here R A , R − A and Λ A are given by: R A = − α ρ c α ρ c − α ρ c α ρ c − ρc ρc , R − A = − ρc
01 0 0 0 − α ρ ρc
00 1 0 0 − α ρ ρc
00 0 0 1 0 00 0 0 0 0 10 0 1 0 ρc , Λ A = u − c u u u u
00 0 0 0 0 u + c . (A.4) R B , R − B and Λ B have similar corresponding forms. Appendix B. Convergence analysis of the incremental-stencil WCNS
The convergence analysis discussed in this section is the extension of the convergenceanalysis of Yan et al. [50] to HCS with WENO-IS interpolation. Assume that we have a1D scalar hyperbolic conservation law of dependent variable u : ∂ t u + ∂ x F = 0 , (B.1)where F = F ( u ) is the flux. Under the assumption that ∂F ( u ) /∂u > F i + = F (cid:16) ˜ u i + ,L (cid:17) , where ˜ u i + ,L is the left-biased WENO interpolatedvalue of u . The subscript “ L ” is dropped in this section for convenience.After Taylor series expansion of the interpolation equations (84)-(87) for the sub-stencils, we obtain:˜ u ki + = (cid:40) u i + + A ki ∆ x + O (cid:0) ∆ x (cid:1) , if k = 0 , ,u i + + A ki ∆ x + O (cid:0) ∆ x (cid:1) , if k = 2 , , (B.2)where A i = u (cid:48)(cid:48) i / A i = − u (cid:48)(cid:48) i / A i = − u (cid:48)(cid:48)(cid:48) i /
16, and A i = − u (cid:48)(cid:48)(cid:48) i /
16. If we replace thenonlinear weights with the corresponding linear weights in equation (91), (cid:88) k =0 d k ˜ u ki + = u i + + B i ∆ x + O (cid:0) ∆ x (cid:1) , (B.3)where B i = − u (5) i / u i + = (cid:88) k =0 ω k ˜ u ( k ) i + = (cid:88) k =0 d k ˜ u ( k ) i + + (cid:88) k =0 ( ω k − d k ) ˜ u ( k ) i + . (B.4)38sing equations (B.2) and (B.3), the equation above becomes:˜ u i + = u i + + B i ∆ x + O (cid:0) ∆ x (cid:1) + (cid:88) k =0 (cid:110) ( ω k − d k ) (cid:104) u i + + A ki ∆ x + O (cid:0) ∆ x (cid:1)(cid:105)(cid:111) + (cid:88) k =2 (cid:110) ( ω k − d k ) (cid:104) u i + + A ki ∆ x + O (cid:0) ∆ x (cid:1)(cid:105)(cid:111) = u i + + B i ∆ x + O (cid:0) ∆ x (cid:1) + u i + (cid:88) k =0 ( ω k − d k )+ ∆ x (cid:88) k =0 (cid:2) A ki ( ω k − d k ) (cid:3) + ∆ x (cid:88) k =2 (cid:2) A ki ( ω k − d k ) (cid:3) + (cid:88) k =0 (cid:2) ( ω k − d k ) O (cid:0) ∆ x (cid:1)(cid:3) + (cid:88) k =2 (cid:2) ( ω k − d k ) O (cid:0) ∆ x (cid:1)(cid:3) . (B.5)The derivation is similar for ˜ u i − . Therefore,˜ u i ± − u i ± = B i ∆ x + u i ± (cid:88) k =0 (cid:0) ω ± k − d k (cid:1) + ∆ x (cid:88) k =0 (cid:2) A ki (cid:0) ω ± k − d k (cid:1)(cid:3) + ∆ x (cid:88) k =2 (cid:2) A ki (cid:0) ω ± k − d k (cid:1)(cid:3) + (cid:88) k =0 (cid:2)(cid:0) ω ± k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) + (cid:88) k =2 (cid:2)(cid:0) ω ± k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) + O (cid:0) ∆ x (cid:1) , (B.6)at midpoints i +1 / i − /
2. The superscript ± is added to ω k to distinguish the valuesat the two different midpoints. Note that no ± is added to d k , A ki and B i since they havethe same values at the two midpoints. Using the fact that (cid:80) k =0 ω ± k = (cid:80) k =0 d k = 1, wefinally get:˜ u i ± − u i ± = B i ∆ x + ∆ x (cid:88) k =0 (cid:2) A ki (cid:0) ω ± k − d k (cid:1)(cid:3) + ∆ x (cid:88) k =2 (cid:2) A ki (cid:0) ω ± k − d k (cid:1)(cid:3) + (cid:88) k =0 (cid:2)(cid:0) ω ± k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) + (cid:88) k =2 (cid:2)(cid:0) ω ± k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) + O (cid:0) ∆ x (cid:1) . (B.7)We will show that the WCNS-IS is fifth order accurate if ω ± k − d k = O (cid:0) ∆ x (cid:1) . If we39ssume that ω ± k − d k = O (cid:0) ∆ x (cid:1) , by Taylor series expansion,˜ F i ± = F i ± + ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ± (cid:16) ˜ u i ± − u i ± (cid:17) + 12 ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ± (cid:16) ˜ u i ± − u i ± (cid:17) + · · · = F i ± + ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ± (cid:16) ˜ u i ± − u i ± (cid:17) + O (cid:0) ∆ x (cid:1) . (B.8)Substituting equations (B.7) and (B.8) into the 1D version of equation (77) (noticing F isthe scalar version of G x here) and using Taylor-sereis expanded equation (78) (assuming ψ = 256 / (cid:100) ∂F∂x (cid:12)(cid:12)(cid:12)(cid:12) i = ∂F∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1) + ψ (cid:40) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:88) k =0 (cid:2)(cid:0) ω + k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:88) k =0 (cid:2)(cid:0) ω − k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3)(cid:41) + ψ (cid:40) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:88) k =2 (cid:2)(cid:0) ω + k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:88) k =2 (cid:2)(cid:0) ω − k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3)(cid:41) + ψ ∆ x B i (cid:32) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:33) + ψ ∆ x (cid:40) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:88) k =0 (cid:2) A ki (cid:0) ω + k − d k (cid:1)(cid:3) − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:88) k =0 (cid:2) A ki (cid:0) ω − k − d k (cid:1)(cid:3)(cid:41) + ψ ∆ x (cid:40) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:88) k =2 (cid:2) A ki (cid:0) ω + k − d k (cid:1)(cid:3) − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:88) k =2 (cid:2) A ki (cid:0) ω − k − d k (cid:1)(cid:3)(cid:41) + O (cid:0) ∆ x (cid:1) . (B.9)Also, by Taylor series expansion, ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + = ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:16) u i + − u i (cid:17) + 12 ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:16) u i + − u i (cid:17) + · · · , (B.10) u i + = u i + ∆ x ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1) . (B.11)40herefore, ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + = ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:20) ∆ x ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21) + 12 ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:20) ∆ x ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21) + O (cid:0) ∆ x (cid:1) = ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:18) ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i (cid:19) + O (cid:0) ∆ x (cid:1) . (B.12)Similarly, by Taylor series expansion, ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − = ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:20) − ∆ x ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21) + 12 ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:20) − ∆ x ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21) + O (cid:0) ∆ x (cid:1) = ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i (cid:18) ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i (cid:19) + O (cid:0) ∆ x (cid:1) . (B.13)41s a result, equation (B.9) is simplified to: (cid:100) ∂F∂x (cid:12)(cid:12)(cid:12)(cid:12) i = ∂F∂x (cid:12)(cid:12)(cid:12)(cid:12) i + ψ (cid:40) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:88) k =0 (cid:2)(cid:0) ω + k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:88) k =0 (cid:2)(cid:0) ω − k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3)(cid:41) + ψ (cid:40) ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:88) k =2 (cid:2)(cid:0) ω + k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3) − ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i − (cid:88) k =2 (cid:2)(cid:0) ω − k − d k (cid:1) O (cid:0) ∆ x (cid:1)(cid:3)(cid:41) + ψ ∆ x B i ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + (cid:20) ψ ∆ x ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21) (cid:88) k =0 (cid:2) A ki (cid:0) ω + k − ω − k (cid:1)(cid:3) + (cid:20) ψ ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21)(cid:40) (cid:88) k =0 (cid:2) A ki (cid:0) ω + k − d k (cid:1)(cid:3) + (cid:88) k =0 (cid:2) A ki (cid:0) ω − k − d k (cid:1)(cid:3)(cid:41) + (cid:20) ψ ∆ x ∂F∂u (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21) (cid:88) k =2 (cid:2) A ki (cid:0) ω + k − ω − k (cid:1)(cid:3) + (cid:20) ψ ∆ x ∂ F∂u (cid:12)(cid:12)(cid:12)(cid:12) i ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) i + O (cid:0) ∆ x (cid:1)(cid:21)(cid:40) (cid:88) k =2 (cid:2) A ki (cid:0) ω + k − d k (cid:1)(cid:3) + (cid:88) k =2 (cid:2) A ki (cid:0) ω − k − d k (cid:1)(cid:3)(cid:41) + O (cid:0) ∆ x (cid:1) . (B.14)It can be seen from equation (B.14) that a sufficient condition for fifth order convergenceis ω ± k − d k = O (cid:0) ∆ x (cid:1) , which is already assumed earlier.To prove that ω ± k − d k = O (cid:0) ∆ x (cid:1) is true, we can perform Taylor series expansion on42he smoothness indicators given by equations (95)–(99), β = ( u (cid:48) i ) ∆ x + u (cid:48) i u (cid:48)(cid:48) i ∆ x + (cid:32) u (cid:48) i u (cid:48)(cid:48)(cid:48) i u (cid:48)(cid:48) i ) (cid:33) ∆ x + (cid:32) u (cid:48) i u (4) i
12 + u (cid:48)(cid:48) i u (cid:48)(cid:48)(cid:48) i (cid:33) ∆ x + O (cid:0) ∆ x (cid:1) , (B.15) β = ( u (cid:48) i ) ∆ x − u (cid:48) i u (cid:48)(cid:48) i ∆ x + (cid:32) u (cid:48) i u (cid:48)(cid:48)(cid:48) i u (cid:48)(cid:48) i ) (cid:33) ∆ x + (cid:32) − u (cid:48) i u (4) i − u (cid:48)(cid:48) i u (cid:48)(cid:48)(cid:48) i (cid:33) ∆ x + O (cid:0) ∆ x (cid:1) , (B.16) β = ( u (cid:48) i ) ∆ x + (cid:32) u (cid:48) i u (cid:48)(cid:48)(cid:48) i u (cid:48)(cid:48) i ) (cid:33) ∆ x + O (cid:0) ∆ x (cid:1) , (B.17) β = ( u (cid:48) i ) ∆ x + (cid:32) − u (cid:48) i u (cid:48)(cid:48)(cid:48) i u (cid:48)(cid:48) i ) (cid:33) ∆ x + (cid:32) − u (cid:48) i u (4) i u (cid:48)(cid:48) i u (cid:48)(cid:48)(cid:48) i (cid:33) ∆ x + O (cid:0) ∆ x (cid:1) , (B.18) β = ( u (cid:48) i ) ∆ x + (cid:32) − u (cid:48) i u (cid:48)(cid:48)(cid:48) i u (cid:48)(cid:48) i ) (cid:33) ∆ x + (cid:32) u (cid:48) i u (4) i − u (cid:48)(cid:48) i u (cid:48)(cid:48)(cid:48) i (cid:33) ∆ x + O (cid:0) ∆ x (cid:1) . (B.19)The Taylor series expansion of reference smoothness indicator (equation (101)) gives: τ = ( u (cid:48)(cid:48)(cid:48) i ) ∆ x + O (cid:0) ∆ x (cid:1) . (B.20)Therefore, τ β k + (cid:15) = O (cid:0) ∆ x (cid:1) , (B.21)provided not at critical points. As explained in Borges et al. [4] and Wang et al. [48],this is sufficient to have ω ± k − d k = O (cid:0) ∆ x (cid:1) . References [1] G. Allaire, S. Clerc, and S. Kokh. A five-equation model for the simulation of interfaces betweencompressible fluids.
Journal of Computational Physics , 181(2):577–616, 2002.[2] M. Aslani and J. D. Regele. A localized artificial diffusivity method to simulate compressiblemultiphase flows using the stiffened gas equation of state.
International Journal for NumericalMethods in Fluids , 88(9):413–433, 2018.[3] P. Batten, N. Clarke, C. Lambert, and D. Causon. On the choice of wavespeeds for the HLLCRiemann solver.
SIAM Journal on Scientific Computing , 18(6):1553–1570, 1997.[4] R. Borges, M. Carmona, B. Costa, and W. S. Don. An improved weighted essentially non-oscillatoryscheme for hyperbolic conservation laws.
Journal of Computational Physics , 227(6):3191–3211,2008.[5] C.-H. Chang and M.-S. Liou. A robust and accurate approach to computing compressible multiphaseflow: Stratified flow model and AUSM+-up scheme.
Journal of Computational Physics , 225(1):840–873, 2007.[6] H. Chen and S. Liang. Flow visualization of shock/water column interactions.
Shock Waves , 17(5):309–321, 2008.
7] J. Cheng and C.-W. Shu. Positivity-preserving Lagrangian scheme for multi-material compressibleflow.
Journal of Computational Physics , 257:143–168, 2014.[8] V. Coralic and T. Colonius. Finite-volume WENO scheme for viscous compressible multicomponentflows.
Journal of computational physics , 274:95–121, 2014.[9] X. Deng. New high-order hybrid cell-edge and cell-node weighted compact nonlinear schemes. In , page 3857, 2011.[10] X. Deng and H. Zhang. Developing high-order weighted compact nonlinear schemes.
Journal ofComputational Physics , 165(1):22–44, 2000.[11] B. Einfeldt, C.-D. Munz, P. L. Roe, and B. Sj¨ogreen. On Godunov-type methods near low densities.
Journal of Computational Physics , 92(2):273–295, 1991.[12] L. Fu, X. Y. Hu, and N. A. Adams. A family of high-order targeted ENO schemes for compressible-fluid simulations.
Journal of Computational Physics , 305:333–359, 2016.[13] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretizationmethods.
SIAM review , 43(1):89–112, 2001.[14] S. Gottlieb, D. I. Ketcheson, and C.-W. Shu. High order strong stability preserving time discretiza-tions.
Journal of Scientific Computing , 38(3):251–289, 2009.[15] J. Gressier, P. Villedieu, and J.-M. Moschetta. Positivity of flux vector splitting schemes.
Journalof Computational Physics , 155(1):199–220, 1999.[16] F. H. Harlow and A. A. Amsden. Fluid dynamics. Monograph LA-4700, Los Alamos National Lab.,NM (United States), 1971.[17] J. A. Housman, C. C. Kiris, and M. M. Hafez. Time-derivative preconditioning methods for multi-component flows—part i: Riemann problems.
Journal of applied mechanics , 76(2), 2009.[18] J. A. Housman, C. C. Kiris, and M. M. Hafez. Time-derivative preconditioning methods for mul-ticomponent flows—part ii: Two-dimensional applications.
Journal of applied mechanics , 76(3),2009.[19] X. Hu, Q. Wang, and N. Adams. An adaptive central-upwind weighted essentially non-oscillatoryscheme.
Journal of Computational Physics , 229(23):8952–8965, 2010.[20] X. Y. Hu, N. A. Adams, and C.-W. Shu. Positivity-preserving method for high-order conservativeschemes solving compressible Euler equations.
Journal of Computational Physics , 242:169–180,2013.[21] S. S. Jain, A. Mani, and P. Moin. A conservative diffuse-interface method for compressible two-phaseflows.
Journal of Computational Physics , page 109606, 2020.[22] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted eno schemes.
Journal of compu-tational physics , 126(1):202–228, 1996.[23] E. Johnsen and T. Colonius. Implementation of WENO schemes in compressible multicomponentflow problems.
Journal of Computational Physics , 219(2):715–732, 2006.[24] E. Johnsen, J. Larsson, A. V. Bhagatwala, W. H. Cabot, P. Moin, B. J. Olson, P. S. Rawat, S. K.Shankar, B. Sj¨ogreen, H. C. Yee, X. Zhong, and S. K. Lele. Assessment of high-resolution methodsfor numerical simulations of compressible turbulence with shock waves.
Journal of ComputationalPhysics , 229(4):1213–1237, 2010.[25] S. Kawai, S. K. Shankar, and S. K. Lele. Assessment of localized artificial diffusivity scheme forlarge-eddy simulation of compressible turbulent flows.
Journal of Computational Physics , 229(5):1739–1762, 2010.[26] C. C. Kiris, J. A. Housman, M. F. Barad, C. Brehm, E. Sozer, and S. Moini-Yekta. Computationalframework for launch, ascent, and vehicle aerodynamics (LAVA).
Aerospace Science and Technology ,55:189–219, 2016.[27] T. Linde, P. Roe, T. Linde, and P. Roe. Robust Euler codes. In , page 2098, 1997.[28] M.-S. Liou. A sequel to AUSM: AUSM+.
Journal of computational Physics , 129(2):364–382, 1996.[29] T. Nonomura and K. Fujii. Effects of difference scheme type in high-order weighted compactnonlinear schemes.
Journal of Computational Physics , 228(10):3533–3539, 2009.[30] T. Nonomura and K. Fujii. Robust explicit formulation of weighted compact nonlinear scheme.
Computers & Fluids , 2013.[31] T. Nonomura, N. Iizuka, and K. Fujii. Increasing order of accuracy of weighted compact nonlinearscheme.
AIAA Paper , 893, 2007.[32] G. Perigaud and R. Saurel. A compressible flow model with capillary effects.
Journal of Computa-tional Physics , 209(1):139–178, 2005.[33] B. Perthame and C.-W. Shu. On positivity preserving finite volume schemes for euler equations.
Numerische Mathematik , 73(1):119–130, 1996.
34] R. Saurel and R. Abgrall. A simple method for compressible multifluid flows.
SIAM Journal onScientific Computing , 21(3):1115–1145, 1999.[35] R. Saurel and C. Pantano. Diffuse-interface capturing methods for compressible two-phase flows.
Annual Review of Fluid Mechanics , 50:105–130, 2018.[36] K. Sebastian and C.-W. Shu. Multidomain WENO finite difference method with interpolation atsubdomain interfaces.
Journal of Scientific Computing , 19(1-3):405–438, 2003.[37] L. I. Sedov.
Similarity and dimensional methods in mechanics . CRC press, 1993.[38] S. Sembian, M. Liverts, N. Tillmark, and N. Apazidis. Plane shock wave interaction with a cylin-drical water column.
Physics of Fluids , 28(5):056102, 2016.[39] H. Shen, C.-Y. Wen, M. Parsani, and C.-W. Shu. Maximum-principle-satisfying space-time con-servation element and solution element scheme applied to compressible multifluids.
Journal ofComputational Physics , 330:668–692, 2017.[40] C.-W. Shu. Total-variation-diminishing time discretizations.
SIAM Journal on Scientific andStatistical Computing , 9(6):1073–1084, 1988.[41] C.-W. Shu. High-order finite difference and finite volume weno schemes and discontinuous galerkinmethods for cfd.
International Journal of Computational Fluid Dynamics , 17(2):107–118, 2003.[42] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturingschemes.
Journal of computational physics , 77(2):439–471, 1988.[43] K.-M. Shyue. A fluid-mixture type algorithm for compressible multicomponent flow with van derWaals equation of state.
Journal of Computational Physics , 156(1):43–88, 1999.[44] A. Subramaniam, M. L. Wong, and S. K. Lele. A high-order weighted compact high resolutionscheme with boundary closures for compressible turbulent flows with shocks.
Journal of Computa-tional Physics , 397:108822, 2019.[45] H. Takahira, T. Matsuno, and K. Shuto. Numerical investigations of shock–bubble interactions inmercury.
Fluid Dynamics Research , 40(7-8):510, 2008.[46] V. A. Titarev and E. F. Toro. Finite-volume WENO schemes for three-dimensional conservationlaws.
Journal of Computational Physics , 201(1):238–260, 2004.[47] B. T. Vu, N. Bachchan, O. Peroomian, and V. Akdag. Multiphase modeling of water injection onflame deflector. In , page 2592, 2013.[48] B. Wang, G. Xiang, and X. Y. Hu. An incremental-stencil WENO reconstruction for simulation ofcompressible two-phase flows.
International Journal of Multiphase Flow , 104:20–31, 2018.[49] M. L. Wong and S. K. Lele. High-order localized dissipation weighted compact nonlinear schemefor shock-and interface-capturing in compressible flows.
Journal of Computational Physics , 339:179–209, 2017.[50] Z. Yan, H. Liu, M. Mao, H. Zhu, and X. Deng. New nonlinear weights for improving accuracy andresolution of weighted compact nonlinear scheme.
Computers & Fluids , 127:226–240, 2016.[51] S. Zhang, S. Jiang, and C.-W. Shu. Development of nonlinear weighted compact schemes withincreasingly higher order accuracy.
Journal of Computational Physics , 227(15):7294–7321, 2008.[52] X. Zhang and C.-W. Shu. On maximum-principle-satisfying high order schemes for scalar conser-vation laws.
Journal of Computational Physics , 229(9):3091–3120, 2010.[53] X. Zhang and C.-W. Shu. On positivity-preserving high order discontinuous Galerkin schemes forcompressible Euler equations on rectangular meshes.
Journal of Computational Physics , 229(23):8918–8934, 2010.[54] X. Zhang and C.-W. Shu. Positivity-preserving high order discontinuous Galerkin schemes forcompressible Euler equations with source terms.
Journal of Computational Physics , 230(4):1238–1248, 2011.[55] X. Zhang and C.-W. Shu. Positivity-preserving high order finite difference WENO schemes forcompressible Euler equations.
Journal of Computational Physics , 231(5):2245–2258, 2012. a) Experiment (b) t = 4 µ s, PP-WCNS-IS(c) Experiment (d) t = 17 µ s, PP-WCNS-IS(e) Experiment (f) t = 40 µ s, PP-WCNS-IS(g) Experiment (h) t = 67 µ s, PP-WCNS-ISFigure 12: Comparision of 2D Mach 2.4 shock water cylinder interaction problem by Sembian et al. [38].Left: shadowgraph images from [38]; right: density gradient from simulation with PP-WCNS-IS. ost-shock air Pre-shock airWater30 L L Lxy
Figure 13: Schematic diagram of 2D Mach 10 shock water cylinder interaction problem. a) t = 1 µ s, HLLC (b) t = 1 µ s, PP-WCNS-IS(c) t = 4 µ s, HLLC (d) t = 4 µ s, PP-WCNS-IS(e) t = 8 µ s, HLLC (f) t = 8 µ s, PP-WCNS-IS(g) t = 16 µ s, HLLC (h) t = 16 µ s, PP-WCNS-ISFigure 14: Numerical schlieren (exp (cid:0) |∇ ρ | / |∇ ρ | max (cid:1) ) of 2D Mach 10 shock water cylinder interactionproblem. a) t = 1 µ s, HLLC (b) t = 1 µ s, PP-WCNS-IS(c) t = 4 µ s, HLLC (d) t = 4 µ s, PP-WCNS-IS(e) t = 8 µ s, HLLC (f) t = 8 µ s, PP-WCNS-IS(g) t = 16 µ s, HLLC (h) t = 16 µ s, PP-WCNS-ISFigure 15: Speed of sound of 2D Mach 10 shock water cylinder interaction problem. a) t = 1 µ s, HLLC (b) t = 1 µ s, PP-WCNS-IS(c) t = 4 µ s, HLLC (d) t = 4 µ s, PP-WCNS-IS(e) t = 8 µ s, HLLC (f) t = 8 µ s, PP-WCNS-IS(g) t = 16 µ s, HLLC (h) t = 16 µ s, PP-WCNS-ISFigure 16: Volume fraction of water of 2D Mach 10 shock water cylinder interaction problem. a) t = 1 µ s, HLLC (b) t = 1 µ s, PP-WCNS-IS(c) t = 2 µ s, HLLC (d) t = 2 µ s, PP-WCNS-IS(e) t = 4 µ s, HLLC (f) t = 4 µ s, PP-WCNS-IS(g) t = 6 µ s, HLLC (h) t = 6 µ s, PP-WCNS-ISFigure 17: Speed of sound of 2D Mach 100 water jet problem. a) t = 1 µ s, HLLC (b) t = 1 µ s, PP-WCNS-IS(c) t = 2 µ s, HLLC (d) t = 2 µ s, PP-WCNS-IS(e) t = 4 µ s, HLLC (f) t = 4 µ s, PP-WCNS-IS(g) t = 6 µ s, HLLC (h) t = 6 µ s, PP-WCNS-ISFigure 18: Numerical schlieren (exp (cid:0) |∇ ρ | / |∇ ρ | max (cid:1) ) of 2D Mach 100 water jet problem. a) t = 1 µ s, HLLC (b) t = 1 µ s, PP-WCNS-IS(c) t = 2 µ s, HLLC (d) t = 2 µ s, PP-WCNS-IS(e) t = 4 µ s, HLLC (f) t = 4 µ s, PP-WCNS-IS(g) t = 6 µ s, HLLC (h) t = 6 µ s, PP-WCNS-ISFigure 19: Volume fraction of water of 2D Mach 100 water jet problem.s, PP-WCNS-ISFigure 19: Volume fraction of water of 2D Mach 100 water jet problem.