# Numerical Approximation of Compressible Two-Phase Six-Equation Model UsingCE/SE and RKDG Schemes

Authors: Omar Rabbani; Saqib Zia; Munshoor Ahmed; Asad Rehman; Ilyas Khan; Mulugeta Andualem
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

In this article, two-phase compressible six equation flow model is numerically investigated. The six-equation model consists of velocity, pressure equations, and also relaxation terms. An extra seventh equation is included describing the total energy of the mixture to ensure the correct treatment of the sharp discontinuities. The model is hyperbolic and poses numerous difficulties for numerical schemes. An efficient and well-balanced scheme can handle the numerical difficulties related to this model. The second order space-time CE/SE scheme is extended to solve the model. This scheme offers an effective numerical method for several continuum mechanics problems. The suggested scheme suppresses the numerical oscillations and dissipation effects. Several numerical test cases have been carried out to reveal the efficiency and performance of the proposed approach. The results are compared with the exact solution and also with Runge-Kutta Discontinuous Galerkin (RKDG) and central (NT) schemes.


---

## Body

## 1. Introduction

A lot of research work has been drafted on multiphase flows caused by their various applications in divergent scientific and engineering fields [1].In many cases, multiphase flow problems arise in technical and advanced heat transfer systems. They have industrial implementation in energy conversion sector, paper creation, food making, and medical applications. Furthermore, an effective and secure evaluation of nuclear reactors accident situation needs an accurate forecasting of two-phase flow structure. Moreover, two-fluid flow has some daily operations in geysers and boilers. Ishii classified two-phase flow into four different groups on the grounds of phase mixture in the flow [2]. The groups contain gas-liquid flow (i.e., gas droplet flows, bubble flows, and separated flows), gas-solid flow (i.e., fluidized flows and flow of gas particle), liquid-solid flow (i.e., flow of slurry and transport of sediment), and immiscible liquid-liquid flow. As two-phase flows have substantial implementations, appropriate mathematical templates are mandatory to be constructed and associated with hypothetical outcomes. The averaging methodology is usually the foundation of multiphase flows as there is intricate connection between the fluids in a multiphase flow approach [3]. Researchers and engineers have proposed a number of multiphase mathematical flow models in the past to study the physical behavior of such flows. One such type of those models is seven-equation models; they use different velocity, density, and pressure for two phases. The seventh equation is obtained by coupling of conservation laws and the convection equation for the interface motion. Baer and Nunziato [4] were the pioneers who formulated a solid gas two-phase seven-equation model; and later, this model was modified and further studied by Abgrall and Saurel [5, 6] and other investigators. Although, such models are treated as the best two-phase flow models. However, there exist a number of physical and numerical difficulties in such models. To overcome those problems present in such flow models, the investigators have suggested models with a smaller number of equations ranging from three to six [7, 8].In this work, we have primarily considered the six-equation two-phase flow model studied by [9]. As mentioned earlier that because of complex nature of the seven-equation model, Kapila et al. [7] extracted a reduced five-equation model from seven-equation model [4]. The extracted model [7] has been used successfully to study the two-phase flows. Despite the fact that this model is simple but the presence of nonconservative term and shocks make it difficult to obtain the convergent physical solution.To resolve difficulties, present in the Kapila’s model, Saurel and coauthors [9] proposed a new model by incorporating the nonequilibrium pressure terms present in the seven-equation model. This new model is the six-equation model with same velocity, separate pressure, and relaxation terms for each phase. To treat the shock correctly in the single phase limit, an additional equation of mixture energy has been incorporated in the model. Importantly, this new model is hyperbolic in nature with three wave propagation speeds. Moreover, the positivity of the volume fraction is an important feature present in the model. Saurel and coauthors [9] applied Godunov-type method and HLLC-type solver to solve the six-equation model. Zia et al. [10] applied KFVS and central-upwind schemes to study the model under discussion.In this article, a space-time CE/SE scheme [11] is extended to solve the reduced six-equation model equations. The space-time CE/SE scheme is different from the existing methods in literature. In CE/SE scheme, the space and time variables are discretized simultaneously, whereas in other methods, the space variable is discretized first, and the resulting system of ODE’s is solved by any other method. Because of the unique idea of discretization, the CE/SE scheme effectively handles the numerical complexities of the two-layer equations. This scheme has been successfully implemented in different areas, e.g., see unsteady flows [12, 13]. The solution of two-phase six-equation model is also available in the literature by upwind scheme and kinetic flux vector splitting (KFVS) scheme [10]. The results obtained by CE/SE scheme are compared with RKDG scheme [14] and central (NT) scheme [15].This article is outlined as follows. The 1-D two-phase compressible six-equation flow model is presented in Section2. A review of RKDG scheme is given in Section 3. Section 4 is devoted to numerical test problems. Conclusions are explained in Section 5.


