A quasi-conservative discontinuous Galerkin method for multi-component flows using the non-oscillatory kinetic flux
AA quasi-conservative discontinuous Galerkin method formulti-component flows using the non-oscillatory kinetic flux
Dongmi Luo , Jianxian Qiu , Jun Zhu , Yibing Chen Abstract
In this paper, a high order quasi-conservative discontinuous Galerkin (DG) methodusing the non-oscillatory kinetic flux is proposed for the 5-equation model of compress-ible multi-component flows with Mie-Gr¨uneisen equation of state. The method mainlyconsists of three steps: firstly, the DG method with the non-oscillatory kinetic flux isused to solve the conservative equations of the model; secondly, inspired by Abgrall’sidea, we derive a DG scheme for the volume fraction equation which can avoid theunphysical oscillations near the material interfaces; finally, a multi-resolution WENOlimiter and a maximum-principle-satisfying limiter are employed to ensure oscillation-free near the discontinuities, and preserve the physical bounds for the volume fraction,respectively. Numerical tests show that the method can achieve high order for smoothsolutions and keep non-oscillatory at discontinuities. Moreover, the velocity and pres-sure are oscillation-free at the interface and the volume fraction can stay in the interval[0,1].
Keywords : DG method, multi-component flows, non-oscillatory kinetic, Mie-Gr¨uneisenequations of state Institute of Applied Physics and Computational Mathematics, Beijing 100088, China. E-mail: [email protected]. School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modelingand High-Performance Scientific Computing, Xiamen University, Xiamen, Fujian 361005, China. E-mail:[email protected]. College of Sciences, Nanjing University of Aeronautics and Astronautics, Nanjing, Jiangsu 210016, China.E-mail: [email protected]. Institute of Applied Physics and Computational Mathematics, Beijing 100088, China, E-mail:chen [email protected]. a r X i v : . [ phy s i c s . c o m p - ph ] J u l Introduction
Numerical simulation of compressible multi-component flows with immiscible interfaceshas been an active research topic in the computational fluid dynamics because of theirapplication to a wide range of field, such as inertial confinement fusion, underwater bubbledynamic and so on. The major difficulty of the simulations of multi-component flows is howto track the material interfaces clearly.The numerical approaches can be split into two major groups with respect to the treat-ment of the material interface in the Eulerian framework. One is the sharp interfacemethod (SIM) and the other is the diffuse interface method (DIM). Sharp interface methods[4, 13, 25–27, 32, 34, 35, 47, 48, 54] view the multi-material interfaces as genuine disconti-nuities. Thus the sharp interfaces are strictly maintained. However, none of these methodsis able to dynamically create interfaces and to solve interfaces separating pure medium andmixtures as stated in [41].In contrast, in the diffuse interface approach the interfaces are viewed as a numericallydiffused zone, and an artificial mixture zone is created. A number of different models [1–3, 18, 40, 43–45] have been developed so far based on this idea, including 4-equation model,5-equation model, 7-equation model and so on. However, these models are usually non-conservative, which leads both theoretical and computational challenging problems [2]. Thespecial strategies are required to handle these non-conservative terms in order to keep thepressure and velocity non-oscillatory at the interfaces. The quasi-conservative approachdeveloped by Abgrall in [1] is an effective means to deal with this problem. Based onGodunov method for solution evolution, Shyue extended the idea to the different equationsof state (EOSs), such as stiffened gas equation [43], Van der Waals [45], Mie-Gr¨uneisen [44].Besides the traditional Godunov method, an alternative is the gas kinetic scheme (GKS) [49],which provides more physical information of the flow and is free from constructing Riemannsolver. In the past decades, the GKS has been well developed to solve for multi-componentflows [21–24, 30, 33, 50]. A second-order gas-kinetic scheme for multicomponent flow was2resented in [23, 50] based on the BGK equation for each component with its own equilibriumstate. Chen and Jiang proposed a non-oscillatory kinetic (NOK) scheme for the ideal gas[6] and stiffened gas [5]. Ni and Sun [30] proposed a γ -DGBGK scheme for compressiblemulticonponent flows simulation. Recently, an improved GKS for multicomponent flows [22]is proposed to increase the computational efficiency. In these papers mentioned above, mostof them are the second-order schemes at most. The works in [24, 30, 33] can achieve highorder, but only are applied to the ideal gas or stiffened gas.For the diffuse interface method, the numerical diffusion may lead to a very bad repre-sentation of the interfaces, especially when long time computations are needed. A way tocircumvent the numerical diffusion is to adopt a high order method, such as spectral volumemethod [24], weighted essentially non-oscillatory (WENO) method [12, 17, 31, 38], discon-tinuous Galerkin (DG) method which we are interested in. DG method has been appliedto solve a variety of different models [7, 14–16, 39]. There exist a few research works inthe 5-equation model. Saleem, Ali and Qamar [39] adopted the second-order Runge-KuttaDG (RKDG) method for solving the reduced 5-equation model [18]. In their work, the Lax-Friedrichs (LF) flux and the local LF flux were used to compute the numerical flux and aWENO limiter was utilized to eliminate oscillations at discontinuities. However, one canobserve that the velocity and pressure produce the oscillations at the interface from the re-sults of the interface only problem since the limiter was applied to the conservative variables.Gryngarten and Menon [15] applied the local DG method [51] to the 5-equation model [3]with Peng-Robinson EOS, where the non-conservative equations were rewritten into con-servative formula with source terms. A moment limiter [19] was applied to the conservedand primitive variables. The numerical flux used in their study is the HLLC approximateRiemann solver for the conservation of mass, momentum and energy. But the additionalequations must be solved which increases the computations.In this paper, a high order quasi-conservative DG method for compressible multi-componentflows with Mie-Gr¨uneisen EOS based on the 5-equation model [3] is developed. The method3an obtain the high order in smooth regions, keep oscillation-free at discontinuities, includingthe material interface, which is different from the work in [39], and guarantee the volumefraction in [0,1]. In addition, we do not need the extra equations to solve and reduce thecomputations compared to the work in [15]. Following the idea of the quasi-conservativemethod introduced by Abgrall [1], the quasi-conservative DG method with NOK flux hasthree steps. Firstly, we adopt DG method to discretize the conservative equations in space.In order to treat Mie-Gr¨uneisen EOS, the NOK flux [5, 24] is utilized to compute the nu-merical flux in our work instead of the traditional numerical flux. Secondly, according tothe discretizations of the conservative equations, the necessary condition that avoids theunphysical oscillation near the material interfaces is derived, which is also the discretizationmethod for the volume fraction equation in (2.1). At last, the new multi-resolution WENOlimiter [55] is employed to prevent the oscillations at discontinuities. In order to keep thepressure and velocity oscillation-free at the interfaces, we applied the limiter to the primitivevariables as in [15] and the maximum-principle-satisfying limiter developed by Zhang andShu [52, 53] is applied to ensure that the volume fraction does not go out of the range. Thus,a high order quasi-conservative discontinuous Galerkin method for multi-component flowswith Mie-Gr¨uneisen EOS using the NOK flux is developed.The organization of the paper is as follows. The governing equations and EOSs aredescribed in Section 2. In Section 3, the DG method for multi-component flows, identificationof troubled cells, and limiters are described in detail. One- and two-dimensional numericalexamples are presented to demonstrate the accuracy and the oscillation-free of the methodin Section 4. In Section 5, the conclusions are given. In one dimension, the 5-equation model for an immiscible two-material compressible flow[3] is considered, which is in the form of: 4 W t + F ( W ) x = 0 , ∂Y∂t + v ∂Y∂x = 0 , (2.1)where W = ( ρ Y , ρ Y , ρv, E ) T and F ( W ) = ( ρ Y v, ρ Y v, ρv + P, v ( E + P )) T ; ρ and ρ are the partial density of the fluids 1 and 2, respectively; P is the pressure, v is the velocityand E = ρe + ρv is the total energy with ρe being the internal energy; Y = Y is thevolume fraction of fluid 1, lies in the interval [0 , Y + Y = 1. The total density,momentum and energy of the mixture are defined as ρ = ρ Y + ρ Y , ρv = Y ρ v + Y ρ v , (2.2) E = Y ρ e + 12 ρ Y ( v ) + Y ρ e + 12 ρ Y ( v ) . (2.3)In order to close the equation (2.1), a mixture EOS is needed. In this work, each of thefluids is modeled by Mie-Gr¨uneisen EOS, i.e. P ( ρ, e ) = Γ( ρ )[ ρe − ρe ref ( ρ )] + P ref ( ρ ) , where Γ is the Gr¨uneisen coefficient, P ref ( ρ ) and e ref ( ρ ) are the reference pressure andinternal energy, respectively. This is a general EOS since it can produce the different typesof EOSs:(1) Ideal gas EOS Γ( ρ ) = γ − ,P ref ( ρ ) = 0 ,e ref ( ρ ) = 0; (2.4)(2) Stiffened gas EOS Γ( ρ ) = γ − ,P ref ( ρ ) = − γB,e ref ( ρ ) = 0; (2.5)(3) Jones-Wilkins-Lee EOS (JWL EOS) Γ( ρ ) = Γ ,P ref ( ρ ) = AR ρ exp( − R ρ ρ ) + BR ρ exp( − R ρ ρ ) − e ,e ref ( ρ ) = A exp( − R ρ ρ ) + B exp( − R ρ ρ ); (2.6)5able 2.1: Material-dependent quantities used in this paper. See [44] for other materialparameters. JWL EOS ρ ( kg/m ) A (GPa) B (GPa) R R Γ α Water 1004 1582 -4.67 8.94 1.45 1.17 0CC EOS ρ ( kg/m ) A (GPa) B (GPa) (cid:15) (cid:15) Γ α Copper 8900 145.67 147.75 2.99 1.99 2 0TNT 1840 12.87 13.42 4.1 3.1 0.93 0Shock EOS ρ ( kg/m ) c ( m/s ) s Γ α P e Molybdenum 9961 4770 1.43 2.56 1 0 0MORB 2660 2100 1.68 1.18 1 0 0where A , R , ρ , B , R and e are the material-dependent parameters.(4) Cochran-Chan EOS (CC EOS) Γ( ρ ) = Γ ,P ref ( ρ ) = − A ρ (1 − (cid:15) ) [( ρ ρ ) − (cid:15) −
1] + B ρ (1 − (cid:15) ) [( ρ ρ ) − (cid:15) − − e ,e ref ( ρ ) = A ( ρ ρ ) − (cid:15) − B ( ρ ρ ) − (cid:15) ; (2.7)where A , (cid:15) , ρ , B , (cid:15) and e are the material-dependent parameters.(5) Shock-Wave EOS (Shock EOS) Γ( ρ ) = Γ ( VV ) α , V = ρ , V = V ,P ref ( ρ ) = P + c ( V − V )[ V − s ( V − V )] ,e ref ( ρ ) = e + [ P ref ( ρ ) + P ]( V − V ); (2.8)where s, c , ρ , α, P and e are the material-dependent parameters.A wide range of real materials can be modeled by these EOSs. A typical set of numericalvalues for some sample materials is listed in Table 2.1. Mie-Gr¨uneisen EOS can be alsorewritten as P ( ρ, ρe ) = ( γ ( ρ ) − ρe − γ ( ρ ) π ( ρ ) , where γ ( ρ ) = Γ( ρ ) + 1 , π ( ρ ) = Γ( ρ ) ρe ref ( ρ ) − P ref ( ρ )Γ( ρ )+1 . To close system (2.1) for the mixing cells, the isobaric closure assumption [3] is sup-plemented here, which assumes that there is no pressure jump across a material interface,i.e. P = P = P. ρe = (cid:88) k Y k ρ k e k = (cid:88) k Y k P k + γ k ( ρ k ) π k ( ρ k ) γ k ( ρ k ) − . (2.9)Using P = P = P , we can obtain P + γπγ − ρe = (cid:88) k Y k P + γ k ( ρ k ) π k ( ρ k ) γ k ( ρ k ) − . Therefore, we have the two following equations:1 γ − (cid:88) k Y k γ k ( ρ k ) − , γπγ − (cid:88) k Y k γ k ( ρ k ) π k ( ρ k ) γ k ( ρ k ) − . (2.10)Finally, the mixing sound speed [3] can be written as follows: c = ( γ − (cid:88) k z k c k γ k − , where c k is the sound speed of the k th material and z k is mass fraction of fluid k defined as z k = ρ k Y k (cid:80) l ρ l Y l . In this section, we describe a quasi-conservative RKDG method for the numerical solutionof compressible multi-components in the form of (2.1) on a uniform mesh, which containsthree steps.Step 1. Discretize the conservative equations in (2.1) using RKDG method [8–11] withNOK flux which is suitable for Mie-Gr¨uneisen EOS.Step 2. Following the idea of Abgrall’s quasi-conservative method, define a numericalscheme for the equation of volume fraction which can prevent the oscillation of pressure andvelocity near the material interfaces.Step 3. Add the limiters, which include the limiters for oscillation-free near the discon-tinuities and the maximum-principle-satisfying limiters for the volume fraction.7 .1 DG method for the conservative equations
For simplicity, we take one dimension for example. For two dimensions, we can im-plement it similarly. Assume the domain Ω is divided into N nonoverlapping cells { I j =( x j − , x j + ) , j = 1 , · · · , N } , and ∆ x = x j + − x j − . The DG finite element space is definedas V kh = { p ( x, t ) : p | I j ∈ P k ( I j ) } , where P k ( I j ) is the space of polynomials of degree ≤ k defined on I j . Notice that P k ( I j )can be expressed as P k ( I j ) = span { ϕ ( x ) , · · · , ϕ L ( x ) } , where L = k + 1 for one dimensional case, and { ϕ ( x ) , · · · , ϕ L ( x ) } is a basis of P k ( I j ). Thefirst three basis functions in one dimension we employ on the cell I j are ϕ ( x ) = 1 , ϕ ( x ) = x − x j ∆ x , ϕ ( x ) = ( x − x j ∆ x ) − , ∀ x ∈ I j . (3.1)The semi-discrete DG approximation for the conservative equations in (2.1) is to find thenumerical solution u h ( · , t ) ∈ V kh , t ∈ (0 , T ] such that (cid:90) I j ∂W h ∂t ψdx + ( F ψ ) | I j − (cid:90) I j F ψ x dx = 0 , ∀ ψ ∈ P k ( I j ) , (3.2)Expressing W h as W h ( x, t ) = L (cid:88) l =1 W ( l ) j ( t ) ϕ l ( x ) , ∀ x ∈ I j , (3.3)applying the Gauss quadrature rule to the third terms in (3.2) and replacing the flux F bythe numerical flux ˆ F , we obtain (cid:82) I j ( L (cid:80) l =1 dW ( l ) j ( t ) dt ϕ l ( x )) ψdx + ˆ F j + ψ − j + − ˆ F j − ψ + j − − (cid:80) G F ( W h ( x G )) w G ψ x ( x G ) | I j | = 0 , ∀ ψ ∈ V kh , (cid:82) I j ( W h ( x, − W ( x, ψdx = 0 , ∀ ψ ∈ V kh , (3.4)where | I j | is the volume of the element I j , and x G and w G represent the Gaussian points andthe weights on I j , respectively. The summation (cid:80) G is taken over the Gauss points on I j . The8umerical flux has the form ˆ F j + = ˆ F ( u − j + , u + j + ), and u − j + = u ( x − j + ) and u + j + = u ( x + j + )are defined as the values from the left and right limit of x j + , respectively.In this work, the NOK flux [5, 24] which can deal with the general EOS, is employed. Itconsists of two parts ˆ F j + = ηF Kj + + (1 − η ) F Ej + , where η ∈ [0 , . The non-equilibrium part is F Kj + = F + j + + F − j + , where F ± j + = < u > j + ,L/R ρ Y ρ Y ρvE ∓ j + + P ∓ j + < u > j + ,L/R P ∓ j + < u > j + ,L/R + P ∓ j + v ∓ j + < u > j + ,L/R , (3.5) < u > j + ,L/R = 12 erfc( ∓ (cid:113) λ j + v ∓ j + ) ,< u > j + ,L/R = v ∓ j + < u > j + ,L/R ±
12 e − λ j + 12 ( v ∓ j + 12 ) (cid:113) πλ j + ,λ j + = min { c − j + ) , c + j + ) } ,c is the speed of sound.In order to avoid oscillations of the pressure and velocity near a contact discontinuity,the equilibrium part should satisfy the consistent condition. The primitive variables arecomputed by ¯ ρ ¯ ρ ¯ v ¯ P ¯ Y j + = ( ρ ) − j + < u > j + ,L +( ρ ) + j + < u > j + ,R ( ρ ) − j + < u > j + ,L +( ρ ) + j + < u > j + ,R < u > j + ,L + < u > j + ,R P − j + < u > j + ,L + P + j + < u > j + ,R ( Y ) − j + < u > j + ,L +( Y ) + j + < u > j + ,R . (3.6)9hen we take F Ej + = ¯ ρ ¯ Y ¯ v ¯ ρ ¯ Y ¯ v ¯ ρ ¯ v + ¯ P ¯ v ( ¯ E + ¯ P ) j + , (3.7)and ¯ E is determined by EOS. Following the procedure of Abgrall [1] to avoid the oscillations of the pressure and velocity,we consider the interface only problem, and assume the velocity v and the pressure P areconstants, i.e. v = v , P = P . Thus, < u > L + < u > R = 1 and < u > L + < u > R = v in the NOK flux. Firstly, we introduce˜ Z j + = Z − j + < u > j + ,L + Z + j + < u > j + ,R ,a l = (cid:90) I j ( ϕ l ( x )) dx, l = 1 , · · · , L for notation. Then from the current spatial discretization, the semi-discretized scheme ofthe continuity equation can be written in the form as following: d ( ρ Y ) ( l ) j dt = − a l [ η ( (cid:94) ( ρ Y ) j + ϕ l ( x − j + ) − (cid:94) ( ρ Y ) j − ϕ l ( x + j − )) + (1 − η ) v (( ρ Y ) j + ϕ l ( x − j + ) − ( ρ Y ) j − ϕ l ( x + j − )) − v (cid:90) I j ( ρ Y )( ϕ l ( x )) x dx ] , l = 1 , · · · , L. (3.8) d ( ρ Y ) ( l ) j dt = − a l [ η ( (cid:94) ( ρ Y ) j + ϕ l ( x − j + ) − (cid:94) ( ρ Y ) j − ϕ l ( x + j − )) + (1 − η ) v (( ρ Y ) j + ϕ l ( x − j + ) − ( ρ Y ) j − ϕ l ( x + j − )) − v (cid:90) I j ( ρ Y )( ϕ l ( x )) x dx ] , l = 1 , · · · , L. (3.9)Since ρ = ρ Y + ρ Y , we can get dρ ( l ) j dt = − a l [ η ( ˜ ρ j + ϕ l ( x − j + ) − ˜ ρ j − ϕ l ( x + j − )) + (1 − η ) v ( ¯ ρ j + ϕ l ( x − j + ) − ¯ ρ j − ϕ l ( x + j − )) − v (cid:90) I j ρ ( ϕ l ( x )) x dx ] , l = 1 , · · · , L. (3.10)10imilarly, the discretization for the momentum equation can be written as d ( ρv ) ( l ) j dt = − a l [ ηv ( ˜ ρ j + ϕ l ( x − j + ) − ˜ ρ j − ϕ l ( x + j − )) + (1 − η ) v ( ¯ ρ j + ϕ l ( x − j + ) − ¯ ρ j − ϕ l ( x + j − )) − v (cid:90) I j ρ ( ϕ l ( x )) x dx ] , l = 1 , · · · , L. (3.11)Based on the equation (3.10) and (3.11), we can derive dv ( l ) j dt = 0. The discretization forenergy is dE ( l ) j dt = − a l [ η ( ˜ E j + ϕ l ( x − j + ) − ˜ E j − ϕ l ( x + j − )) + (1 − η ) v ( ¯ E j + ϕ l ( x − j + ) − ¯ E j − ϕ l ( x + j − )) − v (cid:90) I j E ( ϕ l ( x )) x dx ] , l = 1 , · · · , L. (3.12)Inserting the equation of state to (3.12), we can observe that the pressure keeps constantat the next time if the following condition is satisfied, dκ ( l ) j dt = − a l [ η (˜ κ j + ϕ l ( x − j + ) − ˜ κ j − ϕ l ( x + j − )) + (1 − η ) v (¯ κ j + ϕ l ( x − j + ) − ¯ κ j − ϕ l ( x + j − )) − v (cid:90) I j κ ( ϕ l ( x )) x dx ] , l = 1 , · · · , L, (3.13) dχ ( l ) j dt = − a l [ η ( ˜ χ j + ϕ l ( x − j + ) − ˜ χ j − ϕ l ( x + j − )) + (1 − η ) v ( ¯ χ j + ϕ l ( x − j + ) − ¯ χ j − ϕ l ( x + j − )) − v (cid:90) I j χ ( ϕ l ( x )) x dx ] , l = 1 , · · · , L, (3.14)where κ = γ − and χ = γπγ − . Following the work in [5, 24, 44], if the condition is replacedby dY ( l ) j dt = − a l [ η ( ˜ Y j + ϕ l ( x − j + ) − ˜ Y j − ϕ l ( x + j − )) + (1 − η ) v ( ¯ Y j + ϕ l ( x − j + ) − ¯ Y j − ϕ l ( x + j − )) − v (cid:90) I j Y ( ϕ l ( x )) x dx ] , l = 1 , · · · , L, (3.15)the conditions (3.13)-(3.14) will be satisfied. It is easy to observe that the conditions (3.15)can be viewed as the discretization of the equation Y t + ( vY ) x = 0 . Y t + vY x = Y t + ( vY ) x − Y v x = 0 , and discretized as follows, dY ( l ) j dt = − a l [ η ( ˜ Y j + ϕ l ( x − j + ) − ˜ Y j − ϕ l ( x + j − )) + (1 − η ) v ( ¯ Y j + ϕ l ( x − j + ) − ¯ Y j − ϕ l ( x + j − )) − v (cid:90) I j Y ( ϕ l ( x )) x dx − Y ( x j )(ˆ v j + ϕ l ( x − j + ) − ˆ v j − ϕ l ( x + j − ) − (cid:90) I j v ( ϕ l ( x )) x dx )] ,l = 1 , · · · , L, (3.16)where ˆ v j + = < u > j + ,L + < u > j + ,R . Then (3.16) is degenerated to (3.15) near thecontact discontinuities. Thus, the velocity and the pressure are oscillation-free.Finally the semi-discrete schemes (3.4) and (3.16) are discretized in time. Here, we usean explicit, the third order TVD Runge-Kutta scheme [42]. Casting (3.4) and (3.16) in theform ∂u h ∂t = L h ( u h , t ) , the scheme reads as u ∗ h = u nh + ∆ t n L h ( u nh , t n ) ,u ∗∗ h = 34 u nh + 14 ( u ∗ h + ∆ t n L h ( u ∗ h , t n + ∆ t n )) , (3.17) u n +1 h = 13 u nh + 23 ( u ∗∗ h + ∆ t n L h ( u ∗∗ h , t n + 12 ∆ t n )) . In this subsection, two kinds of limiters are described briefly. One is the limiter tokeep oscillation-free at discontinuities and the other one is the maximum-pricinple-satisfyinglimiter for the volume fraction, since it should satisfy Y ∈ [0 , . .3.1 The limiter to control oscillations As is well known, nonlinear limiters must be applied to control the spurious oscillationsin the numerical solution for strong shocks, which has two steps following the work of Qiuand Shu [37].Step 1. We identify the “troubled cells” using the minmod-type TVB limiter as in[28, 36, 37]. All the primitive variables are taken as the indicator variables.Step 2. We add the nonlinear limiters in the troubled cells. In this paper, the new typeof multi-resolution WENO limiters developed in [55, 56] is adopted. In order to keep thepressure non-oscillatory, we limit the primitive variables component-wisely here.Next, we describe the method to detect “troubled cells” and the new multi-resolutionWENO limiter briefly.For simplicity, we assume u ( x ) ∈ V kh is the primitive variable such as ρ , ρ , v, P, Y on I j .Denote u − j + = u (1) j + ˜ u j , u + j − = u (1) j − ˜˜ u j , ¯ u j = 1∆ x (cid:90) I j udx. These values are modified by the standard minmod limiter˜ u modj = ˜ m (˜ u j , ∆ + ¯ u j , ∆ − ¯ u j ) , ˜˜ u modj = ˜ m (˜˜ u j , ∆ + ¯ u j , ∆ − ¯ u j ) , where ∆ + ¯ u j = ¯ u j +1 − ¯ u j , ∆ − ¯ u j = ¯ u j − ¯ u j − , and ˜ m is defined as (cid:101) m ( a , a , a ) = (cid:26) a , if | a | ≤ M ∆ x ,m ( a , a , a ) , otherwise , (3.18) m ( a , a , a ) = (cid:26) sign( a ) min( | a | , | a | , | a | ) , if sign( a ) = sign( a ) = sign( a ) , , otherwise , (3.19)and M > M depends on the solution of the problem; see, e.g.,[10] for detailed discussion. We use M = 1 in our computation. Finally, I j is marked as atroubled cell for further reconstructions if one of the minmod functions does not return thefirst argument.In order to keep the velocity and pressure non-oscillatory, we limit the primitive vari-ables component-wisely here. Then the limited primitive variables are used to compute the13umerical fluxes and evolve the equation to obtain the new solutions at the next time in theconservative form. Assume I j is a troubled cell. The procedure of the limiting for the scalarcase [55] is given in the following. Step 1.
Define a series of polynomials of different degrees on the troubled cell I j . Step 1.1.
For a second-order spatial approximation, a zeroth degree polynomial q ( x )and a linear polynomial q ( x ) are constructed, which satisfy (cid:90) I j q ( x ) ϕ ( x ) dx = (cid:90) I j u ( x ) ϕ ( x ) dx, and (cid:90) I j q ( x ) ϕ l ( x ) dx = (cid:90) I j u ( x ) ϕ l ( x ) dx, l = 1 , . Step 1.2.
For a third-order spatial approximation, a quadratic polynomial q ( x ) isconstructed which satisfies (cid:90) I j q ( x ) ϕ l ( x ) dx = (cid:90) I j u ( x ) ϕ l ( x ) dx, l = 1 , , . Step 2.
Get equivalent expressions for these constructed polynomials of different degrees.
Step 2.1.
For the second-order approximation, we obtain a polynomial p , ( x ) by p , ( x ) = 1 γ , q ( x ) − γ , γ , p , ( x )with γ , + γ , = 1 and γ , (cid:54) = 0 , where p , ( x ) = q ( x ) . Step 2.2.
For the third-order approximation, we define p , ( x ) = ω , p , ( x )+ ω , p , ( x ) , and obtain a polynomial p , ( x ) through p , ( x ) = 1 γ , q ( x ) − γ , γ , p , ( x )with γ , + γ , = 1 and γ , (cid:54) = 0 . In these expressions, γ l,l and ω l,l for l = l − , l ; l = 1 , · · · , k are the linear weightsand nonlinear weights, respectively. Step 3.
Compute the smoothness indicators β l,l by β l,l = l (cid:88) s =1 (cid:90) I j ∆ x s − ( d s dx s p l,l ( x )) dx, l = l − , l ; l = 1 , . (3.20)14owever, β , can not be computed by (3.20), which is defined below. We first define thelinear polynomial q j − ( x ) on I j − by (cid:90) I j − q j − ( x ) ϕ l ( x ) dx = (cid:90) I j − u ( x ) ϕ l ( x ) dx, l = 1 , . and similarly, the linear polynomial q j +1 ( x ) on I j +1 by (cid:90) I j +1 q j +1 ( x ) ϕ l ( x ) dx = (cid:90) I j +1 u ( x ) ϕ l ( x ) dx, l = 1 , . Then, the smoothness indicators are computed by ζ j − = (cid:90) I j ∆ x ( ddx q j − ( x )) dx, ζ j +1 = (cid:90) I j ∆ x ( ddx q j +1 ( x )) dx. Thus, β , is defined as β , = min( ζ j − , ζ j +1 ) . Step 4.
Compute the nonlinear weights ω l ,l = ¯ ω l ,l l (cid:80) s =1 ¯ ω s,l , ¯ ω l ,l = γ l ,l (1 + τ l υ + β l ,l ) , l = l − , l ; l = 1 , , where τ l = ( β l ,l − β l − ,l ) , l = 1 , , and υ is set to be 10 − in all the computations. Step 5.
Finally the new constructed polynomial u new ( x ) on the cell I j is given by u new ( x ) = l (cid:88) l = l − ω l,l p l,l ( x ) , l = 1 , , for the second-order, third-order, respectively. After limiting described above, the volume fraction Y may still have a non-valid value,such as Y <
Y > Y ( x ) is the polynomial defined on I j and ¯ Y is the cell averageon I j . Then we modify Y ( x ) such that Y ( x ) ∈ [ (cid:15), − (cid:15) ] for all x ∈ S where S is the set of15he Legendre Gauss-Lobatto quadrature points for I j . For all j , assume ¯ Y ∈ [ (cid:15), − (cid:15) ], weuse the modified polynomial ˜ Y ( x ) instead of Y ( x ), i.e.,˜ Y ( x ) = θ ( Y ( x ) − ¯ Y ) + ¯ Y , θ = min (cid:110)(cid:12)(cid:12)(cid:12) − (cid:15) − ¯ YY max − ¯ Y (cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12) (cid:15) − ¯ YY min − ¯ Y (cid:12)(cid:12)(cid:12) , (cid:111) , (3.21)where Y max = max x ∈ S Y ( x ) , Y min = min x ∈ S Y ( x ) . It is clear that the volume fraction ˜ Y ( x ) shouldbe in [ (cid:15), − (cid:15) ] after this limiting. The parameter (cid:15) is set to be 10 − in this work.For two-component flows, it is easy to see that the volume fraction of the fluid 2 alsostays in [ (cid:15), − (cid:15) ] due to Y = 1 − Y if the volume fraction Y = Y ( x ) of the fluid 1 is limitingusing the method described above. However, for more than two fluids, it is a little different.We take three-component flows for example. Assume Y ( x ), Y ( x ) and Y ( x ) are the volumefraction of the fluid 1, 2 and 3, respectively. The limiting procedure is given as follows. Step 1.
Let Y ( x ) = Y ( x ) + Y ( x ) and use the new volume fraction Y ( x ) to define theparameter θ θ = min (cid:110)(cid:12)(cid:12)(cid:12) − (cid:15) − ¯ Y Y ,max − ¯ Y (cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12) (cid:15) − ¯ Y Y ,min − ¯ Y (cid:12)(cid:12)(cid:12) , (cid:111) . Step 2.
Similarly, define the parameters θ i θ i = min (cid:110)(cid:12)(cid:12)(cid:12) − (cid:15) − ¯ Y i − Y i − ,max − ¯ Y i − (cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12) (cid:15) − ¯ Y i − Y i − ,min − ¯ Y i − (cid:12)(cid:12)(cid:12) , (cid:111) , i = 2 , . Step 3.
Finally, use the modified polynomials ˜ Y ( x ) and ˜ Y ( x ) instead of Y ( x ) and Y ( x ),i.e. ˜ Y i ( x ) = θ ( Y i ( x ) − ¯ Y i ) + ¯ Y i , θ = min (cid:110) θ , θ , θ , (cid:111) , i = 1 , . (3.22)From the procedure described above, we need to compute a common coefficient θ , then θ is applied to modify the volume fraction polynomials. Moreover, the procedure can beextended to more than three fluids easily. In this section we present numerical results obtained with the quasi-conservative DGmethod described in the previous sections for a selection of one- and two-dimensional exam-ples. Recall that the method has been described in one dimension. Its implementation in two16able 4.1: Example 4.1: Solution error with periodic boundary conditions and T = 1. k N
10 20 40 80 160 3201 L L L ∞ L L L ∞ P elements,0.15 for P elements. For the examples with Mie-Gr¨uneisen EOS, the material-dependentparameters are given in Table 2.1. Example 4.1
To assess the accuracy of the new method, we first consider a one-dimensionalconvection of change in volume fraction with the equation of state (2.5), which is also studiedin [15]. And the parameters are set to be γ = 1 . , γ = 1 . , B = 1 , B = 0 . The initialcondition is given by ρ ( x,
0) = 1 , v ( x,
0) = 1 , P ( x,
0) = 1 , Y ( x,
0) = 0 . .
499 sin( πx )with a periodic boundary condition. The computational domain is on (0,2). We computethe solution up to T = 1. The error of the volume fraction is listed in Table 4.1, which showsthe convergence of the second order for P elements, the third order for P elements for thequasi-conservative DG method. Example 4.2
To verify the non-oscillation property for the pressure and velocity fields,17n this example we consider the interface only problem with the initial condition given by( ρ, v, P, γ, B ) = (cid:40) (1 , , , . , , x (cid:54) , (0 . , , , . , , x > . The stiffened gas equation of state (2.5) is used and the numerical results with 100 pointsat T = 2 are plotted in Fig 4.1.From the figure, one can observe that the new DG method can preserve the oscillation-free property of the pressure and velocity at the material interface. Moreover, from theclose-up of the density at the interface in Fig 4.1, it is clear that the high order method havebetter resolution than lower order method. Example 4.3
In this example the gas-liquid shock tube test with a strong shock waveis considered and the stiffened gas EOS (2.5) is employed. This is a very challenging testcase with a strong shock wave since the shock and the material interface are close and thepressure ratio is excessively high. The initial condition is( ρ, v, P, γ, B ) = (cid:40) (10 , , , . , × ) , x (cid:54) . , (50 , , , . , , x > . , and the computational domain is (-0.2,1). The final time is T = 0 . P elements have betterresolution than ones with P elements. Example 4.4
In order to show that the quasi-conservative DG method works with Mie-Gr¨uneisen EOS, we test a two-component impact problem, which is also studied in [40, 44].In this problem, to model the material properties of the copper and solid explosive, the sameCC EOSs (2.7) are used, but with a different set of material-dependent quantities for each ofthem. At the beginning, the copper has an initial velocity of 1500 m/s, while the explosiveis at rest. The computational domain is (0 ,
1) and the initial condition is given by( ρ, v, P ) = (cid:40) (8900 , , ) , x (cid:54) . , (1840 , , ) , x > . . r ho
4 2 0 2 40.20.40.60.81 exactP1,N=100P2,N=100 x r ho exactP1,N=100P2,N=100 x u
4 2 0 2 40.60.811.21.4 exactP1,N=100P2,N=100 x P
4 2 0 2 40.60.811.21.4 exactP1,N=100P2,N=100
Figure 4.1: Example 4.2 interface only problem. N = 10019 r ho exactP1,N=5000P2,N=5000 x r ho exactP1,N=5000P2,N=5000 x u exactP1,N=5000P2,N=5000 x P exactP1,N=5000P2,N=5000 Figure 4.2: Example 4.3 N = 500020he boundary conditions are constant states on both the left and right sides of the domain.The integration is stopped at T = 85 µs .The exact solution for this problem consists of a rightward-moving shock, a leftward-moving shock and a material interface in between. The results with 200 uniform points aredemonstrated in Fig 4.3, where the solid line is the fine grid solution computed by ∆ x = with P elements. From that one can observe these nonlinear structures are all resolved well. Example 4.5
In order to test the accuracy in two-dimensional case, similar to Example4.1, we first consider a two-dimensional convection of change in volume fraction with theequation of state (2.5) and the parameters are γ = 1 . , γ = 1 . , B = 1 , B = 0 . The initialcondition is given by ρ ( x,
0) = 1 , µ ( x,
0) = 1 , ν ( x,
0) = 1 , P ( x,
0) = 1 , Y ( x,
0) = 0 . .
499 sin( π ( x + y ))with a periodic boundary condition, where µ and ν are the velocities in the x -direction and y -direction, respectively. The computational domain is taken as ( x, y ) ∈ (0 , × (0 , T = 1. The error of the volume fraction is listed in Table 4.2,which shows the convergence of the second order for P elements, the third order for P elements for the quasi-conservative DG method in two dimensions. Example 4.6
To show the performance of our method with high pressure ratio in twodimensions, we consider the simulation of a model underwater explosion problem [20, 46].In this test, the computation domain is taken as ( x, y ) ∈ ( − , × ( − . , y = 0 and the center of a circular gas bubblewith the radius 0.12 in the water is located at (0 , − . r ho ’exact’P1,N=200P2,N=200 x r ho ’exact’P1,N=200P2,N=200 x u ’exact’P1,N=200P2,N=200 x u ’exact’P1,N=200P2,N=200 x P ’exact’P1,N=200P2,N=200 x P ’exact’P1,N=200P2,N=200 Figure 4.3: Example 4.4 N = 200, Left: P elements; Right: P elements22able 4.2: Example 4.5: Solution error with periodic boundary conditions and T = 1. k N × M ×
10 20 ×
20 40 ×
40 80 ×
80 160 ×
160 320 × L L L ∞ L L L ∞ ρ, µ, ν, P, γ, B ) = (1 . , , , , . , , y > , (1250 , , , , . , , x + ( y + 0 . (cid:54) . , (1000 , , , , . , × ) , else. And the reflecting boundary conditions are employed on the bottom of the domain, whilenon-reflecting boundary conditions are used on the remaining sides [20].From the initial condition, it is obvious that both the gas and water are in a stationaryposition at the beginning, but due to the pressure difference between the fluids, breakingof the bubble results in a circularly outward-going shock wave in water, an inward-goingrarefaction wave in gas, and an interface lying in between that separates the gas and thewater. Soon after, this shock wave is diffracted through the nearby air-water surface, causingthe subsequent deform of the interface topology from a circle to oval-like shape.The contours of the density and pressure are plotted in Figs. 4.4 and 4.6 at four differenttimes T = 0 . , . , . . ×
400 mesh.From the density and pressure plots, one can clearly see that the improvement on the use ofthe high order method to the sharpness near the interfaces. The cross-sections of the densityand pressure for the same run along line x = 0 are shown in Figs. 4.5 and 4.7, which give23ome information about the differences between P and P elements at the selected times. Example 4.7
To show our method works with complex Mie-Gr¨uneisen EOS in twodimensions, we are concerned with interaction of a shock in molybdenum with a block ofencapsulated mid-ocean ridge basalt (MORB) liquid [20, 29, 44]. The computational domainis set to be an unit square. A Mach 1.163 rightward-moving shock is located at x = 0 . . , . × [0 , . ρ, µ, ν, P ) = (2260 , , , , and outside the MORB, the state variables in the preshock region are given by( ρ, µ, ν, P ) = (9961 , , , , and the state variables in the postshock region are( ρ, µ, ν, P ) = (11042 , , , × ) . The contours of the density and pressure at two different selected times T = 50 and T = 100 µs with a uniform 400 ×
400 mesh are illustrated in Figs 4.8 and 4.9. In the densityplot, we can see the incident shock in molybdenum and transmitted shock in MORB with theformer moving faster than the latter at T = 50 µs . And the transmitted shock has not passedthe MORB block completely at T = 100 µs . In addition, the structure of diffraction of theshock by MORB is well captured in the pressure graphs. From the displayed figures, oneis easy to observe that the improved resolution of the numerical solution near the interfacewhen P elements is adopted in the test. Example 4.8
Finally, we consider the three-component impact problem in two dimen-sions [44]. The computation domain is taken as ( x, y ) ∈ (0 , × (0 , T = 0.2, 0.4, 0.8, 1.2ms. Left: P elements; Right: P elements 25 r ho Y r ho Y r ho Y r ho Y r ho Y r ho Y r ho Y r ho Figure 4.5: Cross-sectional plots of the results in Fig. 4.4 along x = 0 . From top to bottom: T = 0.2, 0.4, 0.8, 1.2 ms. Left: P elements; Right: P elements26igure 4.6: Example 4.6 The pressure contours. From top to bottom: T = 0.2, 0.4, 0.8, 1.2ms. Left: P elements; Right: P elements 27 P Y P Y P Y P Y P Y P Y P Y P Figure 4.7: Cross-sectional plots of the results in Fig. 4.6 along x = 0 . From top to bottom: T = 0.2, 0.4, 0.8, 1.2 ms. Left: P elements; Right: P elements28 Y X Y X Y X Y Figure 4.8: Example 4.7 N = 400 × T = 50 µ s; Bottom: T =100 µ s. Left: P elements; Right: P elements29 Y X Y X Y X Y Figure 4.9: Example 4.7 N = 400 × T = 50 µ s; Bottom: T =100 µ s. Left: P elements; Right: P elements30oing copper plate traveling vertically in a shock tube with speed 1500 m/s from right toleft in region x (cid:62) .
6, while in region x < .
6, we have a solid inert explosive on the top anda liquid water on the bottom separated by the interface at y = 0 .
5. The solid inert explosiveand liquid water are at rest and all three fluid components are in the usual atmosphericcondition initially throughout the domain. The copper and explosive are modeled by theCC EOS (2.7) while the water is modeled by JWL EOS (2.6).The numerical results are shown in Figs. 4.10-4.13. Clearly, we observe that the shockspeed in explosive is larger than the one in water from the figures since the acoustic impedanceof explosive is greater than the one for the water.
We have presented a high order quasi-conservative DG method for compressible multi-component flows with Mie-Gr¨uneisen equation of state based on the 5-equation model in theprevious sections. In this paper the NOK flux is used to compute the numerical flux, whichis free from constructing Riemann solver. Then, a DG scheme is defined for the volumefraction equations according to the procedure of the quasi-conservative method, which cankeep the velocity and pressure oscillation-free at the interface. In addition, a maximum-pricinple-satisfying limiter is employed to ensure that the volume fraction does not go out ofthe range. Numerical results in one and two dimensions shown in the paper demonstrate theability of the method to capture shocks and material interfaces and be high order in smoothregions. In the future, we plan to further extend the method to the unstructured mesh. Inorder to reduce the numerical diffusion further, we will extend the quasi-Lagrangian movingDG method [28] to the 5-equation model of multi-component flows.
Acknowledgements
This work is supported by National Natural Science Foundation of China (Grant Nos.11671050), Science Challenge Project (Grant No. TZ2016002), and National Natural Science31igure 4.10: Example 4.8 Schlieren-type images for the density using N = 200 × T = 50 µ s; Bottom: T = 100 µ s. Left: P elements; Right: P elements32igure 4.11: Example 4.8 Schlieren-type images for the pressure using N = 200 × T = 50 µ s; Bottom: T = 100 µ s. Left: P elements; Right: P elements33 r ho X r ho X r ho X r ho Figure 4.12: Example 4.8 The cross-sectional plots of the results shown in Fig. 4.10 along y = 0 .
4. Top: T = 50 µ s; Bottom: T = 100 µ s. Left: P elements; Right: P elements34 P X P X P X P Figure 4.13: Example 4.8 The cross-sectional plots of the results shown in Fig. 4.11 along y = 0 .
4. Top: T = 50 µ s; Bottom: T = 100 µ s. Left: P elements; Right: P elements35oundation of China (Grant Nos. U1630247). References [1] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations:A quasi conservative approach,
J. Comput. Phys.
125 (1996), 150-160.[2] R. Abgrall, S. Karni, Computations of compressible multifluids,
J. Comput. Phys.
J. Comput. Phys.
181 (2002), 577-616.[4] I. L. Chen, J. Glimm, Front tracking for gas dynamics,
J. Comput. Phys.
62 (1986),83-110.[5] Y. Chen, S. Jiang, A non-oscillatory kinetic scheme for multi-component flows with theequation of state for a stiffened gas,
J. Comput. Math.
29 (2011), 661-683.[6] Y. Chen, S. Jiang, Modified kinetic flux vector splitting schemes for compressible flows,
J. Comput. Phys.
228 (2009), 3582-3604.[7] J. Cheng, F. Zhang, T. G. Liu, A discontinuous Galerkin method for the simulation ofcompressible gas-gas and gas-water two-medium flows,
J. Comput. Phys.
403 (2020),109059.[8] B. Cockburn, S. Hou, C.-W. Shu, The Runge-Kutta local projection discontinuousGalerkin finite element method for conservation laws IV: The multidimensional case,
Math. Comp.
54 (1990), 545-581.[9] B. Cockburn, S.-Y. Lin, C.-W. Shu, TVB Runge-Kutta local projection discontinuousGalerkin finite element method for conservation laws III: One dimensional systems,
J.Comput. Phys.
84 (1989), 90-113. 3610] B. Cockburn, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkinfinite element method for conservation laws II: General framework,
Math. Comp.
J. Comput. Phys.
141 (1998), 199-224.[12] V. Coralic, T. Colonius, Finite-volume WENO scheme for viscous compressible multi-component flows,
J. Comput. Phys.
274 (2014), 95-121.[13] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory Eulerian approach tointerfaces in multimaterial flows (the ghost fluid method),
J. Comput. Phys.
152 (1999),457-492.[14] E. Franquet, V. Perrier, Runge-Kutta discontinuous Galerkin method for interface flowswith a maximum preserving limiter,
Computers & Fluids
65 (2012), 2-7.[15] L. Gryngarten, S. Menon, A generalized approach for sub- and super-critical flowsusing the local discontinuous Galerkin method,
Comput. Methods Appl. Mech. Engrg.
253 (2013), 169-185.[16] M. T. Henry de Frahan, S. Varadan, E. Johnsen, A new limiting procedure for discon-tinuous Galerkin methods applied to compressible multiphase flows with shocks andinterfaces,
J. Comput. Phys.
280 (2015), 489-509.[17] E. Johnsen, T. Colonius, Implementation of WENO schemes in compressible multicom-ponent flow problems,
J. Comput. Phys.
219 (2006), 715-732.[18] J. J. Kreeft, B. Koren, A new formulation of Kapila’s five-equation model for com-pressible two-fluid flow, and its numerical treatment,
J. Comput. Phys.
229 (2010),6220-6242. 3719] L. Krivodonova, Limiters for high-order discontinuous Galerkin methods,
J. Comput.Phys.
226 (2007), 879-896.[20] T. S. Lee, J. G. Zheng, S. H. Winoto, An interface-capturing method for resolvingcompressible two-fluid flows with general equation of state,
Comm. Comput. Phys.
Comm. Comput. Phys.
25 (2019), 416-447.[22] Q. Li, An improved gas-kinetic scheme for multimaterial flows,
Comm. Comput. Phys.
27 (2020), 145-166.[23] Y. S. Lian, K. Xu, A gas-kinetic scheme for multimaterial flows and its application inchemical reactions,
J. Comput. Phys.
163 (2000), 349-375.[24] N. Liu, X. Xu, Y. Chen, High-order spectral volume scheme for multi-component flowsusing non-oscillatory kinetic flux,
Computers & Fluids
152 (2017), 120-133.[25] T.G. Liu, B.C. Khoo, C. W. Wang, The ghost fluid method for compressible gas-watersimulation,
J. Comput. Phys.
204 (2005), 193-221.[26] T.G. Liu, B.C. Khoo, K.S. Yeo, Ghost fluid method for strong shock impacting onmaterial interface,
J. Comput. Phys.
190 (2003), 651-681.[27] H. Lu, J. Zhu, D.H. Wang, N. Zhao, Runge-Kutta discontinuous Galerkin method withfront tracking method for solving the compressible two-medium flow,
Comput. Fluids
126 (2016), 1-11.[28] D. Luo, W. Huang, J. X. Qiu, A quasi-Lagrangian moving mesh discontinuous Galerkinmethod for hyperbolic conservation laws,
J. Comput. Phys.
396 (2019), 544-578.[29] G. Miller, E. G. Puckett, A high-order Godunov method for multiple condensed phases,
J. Comput. Phys.
128 (1996), 134-164. 3830] G. Ni, W. Sun, A γ -DGBGK scheme for compressible multi-fluids, Int. J. Numer. Meth.Fluids
66 (2011), 760-777.[31] T. Nonomura, K. Fujii, Characteristic finite-difference WENO scheme for multicompo-nent compressible fluid analysis: Overestimated quasi-conservative formulation main-taining equilibriums of velocity, pressure, and temperature,
J. Comput. Phys.
J. Comput. Phys.
79 (1988), 12-49.[33] L. Pan, J. Cheng, S. Wang, K. Xu, A two-stage fourth-order gas-kinetic scheme forcompressible multicomponent flows,
Comm. Comput. Phys.
22 (2017), 1123-1149.[34] J. X. Qiu, T. G. Liu, B. C.Khoo, Runge-Kutta discontinuous Galerkin methods forcompressible two-medium flow simulations: one-dimensional case,
J. Comput. Phys.
222 (2007), 353-373.[35] J. X. Qiu, T. G. Liu, B. C.Khoo, Simulations of compressible two-medium flow byRunge-Kutta discontinuous Galerkin methods with the ghost fluid method,
Comm.Comput. Phys.
J. Comput. Phys.
SIAM J. Sci. Comput.
26 (2005), 907-929.[38] A. Rehman, S. Qamar, High order finite-volume WENO scheme for five-equation modelof compressible two-fluid flow,
Comput. Math. with Appl.
76 (2018), 2648-2664.3939] M. R. Saleem, I. Ali, S. Qamar, Application of discontinuous Galerkin method for solvinga compressible five-equation two-phase flow model,
Results Phys.
J. Comput. Phys.
150 (1999), 425-467.[41] R. Saurel, F. Petitpas, R. A. Berry, Simple and efficient relaxation methods for interfacesseparating compressible fluids, cavitating flows and shocks in multiphase mixtures,
J.Comput. Phys.
228 (2009), 1678-1712.[42] C.-W. Shu, Total-variation-diminishing time discretizations,
SIAM J. Sci. Stat. Comput.
J. Comput. Phys.
142 (1998), 208-242.[44] K. M. Shyue, A fluid-mixture type algorithm for compressible multicomponent flowwith Mie-Gr¨uneisen equation of state,
J. Comput. Phys.
171 (2001), 678-707.[45] K. M. Shyue, A fluid-mixture type algorithm for compressible multicomponent flowwith Van der Waals equation of state,
J. Comput. Phys.
156 (1999), 43-88.[46] K. M. Shyue, A wave-propagation based volume tracking method for compressible mul-ticomponent flow in two space dimensions,
J. Comput. Phys.
215 (2006), 219-244.[47] C. W. Wang, T. G. Liu, C.-W. Shu, A real ghost fluid method for the simulation ofmultimedium compressible flow,
SIAM J. Sci. Comput.
28 (2006), 278-302.[48] C.-W. Wang, C.-W. Shu, An interface treating technique for compressible multi-mediumflow with Runge-Kutta discontinuous Galerkin method,
J. Comput. Phys.
229 (2010),8823-8843.[49] K. Xu, A kinetic method for hyperbolic-elliptic equations and its application in two-phase flow,
J. Comput. Phys.
166 (2001), 383-399.4050] K. Xu, BGK-based scheme for multicomponent flow calculations,
J. Comput. Phys.
J. Comput. Phys.
230 (2011), 232-244.[52] X. Zhang, C.-W. Shu, Maximum-principle-satisfying and positivity-preserving high or-der schemes for scalar conservation laws: survey and new developments,
Proc. R. Soc.A
467 (2011), 2752-2776.[53] X. Zhang, C.-W. Shu, On maximum-principle-satisfying high order schemes for scalarconservation laws,
J. Comput. Phys.
229 (2010), 3091-3120.[54] J. Zhu, J. Qiu, T. G. Liu, B. C. Khoo, High-order RKDG methods with WENO typelimiters and conservative interfacial procedure for one-dimensional compressible multi-medium flow simulations,
Appl. Numer. Math.
61 (2011), 554-580.[55] J. Zhu, J. Qiu, C.-W. Shu, High-order Runge-Kutta discontinuous Galerkin methodswith a new type of multi-resolution WENO limiters,
J. Comput. Phys.
404 (2020),109105.[56] J. Zhu, C.-W. Shu, J. Qiu, High-order Runge-Kutta discontinuous Galerkin methodswith a new type of multi-resolution WENO limiters on triangular meshes,