## 2. Mathematical Model

Here, we present reduced six-equation compressible flow model which has been deduced from the seven-equation model by using asymptotic limit condition of zero velocity relaxation time in [7]. It is observed that the considered model has the ability to resolve those difficulties which were involved in five-equation model presented by [7, 16]. It is important to mention that we are considering the model without heat and mass transfer mechanisms. The one-dimensional model is given as [10].
(1)∂α1∂t+u∂α1∂x=μp1−p2,(2)∂α1ρ1∂t+∂α1ρ1u∂x=0,(3)∂α2ρ2∂t+∂α2ρ2u∂x=0,(4)∂ρu∂t+∂ρu2+α1p1+α2p2∂x=0,(5)∂α1ρ1e1∂t+∂α1ρ1e1u∂x+α1p1∂u∂x=μp1p1−p2,(6)∂α2ρ2e2∂t+∂α2ρ2e2u∂x+α1p2∂u∂x=μp2p1−p2.Here,u, ρ, and p are the velocity, the mixture density, and the mixture pressure, respectively; and volume fraction for kth phase is given by αk. The mixture energy E and mixture internal energy e are defined in the following equations
(7)E=e+12u2,(8)e=Y1e1ρ1,p+Y2e2ρ2,p.For eachkth phase, mass fraction is Yk=αρk/ρ. Furthermore, the density of the mixture of both phases is defined by ρ=αρ1+αρ2. The equation-of-state (EOS) formula is given by
(9)ek=pk+γkπkρkγk−1+qk.Here,qk, γk, and πk are constants (material) subject to the fluid. Utilizing equation (8), the mixture pressure p can be found as under
(10)p=ρe−α1γ1π1/γ1−1+α2γ2π2/γ2−1α1/γ1−1+α2/γ2−1.Here, the interfacial pressure is denoted bypI and is given by the formula
(11)p1=Z1p2+Z2p1Z1+Z2.The acoustic impedance of phasek is given by the term Zk=ρkck. For the seventh equation in the model, the mass and momentum equations are combined with the internal energy equation, and this additional equation becomes
(12)∂ρe+1/2ρu2∂t+∂uρe+1/2ρu2+α1p1+α2p2∂x=0,where ρe=α1ρ1e1+α2ρ2e2. To handle the nonconservative terms on the RHS of the model and for convergence of the schemes, nonconventional jump conditions are needed which are based on Reimann solver, see [9]. It is worth specifying that the numerical schemes proposed to solve the current model do not need these conditions. Moreover, the frozen sound speed is given as under
(13)cf2=Y1c12+Y2c22,where ckk=1,2 is sound speed of phase and is given as
(14)ck2=pk/ρk2∂ek/∂ρkpk∂ek/∂ρkpk.Now, equation (1) can be reformulated as
(15)∂α1∂t+∂α1u∂x−α1∂u∂x=0.On utilizing equations (1) to (6) along with equation (15), the model takes the following form,
(16)Wt+FWx=S∂u∂x,where
(17)W=ρuα1,FW=α2ρ2uα1ρ1uα1u,’S=α00021.The expressionsα1∂u/∂x and αkpk∂u/∂x for k=1,2 involve derivatives; therefore, these expressions are treated in the similar fashion as the flux vectors. For closure of the model equations (1) to (6) and (12), extra equations are required.Thus, in the present situation, these equations are obtained from the stiffened gas EOS (9). In the next section, we will present a brief derivation of RKDG and CE/SE schemes.


## 3. Proposed Schemes

### 3.1. RKDG Scheme

A short overview of RKDG scheme is presented in this section. The RKDG scheme is applied on1D system of hyperbolic conservation laws of the form
(18)∂tW+∂xFW=Z,t>0,x∈Ω,where ZW,x=S∂u/∂x with initial conditions
(19)Wx,0=W0x.We subdivide the domainΩ into cells Ci=xi−1/2,xi+1/2, where i=1,N. We denote the center of ith cell by xi and defined as 0. 5xi−1/2+xi+1/2. The size of ith cell is defined as △xi=xi+1/2−xi−1/2, for all i. Multiply equation (18) by a smooth function τx and integrate the resulting equation over the ith cell Ci. By integrating, we obtain
(20)∫Ci∂tWx,tτxdx−∫CiFWx,t∂xτxdx+FWxi+1/2,tτxi+1/2−FWxi−1/2,tτxi−1/2=∫CiZW,xτxdx,(21)∫CiWx,0τxdx=∫CiW0xτxdx.Consider for each timet∈0,tfinal, the solution approximation Wh of W belongs to discontinuous Galerkin finite element space
(22)Vh=v:vCi∈PkCi,i=0,1,2⋯N.ThePkC is the polynomial space of degree at most k in cell Ci. We choose the local orthogonal basis ϕil,l=0,1,2,k, named as scaled Legendre polynomials Plx, over cell Ci as defined below
(23)ϕilx=Plx−xi△xi,l=0,1,2,⋯k.Then, the approximate solutionWhx,t∈Vh is defined as
(24)Whx,t=∑l=0kWiltϕilx,forx∈Ci,where Wilt=2l+1/△xi∫CiWhx,tϕilxdx by using the orthogonal property
(25)∫Ciϕilϕimdx=0,l≠m,△xi2l+1,l=m.Now, replace a smooth functionτx by the test function ϕilx∈Vh and flux function FWxi+1/2,t by numerical flux F^ in equations (20) and (21) to obtain a numerical scheme as follows
(26)∫Ci∂tWhx,tϕilxdx−∫CiFWhx,t∂xϕilxdx+F^ϕilxi+1/2−−F^ϕilxi−1/2+=∫CiZwh,xϕilxdx,with
(27)∫CiWhx,0ϕilxdx=∫CiW0xϕilxdx.The fluxF^=┌Whxi+1/2−,t,Whxi+1/2+,t is any monotone numerical flux. Here, we use the LFF as a monotone numerical flux which is defined below
(28)┌Wi+1/2−,Wi+1/2+=12FWi+1/2−+FWi+1/2+−θWi+1/2+−Wi+1/2−,where θ=maxWFW and Whxi+1/2±,t are the limited values of approximate solution Whx,t at the cell interface xi+1/2. By using above definitions, equations (26) and (27) simplify to
(29)dWiltdt=−2l+1△xi∫CiFWhx,t∂xϕilxdx+F^ϕilxi+1/2−−F^ϕilxi−1/2++∫CiZW,xτxdx,(30)Wil0=2l+1△xi∫CiW0xϕilxdx.Next, to obtain nonoscillatory discontinuity propagation and high order accuracy by RKDG method, the limiter procedure is divided into two parts. First, identify the cells, named as “troubled cells,” that may need the limiting procedure. Second, reconstruct the polynomial solutions in these troubled cells by using WENO reconstruction.Here, the troubled cells are identified by total variation bounded (TVB) limiter, which is described as(31)Whxi+1/2+=Wi0+W¯,Whxi−1/2+=Wi0−W¯¯,where
(32)W¯=∑l=1kWilϕilxi+1/2,W¯¯=−∑l=1kWilϕilxi−1/2.Now,W¯ and W¯¯ are modified by the modified minmod (MM) function, as follows
(33)W¯=MMW¯,Wi+10−Wi0,Wi0−Wi−10,(34)W¯¯=MMW¯¯,Wi+10−Wi0,Wi0−Wi−10,where MM is defined as
(35)MMa1,a2,a3=a1,ifa1≤M△x2,Ma1,a2,a3,otherwise,where M is a positive constant and the minmod function M is defined as
(36)Ma1,a2,a3=S·min1≤i≤3ai,ifsigna1=signa2=signa3=S,0,otherwise.Using explicit, nonlinearly stable high order Runge-Kutta time discretization. [Shu and Osher, JCP, 1988].The semidiscrete scheme (26) is written as
(37)Wt=LW,is discretized in time by a nonlinearly stable Runge-Kutta time discretization, e.g., the third order version.


### 3.2. The CE/SE Scheme

In this section, we give a brief derivation ofCE/SE scheme [1, 11].Let the coordinates of a 2-D Euclidean space (E) be x0=t and x1=x. The equivalent integral form can be obtained by employing the Gauss-divergence theorem, on equation (15).
(38)∮SVhi·dS=∫VSidV,i=1−6,in which (a) hi=wi,fiT,i=1,2,3,4, i.e., for all i=1,2,3,4,wi and fi be the elements of the vector hi in the t− and x−directions, respectively, and (b) def dS=dσn in which dσ represents the area whereas n being the typical unit outward vector of a surface component on SV. Over the space-time region CE (conservation element), equation (38) is implemented. The actual numerical combination uses solution components (SEs) in a discreet manner.InE2 space, the arrangement of work focuses j,k is represented by Ω1 where k=0,1/2,1,3/2,⋯ and for each k,j=0,±1/2,±1,±3/2,⋯. There is a SE (solution element) corresponding to each j,k, as depicted in Figure 1 (dashed curve). It includes both a horizontal and a vertical line and the immediate neighborhood.Figure 1
Staggered space-time grid.
(a)(b)(c)The correct size of immediate vicinity does not make a difference. For anyx,t∈SEj,k,wix,t,fi and hi are estimated by wi∗x,t;j,n, fi∗x,t;j,k, and hi∗x,t;j,k separately as takes the following form
(39)wi∗x,t;j,kd=efwijk+witjkt−tk+wixjkx−xj,(40)fi∗x,t;j,kd=effljk+fitjkt−tk+fixjkx−xj.Moreover,(41)τi∗x,t;j,kd=efτwi∗,wix∗,Using chain rule, we acquire(42)fixjk=∑n=14fi,njkwnxjkdef,(43)fitjk=∑n=14fi,njkwntjkdef,where
(44)def∂fifi,n=∂wn¯,i,n=1−6.Here,fi,n,i,k=1−6, where i and k, respectively, represent the column and row indices. Note that wijk,wixjk, and witjk are constants in SE j,k. The numerical analogs of these values are the following wi,∂wi/∂x, and ∂wi/∂t at xj,tk, respectiyely. Since hid=efwi,fiT, therefore,
(45)hi∗x,t:j,kdef=wi∗x,t:j,k,fi∗x,t:j,kT.Using equations (39), (40), (41), and (45), we get
(46)wljk=wij±1/2k−1/2∓△x4wixj±1/2k−1/2+wixjk∓△t△xfij±1/2k−1/2−fijk∓△t24△xfitj±1/2k−1/2+fitjk±Slj±1/2k−1/2−Sljk.The summation gives(47)wiljk=12wij−1/2k−1/2+wij+1/2k−1/2+Qij−1/2k−1/2−Qilj+1/2k−1/2+Slj+1/2k−1/2−Slj−1/2k−1/2,where
(48)Qij±1/2k−1/2=△x4wixj±1/2k−1/2+△t△xfij±1/2k−1/2+△t24△xfitj±1/2k−1/2.Clearly from the source term,S2=S3=S4=0. Moreover,
(49)S1j+1/2k−1/2−S1j−1/2k−1/2=△t2△xα1j−1/2k−1/2+α1j+1/2k−1/2,(50)S5j+1/2k−1/2−S5j−1/2k−1/2=△t2△xα1p1j−1/2k−1/2+α1p1j+1/2k−1/2.(51)S6j+1/2k−1/2−S6j−1/2k−1/2=△t2△xα2p2j−1/2k−1/2+α2p2j+1/2k−1/2.For more details, the reader is referred to [1, 11].


## 3.1. RKDG Scheme

A short overview of RKDG scheme is presented in this section. The RKDG scheme is applied on1D system of hyperbolic conservation laws of the form
(18)∂tW+∂xFW=Z,t>0,x∈Ω,where ZW,x=S∂u/∂x with initial conditions
(19)Wx,0=W0x.We subdivide the domainΩ into cells Ci=xi−1/2,xi+1/2, where i=1,N. We denote the center of ith cell by xi and defined as 0. 5xi−1/2+xi+1/2. The size of ith cell is defined as △xi=xi+1/2−xi−1/2, for all i. Multiply equation (18) by a smooth function τx and integrate the resulting equation over the ith cell Ci. By integrating, we obtain
(20)∫Ci∂tWx,tτxdx−∫CiFWx,t∂xτxdx+FWxi+1/2,tτxi+1/2−FWxi−1/2,tτxi−1/2=∫CiZW,xτxdx,(21)∫CiWx,0τxdx=∫CiW0xτxdx.Consider for each timet∈0,tfinal, the solution approximation Wh of W belongs to discontinuous Galerkin finite element space
(22)Vh=v:vCi∈PkCi,i=0,1,2⋯N.ThePkC is the polynomial space of degree at most k in cell Ci. We choose the local orthogonal basis ϕil,l=0,1,2,k, named as scaled Legendre polynomials Plx, over cell Ci as defined below
(23)ϕilx=Plx−xi△xi,l=0,1,2,⋯k.Then, the approximate solutionWhx,t∈Vh is defined as
(24)Whx,t=∑l=0kWiltϕilx,forx∈Ci,where Wilt=2l+1/△xi∫CiWhx,tϕilxdx by using the orthogonal property
(25)∫Ciϕilϕimdx=0,l≠m,△xi2l+1,l=m.Now, replace a smooth functionτx by the test function ϕilx∈Vh and flux function FWxi+1/2,t by numerical flux F^ in equations (20) and (21) to obtain a numerical scheme as follows
(26)∫Ci∂tWhx,tϕilxdx−∫CiFWhx,t∂xϕilxdx+F^ϕilxi+1/2−−F^ϕilxi−1/2+=∫CiZwh,xϕilxdx,with
(27)∫CiWhx,0ϕilxdx=∫CiW0xϕilxdx.The fluxF^=┌Whxi+1/2−,t,Whxi+1/2+,t is any monotone numerical flux. Here, we use the LFF as a monotone numerical flux which is defined below
(28)┌Wi+1/2−,Wi+1/2+=12FWi+1/2−+FWi+1/2+−θWi+1/2+−Wi+1/2−,where θ=maxWFW and Whxi+1/2±,t are the limited values of approximate solution Whx,t at the cell interface xi+1/2. By using above definitions, equations (26) and (27) simplify to
(29)dWiltdt=−2l+1△xi∫CiFWhx,t∂xϕilxdx+F^ϕilxi+1/2−−F^ϕilxi−1/2++∫CiZW,xτxdx,(30)Wil0=2l+1△xi∫CiW0xϕilxdx.Next, to obtain nonoscillatory discontinuity propagation and high order accuracy by RKDG method, the limiter procedure is divided into two parts. First, identify the cells, named as “troubled cells,” that may need the limiting procedure. Second, reconstruct the polynomial solutions in these troubled cells by using WENO reconstruction.Here, the troubled cells are identified by total variation bounded (TVB) limiter, which is described as(31)Whxi+1/2+=Wi0+W¯,Whxi−1/2+=Wi0−W¯¯,where
(32)W¯=∑l=1kWilϕilxi+1/2,W¯¯=−∑l=1kWilϕilxi−1/2.Now,W¯ and W¯¯ are modified by the modified minmod (MM) function, as follows
(33)W¯=MMW¯,Wi+10−Wi0,Wi0−Wi−10,(34)W¯¯=MMW¯¯,Wi+10−Wi0,Wi0−Wi−10,where MM is defined as
(35)MMa1,a2,a3=a1,ifa1≤M△x2,Ma1,a2,a3,otherwise,where M is a positive constant and the minmod function M is defined as
(36)Ma1,a2,a3=S·min1≤i≤3ai,ifsigna1=signa2=signa3=S,0,otherwise.Using explicit, nonlinearly stable high order Runge-Kutta time discretization. [Shu and Osher, JCP, 1988].The semidiscrete scheme (26) is written as
(37)Wt=LW,is discretized in time by a nonlinearly stable Runge-Kutta time discretization, e.g., the third order version.


## 3.2. The CE/SE Scheme

In this section, we give a brief derivation ofCE/SE scheme [1, 11].Let the coordinates of a 2-D Euclidean space (E) be x0=t and x1=x. The equivalent integral form can be obtained by employing the Gauss-divergence theorem, on equation (15).
(38)∮SVhi·dS=∫VSidV,i=1−6,in which (a) hi=wi,fiT,i=1,2,3,4, i.e., for all i=1,2,3,4,wi and fi be the elements of the vector hi in the t− and x−directions, respectively, and (b) def dS=dσn in which dσ represents the area whereas n being the typical unit outward vector of a surface component on SV. Over the space-time region CE (conservation element), equation (38) is implemented. The actual numerical combination uses solution components (SEs) in a discreet manner.InE2 space, the arrangement of work focuses j,k is represented by Ω1 where k=0,1/2,1,3/2,⋯ and for each k,j=0,±1/2,±1,±3/2,⋯. There is a SE (solution element) corresponding to each j,k, as depicted in Figure 1 (dashed curve). It includes both a horizontal and a vertical line and the immediate neighborhood.Figure 1
Staggered space-time grid.
(a)(b)(c)The correct size of immediate vicinity does not make a difference. For anyx,t∈SEj,k,wix,t,fi and hi are estimated by wi∗x,t;j,n, fi∗x,t;j,k, and hi∗x,t;j,k separately as takes the following form
(39)wi∗x,t;j,kd=efwijk+witjkt−tk+wixjkx−xj,(40)fi∗x,t;j,kd=effljk+fitjkt−tk+fixjkx−xj.Moreover,(41)τi∗x,t;j,kd=efτwi∗,wix∗,Using chain rule, we acquire(42)fixjk=∑n=14fi,njkwnxjkdef,(43)fitjk=∑n=14fi,njkwntjkdef,where
(44)def∂fifi,n=∂wn¯,i,n=1−6.Here,fi,n,i,k=1−6, where i and k, respectively, represent the column and row indices. Note that wijk,wixjk, and witjk are constants in SE j,k. The numerical analogs of these values are the following wi,∂wi/∂x, and ∂wi/∂t at xj,tk, respectiyely. Since hid=efwi,fiT, therefore,
(45)hi∗x,t:j,kdef=wi∗x,t:j,k,fi∗x,t:j,kT.Using equations (39), (40), (41), and (45), we get
(46)wljk=wij±1/2k−1/2∓△x4wixj±1/2k−1/2+wixjk∓△t△xfij±1/2k−1/2−fijk∓△t24△xfitj±1/2k−1/2+fitjk±Slj±1/2k−1/2−Sljk.The summation gives(47)wiljk=12wij−1/2k−1/2+wij+1/2k−1/2+Qij−1/2k−1/2−Qilj+1/2k−1/2+Slj+1/2k−1/2−Slj−1/2k−1/2,where
(48)Qij±1/2k−1/2=△x4wixj±1/2k−1/2+△t△xfij±1/2k−1/2+△t24△xfitj±1/2k−1/2.Clearly from the source term,S2=S3=S4=0. Moreover,
(49)S1j+1/2k−1/2−S1j−1/2k−1/2=△t2△xα1j−1/2k−1/2+α1j+1/2k−1/2,(50)S5j+1/2k−1/2−S5j−1/2k−1/2=△t2△xα1p1j−1/2k−1/2+α1p1j+1/2k−1/2.(51)S6j+1/2k−1/2−S6j−1/2k−1/2=△t2△xα2p2j−1/2k−1/2+α2p2j+1/2k−1/2.For more details, the reader is referred to [1, 11].


## 4. Numerical Test Problems

In the present section, theCE/SE scheme is used to study the single velocity six-equation model along with a pressure relaxation procedure to accommodate the nonconservative terms in the model. It has been mentioned earlier that two-phase flow models inherit serious computational difficulties. Some one-dimensional test problems carried out in this section and their results are compared with each other, to authenticate the results of two proposed schemes. To further authenticate the efficiency of suggested schemes, we have used the RKDG [14] and central (NT) [15] schemes.Problem 1 No reflection test problems.
This test case was also studied in [8] by implementing KFVS scheme. The Riemann left and right data about interface x=0.5 is given as
(52)ρ,p,u,α=3.1748,100,1,9.435,ifx≤0.5,1,1,0,0,otherwise.The problem is simulated at 500 mesh cells, andt=0.02 is the final simulation time. Here, γL=1.667 and γR=1.2 are values of specific heat ratios. Furthermore, πL=πR, and CFL=0.32. Because of significant pressure shifts at the interface, it is one of the tough test problems for a numerical scheme. Choosing the velocity and jumping pressure over the shock makes it difficult to create reflection wave. Figure 2 displays the results obtained from two schemes. Clearly, wiggles can be observed in pressure and velocity plots from all the schemes. These wiggles disappear on refined mesh. Similar type of results is reported in Kreeft [17]. Furthermore, we have also computed the L1 error of CE/SE and RKDG schemes. The errors are given in Table 1. We can clearly observe that CE/SE has less error as compared to RKDG scheme.Figure 2
Problem1: results for alpha α, density (ρ), velocity (u), and pressure (p).Table 1
Comparison ofL1-errors in the schemes.
NRKDGCE/SE400.43200.5850800.89970.96391600.98331.05073200.97210.98356401.18181.2231TheL1-errors for CE/SE, KFVS, and central (NT) schemes at different grid points are listed in Table 1. We can clearly observe from Table 1 that CE/SE scheme has less error as compared to KFVS and central (NT) schemes.Problem 2 Sode’s test problem.
The problem considered is like the Sod problem of single-phase gas dynamics, which is simulated by [10] as well. A thin membrane fixed at x=1/2 separates the two gases which are at rest initially. The gases on both sides of the interface mixes with each other as soon as membrane is removed. The initial state of problem is as under
(53)ρ,p,u,α=310,10,0,1,ifx≤12,0.125,0.1,0,0,otherwise.The specific heat ratioγl=1.4 is on the of interface, and γr=1.6 is on the right of interface. The results at t=0.015 are presented in Figure 3. From the simulated results, a left-going rarefaction wave can be seen. Also, a right-moving two-fluid interface and right-going shock wave can be observed. In Figure 3, the results by proposed schemes are shown at 400 cells. For this comparison, the reference solution is simulated for CE/SE scheme at 2000 meshes. Both the schemes are in close agreement capturing the shocks effectively.Problem 3 Two fluid mixture problem.
This two-fluid mixture problem has been considered in [10] and was simulated by five-equation model. The values of the other physical variables are given as
(54)ρ,p,u,α=2.0,1000,0,1,ifx≤12,1,0.01,0,0,otherwise.Figure 3
Results for alpha (a), density (b), velocity (c), and pressure (d).
(a)(b)(c)(d)In this case study,γL=1.2 and γR=1.4. In the figure, one can see a right moving shock wave, a left going rarefaction wave, and a contact discontinuity. The right going shock hits the interface at x=0.5. The shock continues to move towards right, and a rarefaction wave is created which is moving towards left. For simulation of the results, we take 400 cells, and t=0.012 is the final time. The results are presented in Figure 4. The results reveal that two schemes give comparable results. However, the results of CE/SE scheme are better than the other schemes.Problem 4 Comparison test problem.
The present test of water-air mixture is presented here withCE/SE and RKDG schemes, respectively. The initial values of the problem are [10]
(55)ρ,p,u,α=2.0,1000,0,1,ifx≤12,1,0.01,0,0,otherwise,where γR=4.4,γL=1.4,πL=6x108,CFL=0.4, and α=1. The results at t=200μs are depicted in Figure 5. We can observe that all the schemes are in close agreement with each other. However, among all CE/SE performance is better.Figure 4
Results of Problem3.Figure 5
Results of Problem4.


## 5. Conclusions

In this article, reduced six-equation flow model was numerically investigated by the space-timeCE/SE and RKDG schemes. The numerical complexities related to the model include accurate discretization of source term. The proposed methods have the ability to resolve sharp discontinuous profiles. The numerical test problems show the performance of the suggested method to handle a wide range of flow conditions. The results of suggested schemes are compared with those obtained from the central (NT) scheme and exact solutions. A good agreement was observed among the results of CE/SE and RKDG schemes. However, it was found that our suggested CE/SE scheme effectively resolves sharp discontinuities better than central (NT). The CE/SE scheme generally proposes a good strategy for solving both nonconservative and conservative system of equations.



---

*Source: 1697181-2022-07-15.xml*