A staggered-projection Godunov-type method for the Baer-Nunziato two-phase model
aa r X i v : . [ phy s i c s . c o m p - ph ] J u l A staggered-projection Godunov-type method forthe Baer-Nunziato two-phase model
Xin Lei a , Jiequan Li b,c,d, ∗ a School of Science, China University of Geosciences, Beijing, 100083, P. R. China b Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing, 100088,P. R. China c Center for Applied Physics and Technology, Peking University, 100871, P. R. China d State Key Laboratory for Turbulence Research and Complex System, Peking University, 100871, P. R. China
Abstract
When describing the deflagration-to-detonation transition in solid granular explosives mixed withgaseous products of combustion, a well-developed two-phase mixture model is the compressible Baer-Nunziato (BN) model, containing solid and gas phases. If this model is numerically simulated by aconservative Godunov-type scheme, spurious oscillations are likely to generate from porosity interfaces,which may result from the average process of conservative variables that violates the continuity of Rie-mann invariants across porosity interfaces. In order to suppress the oscillations, this paper proposesa staggered-projection Godunov-type scheme over a fixed gas-solid staggered grid, by enforcing thatsolid contacts with porosity jumps are always inside gaseous grid cells and other discontinuities appearat gaseous cell interfaces. This scheme is based on a standard Godunov scheme for the Baer-Nunziatomodel on gaseous cells and guarantees the continuity of the Riemann invariants associated with thesolid contact discontinuities across porosity jumps. While porosity interfaces are moving, a projectionprocess fully takes into account the continuity of associated Riemann invariants and ensure that poros-ity jumps remain inside gaseous cells. This staggered-projection Godunov-type scheme is well-balancedwith good numerical performance not only on suppressing spurious oscillations near porosity interfacesbut also capturing strong discontinuities such as shocks.
Keywords:
Compressible two-phase flow, Baer-Nunziato (BN) model, Godunov-type scheme,non-conservative products, Riemann invariants
1. Introduction
In the two-phase flows containing solid and gas phases, dispersed solid particles in the gas can beconsidered as a continuous fluid that penetrates the gas phase. Each phase in the gas-solid two-phase ✩ The second author is supported by NSFC (nos. 11771054, 91852207), the Sino-German Research Group Project,National key project(GJXM92579) and Foundation of LCP. ∗ Corresponding author
Email addresses: [email protected] (Xin Lei), [email protected] (Jiequan Li)
Preprint submitted to Journal of L A TEX Templates July 29, 2020 ow is usually non-equilibrium and has individual state variables, including the density, velocity, andpressure. In 1986, Baer and Nunziato proposed a two-velocity two-pressure model for the compressibletwo-phase flow [1], commonly known as the Baer-Nunziato (BN) model. This model provides a quan-titative analysis of the deflagration-to-detonation transition in porous granular explosives. Neglectingthe non-differential source terms due to combustion, drag, heat transfer and chemical reaction in thecomplete BN model, the governing equations remain balance laws of mass, momentum and energyfor each phase, coupled with an evolution equation for the porosity, i.e., the gaseous volume fraction[2]. This simplified model is usually called the homogeneous BN model, and it consists of a hyperbolicsystem of compressible Euler equations plus nozzling terms of non-conservative products. The studyof the homogeneous BN model plays an important role in simulating the complete BN model; however,the homogeneous BN model cannot be written in conervative form so that there is no well-acceptableunderstanding of the non-conservative products when dealing with strong discontinuities, large defor-mation of interfaces, mixing of the two phases and interphase exchange [3]. The goal of this paper isto construct an efficient numerical scheme to compute the homogeneous BN model with good fidelity,through a careful analysis of the behavior of solid contacts at porosity interfaces.Since shocks may exist in the flow, a numerical scheme should be conservative and thus Godunov-type schemes [4] are preferable, for which the associated Riemann problems are solved numerically atcell interfaces to resolve wave structures properly. In the Eulerian framework, many Godunov-typemethods computing the homogeneous BN two-phase flow model have been proposed, including theoperator splitting method [5, 6], the unsplit Roe-type wave-propagation method [6, 7, 8, 9] and thepath-conservative scheme [10, 11]. However, there are some flaws that need to be fixed. For exam-ple, Lowe showed in [6] the performance of various conservative shock-capturing schemes, exhibitingspurious oscillations (or visible errors) around porosity interfaces. In [9], Karni proposed a hybrid al-gorithm composed of a non-conservative Roe-type scheme across the porosity jump and a conservativeRoe-type method away from the porosity jump. Here, the non-conservative Roe-type scheme utilizesthe Riemann invariants across the porosity jump to avoid spurious oscillations. However, the hybridscheme is not fully conservative and needs to track the porosity interface. In addition to Godunov-typemethods, a non-conservative finite-volume approach that uses residuals instead of fluxes was designed[12], and it is compatible with local conservation.As far as a Godunov-type method is applied to a non-conservative hyperbolic system, a key pointis to compute numerical fluxes by solving the associated Riemann problem at each cell interface. Ifthe Riemann problem contains a porosity jump, the solution consists of complicated wave structureson account of the nozzlling terms, especially because of the resonance phenomenon meaning waves ofdifferent families coincide [13, 14]. In [15], Andrianov and Warnecke analyzed various wave structuresof the homogeneous BN model and designed an exact Riemann solver for given wave patterns, which2s called the inverse Riemann solver . Schwendeman et al. [16] proposed an exact Riemann solverfor solving the Riemann problem directly. In the subsequent paper [17], a different exact Riemannsolver, also dealing with resonance, was developed. Tokareva and Toro [18] designed a HLLC-typeapproximate Riemann solver to treat the solid contact. However, for some initial data, the Riemannsolution may be not unique or even do not exist [15, 17]. This paper is mainly restricted to the casethat the Riemann solution exists uniquely.Another key point is about the numerical approximation of non-conservative products. As intro-duced in [19, 20, 15, 16] and summarized in [17], there are various numerical methods for discretizingnozzlling terms and corresponding numerical results satisfy the Abgrall criterion (or free-streamingcondition) generalized from [21], which requires uniform velocity and pressure be preserved for two-phase flows. Nevertheless, the improper numerical integration of non-conservative products acrossporosity jumps could cause spurious oscillations, which will be analyzed in this paper. Across theporosity interface, the Riemann invariants for a 0-characteristic field keep constant across the asso-ciated solid contact. However, this is not the case numerically. It is observed that these Riemanninvariants cannot always keep constant with visible errors and spurious oscillations when the conser-vative scheme is adopted for an isolated moving porosity interface. Hence, in order to eliminate theoscillations, the conservation requirement of numerical solutions across porosity interfaces should beappropriately liberated. Re and Abgrall built a non-conservative pressure-based method for weaklycompressible BN-type model [22, 23]. This method utilizes a staggered description of the flow vari-ables to avoid spurious pressure oscillations. Furthermore, in this paper, we concentrate on gettingoscillation-free solutions of general compressible BN flows that contain shock waves, contact disconti-nuities, and strong rarefactions. The robust approximation of non-conservative products also dependson the positivity and well-balanced property. Recent work in this aspect includes a entropy-satisfyingscheme in [24] and a well-balanced scheme in [25] proven to preserve positive phase volume fractions,densities and internal energies.With the above thinkings in mind, a staggered-projection Godunov-type scheme is designed tosuppress spurious oscillations near porosity interfaces. For such a scheme with piecewise constantinitial data, solid contacts, which are associated with the 0-characteristic field, are assumed to alwaysexist inside fixed grid cells, named as “gaseous cells” , whereas other waves in the Riemann solutionare generated from gaseous cell interfaces. We take these solid contacts as interfaces of “solid cells” ,which are staggered with the gaseous cells. In each gaseous cell, the standard Godunov scheme is usedto update all conservative variables, and the states on both sides of the solid contact are computedin terms of the same family of Riemann invariants associated with the 0-characteristic field. Thenumerical integration of nozzling terms are computed precisely so that this scheme is well-balanced.Moreover, since there is no porosity jump at any gaseous cell interface, it becomes routine to solve the3iemann problem at each cell interface for decoupled two phases. Technically, in order to ensure thatmoving porosity jumps are always located inside fixed gaseous cells, the porosity discontinuities areprojected back to their original locations at each time step. The principle of this projection process is tokeep constant Riemann invariants for the 0-characteristic field, and thus the total energy conservationhas to be given up across porosity jumps for avoiding the contradiction with the former. Afterwards,the staggered-projection Godunov-type scheme is extended to space-time second-order accuracy usingthe generalized Riemann problem (GRP) method for hyperbolic balance laws [26], and extended to amulti-dimensional regular grid by the dimensional splitting method [27] accordingly.This paper is organized as follows. In Section 2, we summarize the characteristic analysis of thehomogeneous BN model in [1, 28] and the properties of the solid contact in the exact Riemann solutionare reviewed. A brief introduction of standard Godunov methods for the homogeneous BN model isdescribed in Section 3 with detailed analysis of visible numerical errors across porosity interfaces. Thestaggered-projection Godunov-type scheme with the second-order GRP version and two-dimensionalextension are presented in Section 4. In Section 5, several examples are presented to demonstrate theperformance. Finally, conclusions are made in Section 6. The technical representation of BN modeland approximate solutions of the nonlinear equations related to Riemann invariants is put in Appendix.
2. The Riemann problem for homogeneous BN model
In this section, we provide some properties of the homogeneous BN model and discuss the Riemannproblem for the later design of Godunov-type schemes. In particular, we investigate the property ofsolid contacts associated with porosity jumps.
A detailed characteristic analysis of the compressible BN two-phase model was presented in [28],and those to be used in the following sections are summarized below. For this model, the gaseousproducts and the dispersed solid particles are considered as two different phases, referring to the gasphase and the solid phase, respectively. In the one-dimensional (1-D) case, the state of each phase isdetermined by a volume fraction α k , density ρ k , velocity u k , pressure p k , and a total specific energy E k , where k = s and g represent the solid phase and the gas phase, respectively. The gaseous volumefraction α g measures the porosity of the solid phase and satisfies the saturation constraint α g = 1 − α s .Ignoring the non-differential source terms, the system of governing equations for the homogeneous BNmodel is written as u t + f ( u ) x = h ( u )( α s ) x , (1)4ith u = α s α s ρ s α s ρ s u s α s ρ s E s α g ρ g α g ρ g u g α g ρ g E g , f = α s ρ s u s α s ρ s u s + α s p s α s u s ( ρ s E s + p s ) α g ρ g u g α g ρ g u g + α g p g α g u g ( ρ g E g + p g ) , h = − u s p g p g u s − p g − p g u s , where E k = e k + u k , e k = e k ( ρ k , p k ) is the specific internal energy specified by the equations ofstate (EOS) for each phase. The non-conservative product h ( u )( α s ) x is called the nozzling termsince it resembles the non-conservative term of the quasi-1-D compressible duct flow model (12) ina converging-diverging nozzle that we will discuss briefly later on. The advection equation for thevolume fraction, ( α s ) t + u s ( α s ) x = 0, represents that the porosity is passively advected with the localsolid velocity u s . We assume that each isolated phase satisfies the thermodynamic (Gibbs) relation,d e k = T k d η k + p k ρ k d ρ k , (2)where η k and T k are the entropy and temperature for the phase k , respectively.For smooth solutions, system (1) can be expressed in terms of primitive variables v t + A ( v ) v x = , with v = α s ρ s u s p s ρ g u g p g , A = u s u s ρ s p s − p g ) / ( α s ρ s ) 0 u s /ρ s ρ s c s u s ρ g ( u s − u g ) /α g u g ρ g
00 0 0 0 0 u g /ρ g ρ g c g ( u s − u g ) /α g ρ g c g u g , where c k is the sound speed related to the entropy η k through the thermodynamic relation (2), c k = s ∂p k ( ρ k , η k ) ∂ρ k , the matrix A has seven real eigenvalues λ = u s ,λ ,s = u s − c s , λ ,s = u s , λ ,s = u s + c s ,λ ,g = u g − c g , λ ,g = u g , λ ,g = u g + c g . λ i,s -fields, i = 1 , ,
3, are the same as the threeeigen-fields of the one-dimensional Euler equations for the solid phase, and so are the λ i,g -characteristicfields for the gas phase. Linearly degenerate fields λ , λ ,s and λ ,g define contact discontinuities (solidcontacts and gas contacts); other genuinely nonlinear fields define nonlinear waves (rarefaction wavesand shocks). Note that λ ,s < λ ,s < λ ,s and λ ,g < λ ,g < λ ,g . However, there are no orderingsamong λ i,s and λ i,g , i = 1 , ,
3, which leads to the occurrence of resonance. Typically, under the soniccondition u s = u g ± c g , a solid contact is located inside a rarefaction wave or coincides with a shockof the gas phase and thus the resonance occurs. Only as the following conditions hold,( u g − u s ) = c g , α s = 0 , α g = 0 , the coefficient matrix A has a complete set of eigenvectors and the system (1) is hyperbolic. Note that λ = λ ,s = u s . The resonance may occur, which is analogous to the resonant case of the duct flowthat we will discuss later on. The Riemann problem is the building block of Godunov-type schemes. Assume, for the 1-D ho-mogeneous BN model (1), that piecewise constant data at initial time t = t n that can be shifted to t = 0, u ( x, t = 0) = u L , x < x , u R , x > x , where u L and u R are constant states on both sides of certain position x = x . The solution of thisRiemann problem has self-similarity, u ( x, t ) = v ( ξ ) , ξ = x − x t , t > , denoted as RP ( ξ ; u L , u R ) and this solution consists of shock waves S k , rarefaction waves R k , andcontact discontinuities C k ( k = s, g ) for each phase. For the time being, we leave aside the discussionon the coalescence of these waves for different phases to the next section, and pay attention to solidcontacts λ = λ ,s . A solid contact is associated with λ -field, and the corresponding five Riemanninvariants are u s , η g ,Q := α g ρ g ( u g − u s ) ,P := α s p s + α g p g + α g ρ g ( u g − u s ) ,H := h g + ( u g − u s ) , where h g = e g + p g /ρ g is the enthalpy of the gas phase. Note that the state variables of the gas phaseand the solid pressure p s do not remain constant across the solid contact, except the situation that6he porosity on both sides of the contact is identical, i.e., ( α s ) L = ( α s ) R . Away from the solid contact,the porosity is constant and the nozzling term is zero, thus the governing equations (1) are decoupledand reduce to the Euler equations for the two individual phases. The state variables of one phase areconstant across the waves of the other phase.Intermediate states of the solid phase Intermediate states of the gas phase Figure 1: The wave configuration to the Riemann problem for a subsonic configuration λ ,g < λ < λ ,g . Upon the position of the solid contact relative to the gas phase, wave configurations of the Riemannproblem are classified into two categories: (i) the solid contact is outside of gas waves , λ < λ ,g , or λ > λ ,g ; (ii) the solid contact is inside the gas waves , λ ,g < λ < λ ,g or λ ,g < λ < λ .g . Forcertain initial data, there may be more than one Riemann solution from different configurations; andfor certain other initial data, the Riemann solution does not exist [15]. This ill-posedness issue mayreflect the invalidity of the Riemann problem for the homogeneous BN model, and it has not beeneffectively resolved [17]. Therefore, in this paper, we are chiefly concerned with the situation that theRiemann solution exists uniquely, and try to get around this issue by solving more simplified Riemannproblems. For simplicity, we take the configuration (ii) as an example for analysis, as shown in Figure1. In this configuration, the gas phase is subsonic relative to the solid phase, i.e., ( u g − u s ) < c g . Itholds in the context of granular explosives. For the solid phase in this configuration, the intermediatestates on the left and right sides of the solid contact are labelled by subscripts 1 and 2, respectively.The gas phase is labelled similarly, in addition to indicating the intermediate state between the gascontact and the solid contact by the subscript 0.A schematic diagram for the complete wave configuration above λ ,g < λ s < λ ,g is shown in Figure2. There are two rays ξ = ξ L and ξ = ξ R , ξ L = ( u s ) ∗ − ǫ , ξ R = ( u s ) ∗ + ǫ and ǫ >
0, bounding asectorial region S ( ξ L , ξ R ) by recalling that the similarity variable ξ = ( x − x ) /t and the solid velocity( u s ) = ( u s ) =: u ∗ s . The solution can only be discontinuous across the solid contact ξ = ( u s ) ∗ . Theporosity is constant outside the region S ( ξ L , ξ R ), namely ( α s ) x ≡ igure 2: The solution in a small sector covering a solid contact propagating to the right. effect only along the solid contact. Since the Riemann solution u ( x, t ) is self-similar locally, we makethe change of variables ( x, t ) → ( ξ, t ) to recast system (1) in the form − ξ v ξ + f ( v ) ξ = h ( v )( α s ) ξ . (3)Integrating system (3) from ξ L to ξ R , we obtain Z ξ R ξ L [ − ξ v ξ + f ( v ) ξ ] d ξ = Z ξ R ξ L h ( v )( α s ) ξ d ξ. Specified to the solid phase, we have Z ξ R ξ L [ − ξ ( α s ρ s ) ξ + ( α s ρ s u s ) ξ ] d ξ = 0 , (4a) Z ξ R ξ L (cid:2) − ξ ( α s ρ s u s ) ξ + ( α s ρ s u s + α s p s ) ξ (cid:3) d ξ = Z ξ R ξ L p g ( α s ) ξ d ξ, (4b) Z ξ R ξ L [ − ξ ( α s ρ s E s ) ξ + ( α s ρ s u s E s + α s u s p s ) ξ ] d ξ = Z ξ R ξ L p g u s ( α s ) ξ d ξ, (4c)and Z ξ R ξ L − u s ( α s ) ξ d ξ = − u ∗ s [( α s ) R − ( α s ) L ] . Since the Riemann invariant u s is a constant across the solid contact, (4a) always holds, and theequations (4b) and (4c) hold under the condition Z ξ R ξ L p g ( α s ) ξ d ξ = Z ξ R ξ L ( α s p s ) ξ d ξ = ( α s ) R ( p s ) − ( α s ) L ( p s ) . (5)This provides a way to discretize the nozzling term.8 . Analysis of spurious oscillations by standard Godunov schemes In this section we review the standard Godunov method for the homogeneous BN model and analyzespurious oscillations in the solution. For the 1-D model, the computational domain [0 , L ] are dividedinto M fixed grid cells I i = [ x i − , x i + ] , i = 1 , , . . . , M , ∆ x = x i + − x i − = L/M , the cell interface x i + = i ∆ x , and the cell center x i = ( i − / x , respectively. The Godunov scheme assumes theinitial data to be piece-wise constant at time t = t n , u ni = 1∆ x Z x i + 12 x i − u ( x, t n )d x. The next step is to compute numerical fluxes along all cell interface x = x i + by solving the localRiemann problem RP (cid:0) ξ ; u ni , u ni +1 (cid:1) , ξ = x − x i + 12 t − t n . For a time increment ∆ t restricted by the CFLcondition, we can evolve the solution to the next time t n +1 = t n +∆ t by the finite-volume discretization u n +1 i − u ni + ∆ t ∆ x (cid:16) f ni + − f ni − (cid:17) = ∆ t ∆ x S ni , (6)where the numerical flux along the cell interface x = x i + is f ni + = 1∆ t Z t n +1 t n f (cid:16) u ( x i + , t ) (cid:17) d t = f (cid:16) u ni + (cid:17) , with the Riemann solution u ni + = RP (0; u ni , u ni +1 ) and the time average of the numerical integral forthe nozzling term S ni = 1∆ t Z t n +1 t n Z x i + 12 x i − h ( u )( α s ) x d x d t. The numerical fluxes can be evaluated by exact Riemann solvers [16, 17] or approximate Riemannsolvers, such as a Roe solver [9], a Rusanovtype solver [29] or a more precise HLLC solver [18]. Inaddition, the non-conservative product in (1) entails extra difficulties for the numerical simulation. Assummarized in [17], there are several approaches to approximate the integral of the nozzling term S ni ,two of which are listed below.The first follows the finite-volume method [20] and the integral of the nozzling term is approximatedas S ni ≃ h ( u ni ) (cid:16) ( α s ) ni + − ( α s ) ni − (cid:17) , (7)where ( α s ) ni + is the exact value of α along the cell interface. This approach assures that the numericalsolutions satisfy the free-streaming condition proposed by Abgrall [21]. That is, the uniformity ofpressure and velocity remains during their time evolution.The other approach was first proposed in [16], which precisely evaluates the integral of the nozzlingterm and seems more reasonable. According to the interpretation of the Riemann solution in Section2.2, ( α s ) x ≡ u s > igure 3: Solid contacts in the local Riemannian problems at cell interfaces. the computational domain under consideration. Then the solid contact propagates into a cell fromits left interface. As shown in Figure 3, the integral of the nozzling term over the control volume[ x i − , x i + ] × [ t n , t n +1 ] can be divided into left and right parts, S ni = S ni − ,i + S ni,i + := 1∆ t Z t n +1 t n Z x i x i − h ( u )( α s ) x d x d t + 1∆ t Z t n +1 t n Z x i + 12 x i h ( u )( α s ) x d x d t. The CFL restriction of the time step ∆ t < ∆ x/ (2 u s ) ensures that the solid contact does not enter intothe interval [ x i , x i + ] during the time interval [ t n , t n +1 ]. Then ( α s ) x ≡ S ni,i + = 1∆ t Z t n +1 t n Z x i + 12 x i h ( u )( α s ) x d x d t = 0 . According to the essential condition (5), the integral of the non-conservative product p g ( α s ) x in[ x i − , x i ] satisfies 1∆ t Z t n +1 t n Z x i x i − p g ( α s ) x d x d t = 1∆ t Z t n +1 t n Z ξ R ξ L p g ( α s ) ξ d ξ d t = 1∆ t Z t n +1 t n Z ξ R ξ L ( α s p s ) ξ d ξ d t = ( α s ) i ( p s ) − ( α s ) i − ( p s ) . (8)Because ( u s ) = ( u s ) , the integrals of other non-conservative products in the nozzling terms in[ x i − , x i ] are obtained as1∆ t Z t n +1 t n Z x i x i − p g u s ( α s ) x d x d t = ( u s ) [( α s ) i ( p s ) − ( α s ) i − ( p s ) ] , t Z t n +1 t n Z x i x i − − u s ( α s ) x d x d t = − ( u s ) [( α s ) i − ( α s ) i − ] . From the formula (8), it is observed that the non-conservative product p g ( α s ) x is approximated by( p s α s ) x across the solid contact. This numerical approach is not consistent with the governing equations(1). In fact, neither of these approaches is robust enough. Spurious oscillations may arise in the vicinityof porosity jumps (cf. numerical results of Test 3 in [17]), especially when the solid contact approachesanother wave. 10 .1. Well-balancing in capturing almost stationary solid contacts
A well-balanced scheme refers to a scheme that accurately preserves specific steady-state solutions.For (1) a steady-state solution satisfies, f ( u ) x = h ( u )( α s ) x . (9)An individual stationary solid contact, with the solid velocity u s ≡ Figure 4: The Riemann solution in a small sector covering a stationary solid contact.
As shown in Figure 4, we consider a stationary solid contact. In the sectorial region S ( ξ L , ξ R )bounded by these two rays close to the stationary solid contact, we integrate equation (9) to get Z x R x L h ( u )( α s ) x d x = f ( u ) − f ( u ) , (10)where u and u are the states on the left and right sides of the solid contact, respectively. Ingeneral, f ( u ) = f ( u ) across the porosity jump coinciding with the cell interface. The inequalityof the numerical flux is a difficulty that needs to be settled down in the numerical simulation of thestationary solid contact. By integrating system (1) over the interval ( x i − , x i + ) from time t n to t n +1 ,the Godunov scheme for the stationary solid contact becomes u n +1 i − u ni + ∆ t ∆ x (cid:16) f n, − i + − f n, + i − (cid:17) = ∆ t ∆ x S ni = 0 , (11)where the numerical fluxes takes limiting values along the inner sides of the interface of cell I i , f n, + i − = f (cid:16) u ( x i − + 0 , t n ) (cid:17) , f n, − i + = f (cid:16) u ( x i + − , t n ) (cid:17) , S ni = 0 due to the fact that ( α s ) x = 0 in ( x i − , x i + ). This method is essentially equivalent tothe unsplit Roe-type wave propagation method in [6], which performs a wave decomposition based onthe flux difference plus the non-conservative product, inspired by (10). That is, f n, + i − = f n, − i − + Z x i x i − h ( u )( α s ) x d x := f n, − i − + S ni − , f n, − i + = f n, + i + − S ni + . Hence, (11) can be symbolically written in the form the Godunov scheme (6), u n +1 i − u ni + ∆ t ∆ x (cid:16) f n, ± i + − f n, ± i − (cid:17) = ∆ t ∆ x S ni ± . This scheme is not robust when the solid velocity u s is extremely close but not equal to 0, since itis hard to take out which numerical fluxes f n, ± i + and integrals of the nozzling term S ni ± . Later on,we design a staggered method to enhance robustness. The principle of this method is to set porosityjumps inside cells rather than at cell interfaces. Then, numerical fluxes must be invariant across cellinterfaces, and the integral of the nozzling term can be evaluated precisely. Alternatively, we may usethe residual framework to resolve this issue, as done in [30], which is left for the future work. Remark.
The 1-D homogeneous BN two-phase model is readily compared with the quasi-1-D compress-ible duct flow model [15, 31]. If the solid phase is assumed incompressible ( ρ s ) t = 0 and stationary u s ≡
0, which means that the granular bed does not move, the porosity becomes a function only asso-ciated with x . The solid phase fits the property of a duct of variable cross-section, and the gas phaseformally satisfies a duct flow model. Here, the porosity α g ( x ) acts as the variable cross-section A ( x ) inthe duct flow. Using a one-to-one correspondence ( α g , ρ g , u g , p g , E g ) ←→ ( A, ρ, u, p, E ) by removingthe subscripts g for specification compliance, system (1) reduces to the system of compressible Eulerequations in a duct u t + f ( u ) x = h ( u ) A x , (12)where notations u , f ( u ) and h ( u ) represent u = AAρAρuAρE , f ( u ) = AρuA ( ρu + p ) Au ( ρE + p ) , h ( u ) = p . Eigenvalues of this reduced system are λ = 0 , λ = u − c, λ = u, λ = u + c. In the solution to the Riemann problem for the duct flow model, the 0-characteristic field correspondsto the stationary solid contact, called 0-contact. The Riemann invariants across the 0-contact are
Aρu, η, h + u , Some numerical experiments in [15, 6, 9] revealed the fact that spurious oscillations are likelyto arise near porosity jumps by using conservative shock-capturing schemes and cannot eliminatedcompletely. In this subsection, we discuss the mechanism of these spurious oscillations. For a singlesolid contact, the Riemann invariants ψ , ψ = u s , η g , Q, P, H are constant across porosity jumps. TheGodunov average of conservative variables may yield wrongly the computation of all of these Riemanninvariants. Therefore, visible errors inevitably generate from porosity interfaces. This observation isverified in the following. solid contact Figure 5: Propagation of a solid contact in a time step.
As shown in Figure 5, we consider a single solid contact within the cell I i propagating to theright from time t n to t n +1 . Assume the solid velocity remains constant ( u s ) > α s ρ s and α g ρ g , total momentum M := α s ρ s u s + α g ρ g u g andtotal energy E := α s ρ s E s + α g ρ g E g . Specifically, we have( α s ρ s ) n +1 i = ( α s ρ s ) ni − ∆ t ∆ x δ x ( α s ρ s u s ) ni , ( α g ρ g ) n +1 i = ( α g ρ g ) ni − ∆ t ∆ x δ x ( α g ρ g u g ) ni , M n +1 i = M ni − ∆ t ∆ x δ x ( P + Qu s + M u s ) ni , where the notation of the central difference δ x ( • ) ni = ( • ) ni + − ( • ) ni − is used. Since the Riemanninvariants u s , P and Q are constant across the solid contact, δ x ( P + Qu s ) ni = 0 holds. Then Q n +1 i =[ α g ρ g ( u g − u s )] n +1 i = M n +1 i − (cid:2) ( α s ρ s ) n +1 i + ( α g ρ g ) n +1 i (cid:3) ( u s ) = M ni − [( α s ρ s ) ni + ( α g ρ g ) ni ] ( u s ) = [ α g ρ g ( u g − u s )] ni = Q ni , (13)13herefore the Godunov scheme maintains the Riemann invariant Q constant across the solid contact.In addition to u s and Q , other three Riemann invariants should remain constant across the solidcontact too, ( η g ) n +1 i = ( η g ) ni , P n +1 i = P ni , H n +1 i = H ni . (14)By the conservation property of the Godunov scheme, we determine four conservative variables at time t n +1 , ( α s ρ s ) n +1 i , ( α g ρ g ) n +1 i , M n +1 i , E n +1 i . (15)These seven quantities in (14) and (15) are independent, which will be illustrated later. However,except for the aforementioned solid velocity u s , there are only six undetermined primitive variables in(1) with the closure of the EOS for each phase and the saturation constraint,( α s ) n +1 i , ( ρ s ) n +1 i , ( p s ) n +1 i , ( ρ g ) n +1 i ( u g ) n +1 i , ( p g ) n +1 i . (16)The number of undetermined variables, six, is less than the number of independent determined quan-tities, seven, which results in an overdetermined problem. This dilemma shows that it is impossible toeliminate numerical errors by an conservative scheme in the vicinity of porosity jumps completely.Take Example 2 in Section 5 as an example. Assume the initial data to generate the solution ofa single solid contact propagating to the right and initially resting on the left boundary of the cell I i , i = 151. We set the spatial step ∆ x = 1 /
300 and time step ∆ t = 1 / t = ∆ t , we canwork out primitive variables in (16) based on all constant Riemann invariants in (14) and conservativevariables ( α s ρ s ) n +1 i , ( α g ρ g ) n +1 i , M n +1 i , and the total energy E n +1 i = ( α s ρ s E s + α g ρ g E g ) n +1 i = 10 . E n +1 i = E ni − ∆ t ∆ x δ x [ α s u s ( ρ s E s + p s ) + α g u g ( ρ g E g + p g )] ni , which updates the total energy to be E n +1 i = 10 . λ -field remain unchangedacross the solid contact. Spurious oscillations near porosity jumps cannot be fully suppressed byadjusting the numerical integration approaches for the nozzling term [17]. The only possible way toeliminate oscillations is to release the conservation, most reasonably the conservation of total energy.In Section 4, a projection method based on Riemann invariants is developed to implement the releaseof conservation. It is observed that spurious oscillations become more violent as the solid contact approaches otherstrong discontinuities, which corresponds to a resonance phenomenon. Example 4 in Section 5 shows14n instance of contiguous gaseous shock and solid contact with porosity jump. The numerical resultsin [15] displayed oscillations in the vicinity of the porosity interface. The proximity of the solid contactand gaseous shock is likely to be the cause.
Figure 6: The Riemann solution with the coalescence of gaseous shock and solid contact.
Let us consider a Riemann problem related to the coalescent gaseous shock and solid contact. Asshown in Figure 6, the Riemann solution includes a transonic solid contact coinciding with a gaseousshock, i.e., a resonant wave associated with λ ,g = λ . Inspired by the research on duct flows in [32],the resonant wave consists of three connected waves: a subsonic solid contact C s ( u , u − ) with the leftstate u and right state u − , a gaseous shock S g ( u − , u + ) and a supersonic solid contact C s ( u + , u ).These three waves coalesce on the ray ξ = ( u s ) , viewed as a region S ( ξ = u s ), where the intermediatesolid volume fraction ( α s ) − = ( α s ) + represents the porous location of the gaseous shock. Recall inSection 2.2 that the expression of the integral condition (5) is the same. However, we have to graspthe resonant wave with great accuracy; otherwise probably visible errors arise in the intermediate solidpressure ( p s ) or ( p s ) . As a matter of fact, this difficulty is hard to be avoided because of round-offerrors. It is almost impossible to judge whether the wave speeds are idential. In other words, the solidcontact is transonic.As indicated in [33], the resonance phenomenon induces instability of the flow field. Such a Riemannproblems are extremely difficult to solve. A straightforward idea is to divide these resonant wavesthrough grid cells and resolve them separately. In the next section, we implement this idea via thestrategy of staggered gas-solid grids, as illustrated in Figure 8.15 . Staggered-projection Godunov-type schemes In 2007, Saurel et al. developed a Lagrange-projection method in [34]. This method is utilizedto suppress spurious pressure oscillations that appear at contact discontinuities for compressible flowsgoverned by complex gas EOS. Inspired by the similarity between this problem and the one we studied,we explore the application of the Lagrange-projection method in order to suppress spurious oscilla-tions that appear at solid contacts. The first step of the Lagrange-projection method rests upon theLagrangian formulation. The second step is the projection of the solution on a fixed (Eulerian) grid.The homogeneous BN model contains two phases, which undoubtedly make it difficult to solve inLagrangian coordinates. To make the numerical method as simple as possible, at the first step, wecompute the solid phase in Lagrangian coordinates and the gas phase in Eulerian coordinates. Then,the question is how to design a gas-solid grid as the basis for the projection on the grid.For the coalescence of solid contacts and other gaseous waves discussed in the previous section,the difficulty of the resonance can be circumvented by separating solid contacts from other gaseouswaves. By adjusting the Godunov scheme, we assume solid contacts with porosity jumps only appearat cell centers of fixed grid cells, while other waves generate from the boundaries of these cells. Suchfixed cells are identified as gaseous grid cells. The cell centers of these gaseous cells can be regardedas boundaries of a column of solid grid cells. Then this hypothesis happens to construct a staggeredgas-solid grid. The Riemann problems at gaseous cell interfaces without porosity jumps become mucheasier to solve. Meanwhile, discontinuous interfacial fluxes commented in Subsection 3.1 no longerarise and the integrals of the nozzling term over gaseous cells can be evaluated accurately based on theformula (5). Therefore such a modified Godunov method is well-balanced. We improve the numericalintegration approach to be consistent with the governing equations (1).In [9], Karni and Hern´andez-Due˜nas had exemplified that the non-conservative scheme can suppressspurious oscillations near porosity interfaces. However, the conservation property is necessary forschemes to capture shocks. It is difficult to distinguish shocks and porosity jumps, mainly for the twowaves coinciding or interacting with each other. Therefore, it is not easy to choose the conservativeor non-conservative scheme in some parts of the computational domain. The Karni method is todetermine locations of porosity jumps by using the gradient of the porosity. The non-conservativemethod is applied near porosity jumps. In this section, we design a projection method automaticallyswitching between conservative and non-conservative formulae. It maintains a unified formulationover the entire computational domain, no matter whether the porosity is discontinuous or not. Thismethod projects solid contacts in Lagrangian coordinates back to Eulerian coordinates. The massof each phase, total momentum, and Riemann invariants for the λ - field keep constant during theprojection. The correct use of Riemann invariants suppress spurious oscillations. The total energy16s non-conservative only at porosity jumps, and conservation errors disappear gradually as jumps getsmaller. Thus the conservation of the scheme can be well guaranteed in cells containing weak solidcontacts but strong shocks. In this section, a staggered-projection Godunov-based method is proposed with detailed justifica-tion. It consists of four steps that will be specified below. For convenience of numerical implementation,the algorithm is diagramed at the end of this section. solid cellgaseous cell
Figure 7: Staggered gas-solid grid cells.
As shown in Figure 7, a black regular grid, a gaseous grid , is fixed. It is assumed that a porosity jumpis only present at the center x = x i of each gaseous grid cell I i = [ x i − , x i + ]. Red points of porosityjumps inside gaseous cells can be treated as interfaces of solid grid cells. For the Godunov scheme,the initial data are assumed to satisfy a piece-wise constant distribution at time t n . The solid volumefraction α s takes constants ( α s ) ni − and ( α s ) ni + in solid cells [ x i − , x i ] and [ x i , x i +1 ], respectively. Thevector u is constant u ni − ,i and u ni,i + in the smaller intervals [ x i − , x i ] and [ x i , x i + ], respectively.Then, the cell average of u over I i is denoted as u ni = (cid:16) u ni − ,i + u ni,i + (cid:17) . We suppose that thereis only one wave at the porosity jump, i.e., the solid contact emanating at x = x i . The Riemanninvariants ψ , ψ = u s , η g , Q, P, H , associated with λ -field, are constant across the solid contact. Then ψ ni − ,i = ψ ni,i + = ψ ni in the gaseous cell I i . In addition, we make another hypothesis that the soliddensity ρ s is constant in each gaseous cell. It means that ( ρ s ) ni − ,i = ( ρ s ) ni,i + = ( ρ s ) ni and initialdiscontinuities of the solid density only appear at gaseous cell interfaces. Based on this hypothesis,solid contact waves associated with the λ , λ ,s - fields are separated. The solid contact λ -wavesgenerate at centers of gaseous cells, and the solid contact λ ,s -waves generate at the boundaries ofgaseous cells. Hence, system (1) is decoupled into two Euler equations for the gas and solid phasesat gaseous cell interfaces. The solid phase does not affect the gas phase, and gaseous parameters donot change across solid contacts. Thus, the local Riemann problem RP (cid:16) ξ ; u ni,i + , u ni + ,i +1 (cid:17) at thegaseous cell interface x i + is easier to be solved than the Riemann problem in Section 2.2.17 igure 8: Coalescence of gaseous shock and solid contact in staggered cells. Remark.
For the coalescence of gaseous shock and solid contact in Figure 6, the staggered gas-solid gridcan decompose the resonance wave. Through an ideal schematic diagram in Figure 8, the resonantwave is divided into three waves: C s ( u , u − ) at the solid cell interface x = x i , S g ( u − , u + ) at thegaseous cell interface x = x i + and C s ( u + , u ) at x = x i +1 . In this way, these three waves can besolved individually. The Godunov scheme (6) for system (1) is utilized to update the cell average u ni over the cell I i .The numerical flux f ni + is obtained by solving the local Riemann problem RP (cid:16) ξ ; u ni,i + , u ni + ,i +1 (cid:17) atthe cell interface x i + . The time step ∆ t is restricted by the CFL condition to avoid wave interactions.Next the integral of the nozzling term is evaluated. Note that solid contacts at centers of gaseouscells propagate along fluid trajectories with the velocity ( u s ) ni . We set( p g ) max = max n ( p g ) ni − ,i , ( p g ) ni,i + o , ( p g ) min = min n ( p g ) ni − ,i , ( p g ) ni,i + o as the maximum and minimum values of gaseous pressure in I i , respectively. By the mean-valuetheorem, there exists a ( p g ) ni ∈ [( p g ) min , ( p g ) max ] such that there holds1∆ t Z t n +1 t n Z x i + 12 x i − p g ( α s ) x d x d t = 1∆ t Z t n +1 t n ( p g ) ni Z ξ R ξ L ( α s ) ξ d ξ d t = ( p g ) ni (cid:16) ( α s ) ni + − ( α s ) ni − (cid:17) , (17)where ξ = ( x − x i ) / ( t − t n ), ξ L = ( u s ) ni − ǫ and ξ R = ( u s ) ni + ǫ . For smooth solutions, this approximationis consistent with first-order accuracy.On the other hand, recall that local solutions across solid contacts satisfy the system (3). Then wehave 1∆ t Z t n +1 t n Z x i + 12 x i − p g ( α s ) x d x d t = 1∆ t Z t n +1 t n Z ξ R ξ L p g ( α s ) ξ d ξ d t = 1∆ t Z t n +1 t n Z ξ R ξ L ( α s p s ) ξ d ξ d t = ( α s ) ni + ( p s ) ni,i + − ( α s ) ni − ( p s ) ni − ,i . Comparison between these two evaluations yields( p g ) ni = ( α s ) ni + ( p s ) ni,i + − ( α s ) ni − ( p s ) ni − ,i ( α s ) ni + − ( α s ) ni − =: P g , α s ) ni + = ( α s ) ni − . Since ( p g ) ni ∈ [( p g ) min , ( p g ) max ] holds, a computation method for p g is proposedas ( p g ) ni = h ( p g ) ni − ,i + ( p g ) ni,i + i , if (cid:12)(cid:12)(cid:12) ( α s ) ni + − ( α s ) ni − (cid:12)(cid:12)(cid:12) < ε, ( p g ) min , if P g < ( p g ) min , ( p g ) max , if P g > ( p g ) max , P g , otherwise, (18)for small positive number ε , typically ε = 10 − . Then the full set of integrals for non-conservativeproducts are summarized as, in addition to (17),1∆ t Z t n +1 t n Z x i + 12 x i − p g u s ( α s ) x d x d t = ( p g ) ni ( u s ) ni h ( α s ) ni + − ( α s ) ni − i . In addition, 1∆ t Z t n +1 t n Z x i + 12 x i − − u s ( α s ) x d x d t = − ( u s ) ni h ( α s ) ni + − ( α s ) ni − i . In the case of large porosity jumps, the integral of the nozzling term S ni is computed almost exactly.In the absence of porosity jumps, ( α s ) ni + = ( α s ) ni − , the numerical integral S ni is identical to . Inthis case, this Godunov method reduces to the standard Godunov scheme about the Euler equationsfor the gas or solid phase. Figure 9: Moving solid contacts and solution distribution over time level [ t n , t n +1 ]. Denote a solid contact inside I i as( x s ) i ( t ) = x i + ( u s ) ni ( t − t n ) , t n ≤ t ≤ t n +1 , (19)as shown in Figure 9, in which β n +1 i − ,i = ( x s ) n +1 i − x i − x i + − x i − , β n +1 i,i + = x i + − ( x s ) n +1 i x i + − x i − , (20)19epresent the proportions of the interval [ x i − , ( x s ) n +1 i ] and [( x s ) n +1 i , x i + ] relative to the length ofthe whole cell I i , respectively. The average of u over the intervals [ x i − , ( x s ) n +1 i ] and [( x s ) n +1 i , x i + ]are denoted as u n +1 , ∗ i − ,i and u n +1 , ∗ i,i + , respectively. Then the cell average of u over the cell I i at time t n +1 is u n +1 , ∗ i = β n +1 i − ,i u n +1 , ∗ i − ,i + β n +1 i,i + u n +1 , ∗ i,i + , (21)which is the same as that obtained directly by the Godunov scheme (6), u n +1 , ∗ i = u ni − ∆ t ∆ x (cid:16) f ni + − f ni − (cid:17) + ∆ t ∆ x S ni . (22)The solid volume fraction α s in the scheme needs to be solved by a Lagrangian method. The advectionof the solid volume fraction indicates( α s ) n +1 , ∗ i − ,i = ( α s ) ni − , ( α s ) n +1 , ∗ i,i + = ( α s ) ni + . (23)In light of the continuity of Riemann invariants associated with the λ - field across the solid contact,we have ψ n +1 , ∗ i − ,i = ψ n +1 , ∗ i,i + , ψ = u s , η g , Q, P, H. (24)It is observed that the density ρ s is not among these Riemann invariants, thus the solid density actsas a free variable. Then we assume that the discontinuity of the solid density only appears at cellinterfaces of the cell I i , implying( ρ s ) n +1 , ∗ i − ,i = ( ρ s ) n +1 , ∗ i,i + = ( ρ s ) n +1 , ∗ i = ( α s ρ s ) n +1 , ∗ i / ( α s ) n +1 , ∗ i . (25)The Riemann invariant u s is calculated as( u s ) n +1 , ∗ i − ,i = ( u s ) n +1 , ∗ i,i + = ( u s ) n +1 , ∗ i = ( α s ρ s u s ) n +1 , ∗ i / ( α s ρ s ) n +1 , ∗ i . (26)In view of (13), the Riemann invariant Q remains constant across the solid contact and is obtainedthrough Q n +1 , ∗ i − ,i = Q n +1 , ∗ i,i + = Q n +1 , ∗ i = M n +1 , ∗ i − h ( α s ρ s ) n +1 , ∗ i + ( α g ρ g ) n +1 , ∗ i i ( u s ) n +1 , ∗ i . (27)Therefore, the Godunov scheme resolves variables α s , ρ s , u s , Q in intervals [ x i − , ( x s ) n +1 i ] and [( x s ) n +1 i , x i + ].Recalling α g = 1 − α s and making a variable substitution u g = Q/ ( α g ρ g ) + u s , it remains to solve theequations (21) and (24). Some of them are U n +1 , ∗ i = β n +1 i − ,i U n +1 , ∗ i − ,i + β n +1 i,i + U n +1 , ∗ i,i + , U = α g ρ g , α g ρ g E g , (28a) ψ n +1 , ∗ i − ,i = ψ n +1 , ∗ i,i + , ψ = η g , H, (28b)20ith four unknown thermodynamic variables of the gas phase ( ρ g ) n +1 , ∗ i − ,i , ( p g ) n +1 , ∗ i − ,i and ( ρ g ) n +1 , ∗ i,i + , ( p g ) n +1 , ∗ i,i + .The nonlinear four-equation algebraic system (28) may be solved, say, by the Newton-Raphson itera-tive method. Sometimes, the Newton-Raphson method does not converge or the nonlinear system hasno solution. This situation is closely related to the case that the Riemann problem has no solution.As described in [17], the conditions for which there is no solution are not special cases but do ariseoften. For the situations that the Newton-Raphson method is not applicable, the least-square methodis applied to solve the nonlinear system approximately. See Appendix B for details. Next, we solvesolid pressures ( p s ) n +1 , ∗ i − ,i and ( p s ) n +1 , ∗ i,i + by the two equations( α s ρ s E s ) n +1 , ∗ i = β n +1 i − ,i ( α s ρ s E s ) n +1 , ∗ i − ,i + β n +1 i,i + ( α s ρ s E s ) n +1 , ∗ i,i + , (29a) P n +1 , ∗ i − ,i = P n +1 , ∗ i,i + . (29b)For the special case where both phases meet the EOS for polytropic gases, these two equations are linearand easy to solve. Through the above processes, the Riemann invariants P n +1 , ∗ i − ,i , H n +1 , ∗ i − ,i , ( η g ) n +1 , ∗ i − ,i are evaluated.Across an individual solid contact, all Riemann invariants for the λ -field are computed exactly,suppressing spurious oscillations. This step not only relies on the conservative formulation but alsolegitimately determines the distribution of Riemann invariants. At this step, a basic requirement is that the solid contact does not cross gaseous cell interfaces.For this purpose, solid contacts are fixed in gaseous cell centers and the solid grid is an Euleriangrid, inspiring a projection procedure on the gas-solid Eulerian grid. In order to suppress spuriousoscillations in the vicinity of porosity interfaces, all Riemann invariants associated with the λ -fieldremain unchanged after the projection. In other words, the projection is based on a target that allRiemann invariants ψ = u s , η g , Q, P, H keep constant ψ n +1 i = ψ n +1 , ∗ i − ,i in the cell I i . According to theassumption (25) on the solid density ρ s , the projected solid density is constant ( ρ s ) n +1 i = ( ρ s ) n +1 , ∗ i .In fact, the projected solution can be regarded as the authentic solution at time t n +1 , as shown inFigure 10. We have obtained six independent variables ( ρ s ) n +1 i and ψ n +1 i . It remains to calculate theseventh independent variable α g or α s , which is constant in the solid cell [ x i , x i +1 ].There is a convenient method to evaluate α s . Recall in system (1) that the solid density satisfiesa conservation law ( ρ s ) t + ( ρ s u s ) x = 0 . (30)Application of the Godunov scheme over the solid cell [ x i , x i +1 ] yields( ρ s ) n +1 i + = ( ρ s ) ni + − ∆ t ∆ x (cid:2) ( ρ s u s ) ni +1 − ( ρ s u s ) ni (cid:3) , (31)21 igure 10: Projection of the solution and variable distribution in solid cells at time t n . where ( ρ s ) ni + = ( ρ s ) ni + ( ρ s ) ni +1 is the integral average of the solid density in [ x i , x i +1 ] at time t n .The Godunov scheme for the solid phase by recalling system (1) again takes( α s ρ s ) n +1 i + = ( α s ρ s ) ni + − ∆ t ∆ x (cid:2) ( α s ρ s u s ) ni +1 − ( α s ρ s u s ) ni (cid:3) (32)over the cell [ x i , x i +1 ], where the integral average ( α s ρ s ) ni + = ( α s ) ni + ( ρ s ) ni + . In the numerical fluxalong the solid cell interface x i , ( α s ) ni = ( α s ) ni − , if ( u s ) ni > , ( α s ) ni + , if ( u s ) ni ≤ . Thus the solid volume fraction at time t n +1 is updated as( α s ) n +1 i + = ( α s ρ s ) n +1 i + (cid:14) ( ρ s ) n +1 i + . (33)With ( α s ) n +1 i ± , ψ n +1 i and ( ρ s ) n +1 i available, the distribution u n +1 i − ,i and u n +1 i,i + at time t n +1 are con-structed. Here the computation of the vector u involves a root-finding process of an equation for ρ g as indicated in [9]. The specific manipulation is described detailedly in Appendix B.The whole procedure is summarized in Algorithm 1 below. In the computational region wherethe porosity is constant, the staggered-projection Godunov-type scheme completely degenerates to thestandard Godunov scheme. In the region containing porosity jumps, this scheme is conservative beforethe projection step, and during the projection process, six independent variables remain unchangedexcept for the porosity. Therefore, the projection rarely impairs the conservation of the Godunov-type scheme. In the later section, numerical experiments corroborate that this scheme significantlysuppresses spurious oscillations near porosity interfaces and accurately captures shock waves of eachphase.
Remark.
It’s worth noting that the computation of the volume fraction (33) causes slight conservationerrors of the individual-phase mass across porosity jumps, for which, however, there is a more sophisti-cated evaluation method of the volume fraction in conform with the conservation of mass. We projectthe conservative variables ( α g ρ g ) onto solid cells and then evaluate the integral average ( α g ρ g ) n +1 i + in22 lgorithm 1 First-order staggered-projection Godunov-type scheme
Input:
Cell averages ( ρ s ) ni and ψ ni in [ x i − , x i + ] for ψ = u s , η g , Q, P, H ; ( α s ) ni + in [ x i , x i +1 ]. Output: ( ρ s ) n +1 i and ψ n +1 i in [ x i − , x i + ]; ( α s ) n +1 i + in [ x i , x i +1 ]. Procedure: Solve the Riemann problem RP (cid:16) ξ ; u ni,i + , u ni + ,i +1 (cid:17) at x i + to obtain the Riemann solution u ni + = RP (cid:16) ξ = 0; u ni,i + , u ni + ,i +1 (cid:17) . Evaluate the nozzling term S ni = h (( p g ) ni , ( u s ) ni ) h ( α s ) ni + − ( α s ) ni − i with ( p g ) ni in (18). Utilize the Godunov scheme (22) to get the updated cell average u n +1 , ∗ i . Calculate the position of the solid contact ( x s ) n +1 i according to (19) and obtain the proportion β n +1 i − ,i = 1 − β n +1 i,i + in (20). Determine υ n +1 , ∗ i − ,i and υ n +1 , ∗ i,i + for υ = α s , ρ s , u s , Q by (23), (25), (26) and (27) and for υ = ρ g , p g , p s by solving the algebraic systems (28) and (29). Then determine ( ρ s ) n +1 i = ( ρ s ) n +1 , ∗ i − ,i and ψ n +1 i = ψ n +1 , ∗ i − ,i for ψ = u s , η g , Q, P, H . By the projection ( x s ) n +1 i onto x i , compute ( α s ) n +1 i + represented in (33) with ( ρ s ) n +1 i + in (31) and( α s ρ s ) n +1 i + in (32).[ x i , x i +1 ]. In fact, by assuming u s > u s < α g ρ g ) n +1 i + = (cid:18) − β n +1 i,i + (cid:19) ( α g ρ g ) n +1 , ∗ i − ,i + β n +1 i,i + ( α g ρ g ) n +1 , ∗ i,i + + 12 ( α g ρ g ) n +1 , ∗ i + ,i +1 . After the projection, the integral average is denoted as( α g ρ g ) n +1 i + = 12 ( α g ρ g ) n +1 i,i + + 12 ( α g ρ g ) n +1 i + ,i +1 . Together with the available Riemann invariants ψ n +1 i in [ x i , x i + ] and ψ n +1 i +1 in [ x i + , x i +1 ], ψ = Q, H, η g , the porosity in the solid cell [ x i , x i +1 ]( α g ) n +1 i + = ( α g ) n +1 i,i + = ( α g ) n +1 i + ,i +1 is obtained by using the Newton-Raphson iteration. Then, the solid density in the cell I i is computedby the conservation of α s ρ s ( ρ s ) n +1 i − ,i = ( α s ρ s ) n +1 , ∗ i / ( α s ) n +1 i , where ( α s ) n +1 i = 1 − ( α g ) n +1 i = 1 − h ( α g ) n +1 i − + ( α g ) n +1 i + i . Considering constant u s across porosityjumps, we may also notice that the total momentum M = ( α g ρ g + α s ρ s ) u s + Q is conservative since α g ρ g , α s ρ s and Q are conservative. The generalized Riemann problem (GRP) approach is adopted to extend the above projectionmethod to second order accuracy [35, 36]. The version we will use is that for general hyperbolic laws2326]. The resulting scheme, when combined with the above staggered projection, is a temporal-spatialcoupling staggered second order method [37]. In the vicinity of the porosity jump, algorithmic stepsthat are different from the first-order version are listed below.In the first-order staggered-projection scheme, the Riemann invariants ψ = u s , η g , Q, P, H and soliddensity ρ s are constant in the gaseous cell I i , while the solid volume fraction α s is constant in thesolid cell [ x i , x i +1 ]. In order to design a second-order version, we assume that φ and ρ s are piece-wiselinear in gaseous cells, and α s are piece-wise linear in solid cells. Since the Riemann invariants arecontinuous across solid contacts, rather than the primitive variables v or u , the local reconstructionreasonably applies for the Riemann invariants, ω := [ ρ s , u s , P, Q, H, η g ], e ω i ( x ) = ω ni + ( ω x ) ni ( x − x i ) , (34)where ω ni is calculated by the cell average u ni over I i , ( ω x ) ni is the slope of the vector ω . In the solidcell [ x i , x i +1 ], α s has a piece-wise linear representation( f α s ) i + ( x ) = ( α s ) ni + + (( α s ) x ) ni + ( x − x i + ) , (35)where ( α s ) ni + is the cell average and (( α s ) x ) ni + is the slope.With such “initial” data at time level t = t n , the GRP needs high order numerical fluxes. De-note w = [ α s , ω ], e w i,i + ( x ) = h ( f α s ) i + ( x ) , e ω i ( x ) i ⊤ and e w i + ,i +1 ( x ) = h ( f α s ) i + ( x ) , e ω i +1 ( x ) i ⊤ withlimiting states w ni + , − = lim x → x i + 12 , − e w i,i + ( x ) , w ni + , + = lim x → x i + 12 , + e w i + ,i +1 ( x ) . This GRP is written as
GRP (cid:16) e w i,i + ( x ) , e w i + ,i +1 ( x ) (cid:17) . In this paper, an acoustic GRP solver [36, 38]is utilized to approximately solve the GRP, which depends on the associated Riemann problem withthe initial data u ( x, t n ) = u (cid:16) w ni + , − (cid:17) , x < x i + , u (cid:16) w ni + , + (cid:17) , x > x i + . Denote again the associated Riemann solution as w ni + = w ( RP ( ξ = 0)). Here the recovery of thevector u from the Riemann invariants is described in detail in Appendix B. Using the acoustic approx-imation w ni + ≈ w ni + , − ≈ w ni + , + , the instantaneous temporal derivative is determined according to( w t ) ni + = R ( z t ) ni + = − [ R Λ + R − ( w x ) ni + R Λ − R − ( w x ) ni +1 ] , (36)where R = R ( w ni + ), z = R − ω , Λ = Λ ( w ni + ), Λ + = diag(max( λ i , Λ − = diag(min( λ i , R , Λ .24hen the GRP scheme updates the solution of the homogeneous BN model (1) in the followingformula, u n +1 , ∗ i = u ni − ∆ t ∆ x (cid:16) f n + i + − f n + i − (cid:17) + ∆ t ∆ x S n + i , (37)where the numerical flux is taken as f n + i + = f (cid:16) w n + i + (cid:17) , w n + i + = w ni + + ∆ t w t ) ni + . (38)The nozzling term is approximated as S n + i = h (cid:16) ( p g ) n + i , ( u s ) n + i (cid:17) h ( α s ) ni + − ( α s ) ni − i , (39)and ( p g ) n + i = h ( p g ) n + i, − + ( p g ) n + i, + i , if (cid:12)(cid:12)(cid:12) ( α s ) ni + − ( α s ) ni − (cid:12)(cid:12)(cid:12) < ε, ( p g ) min , if P g < ( p g ) min , ( p g ) max , if P g > ( p g ) max , P g , otherwise, (40)with ( p g ) max = max n ( p g ) n + i, − , ( p g ) n + i, + o , ( p g ) min = min n ( p g ) n + i, − , ( p g ) n + i, + o , P g = ( α s ) n + i + ( p s ) n + i + , + − ( α s ) n + i − ( p s ) n + i − , − ( α s ) ni + − ( α s ) ni − . The mid-point values ( p g ) n + i, − , ( p g ) n + i, + at the solid cell interface x i are given by w n + i, − = w ni, − + ∆ t w t ) ni, − , w n + i, + = w ni, + + ∆ t w t ) ni, + , which achieves second-order accuracy of the numerical integration. Here, the limiting states are w ni, − = lim x → x i − e w i − ,i ( x ) , w ni, + = lim x → x i + e w i,i + ( x ) , ( w t ) ni, − = − B ( w ni, − )( w x ) ni − , ( w t ) ni, + = − B ( w ni, + )( w x ) ni + . Then the states on both sides of solid contacts are evaluated by using the algebraic systems (28) and(29). So we obtain the Riemann invariants for the λ -field and the solid density. Then the projectionstep is the same as that for the first-order version.In order to compute the solid fractions, we first take the mid-point values on the solid interface, w n + i = w n + i, − , if ( u s ) ni > , w n + i, + , if ( u s ) ni ≤ . (41)The position of the solid contact at next time level t = t n +1 is given as( x s ) n +1 i = x i + ( u s ) n + i ∆ t. (42)25he solid volume fractions in intervals [ x i − , ( x s ) n +1 i ] and [( x s ) n +1 i , x i + ] are given in a Lagrangianstep ( α s ) n +1 , ∗ i − ,i = ( α s ) ni − − ( u s ) n + i − ∆ t ( x s ) n +1 i − x i − h ( α s ) ni − − ( α s ) n + i − i , ( α s ) n +1 , ∗ i,i + = ( α s ) ni + − ( u s ) n + i + ∆ tx i + − ( x s ) n +1 i h ( α s ) n + i + − ( α s ) ni + i , (43)respectively. The density of solid phase over the solid cell [ x i , x i +1 ] is approximated as( ρ s ) n +1 i + = ( ρ s ) ni + − ∆ t ∆ x h ( ρ s u s ) n + i +1 − ( ρ s u s ) n + i i , (44)and ( α s ρ s ) n +1 i + = ( α s ρ s ) ni + − ∆ t ∆ x h ( α s ρ s u s ) n + i +1 − ( α s ρ s u s ) n + i i , (45)where ( ρ s u s ) n + i and ( α s ρ s u s ) n + i are obtained by the mid-point value at x = x i Finally, the data need updating at next time level t = t n +1 by using the minmod limiter to suppressoscillations, particularly the slope ( ω x ) n +1 i [39, 26],( ω x ) n +1 i = minmod φ ω n +1 i +1 − ω n +1 i ∆ x , ω n +1 , − i + − ω n +1 , − i − ∆ x , φ ω n +1 i − ω n +1 i − ∆ x , (46)where ω n +1 , − i + is obtained by w n +1 , − i + = w ni + + ∆ t ( w t ) ni + andminmod( a, b, c ) = min( | a | , | b | , | c | ) , if a, b, c > , − min( | a | , | b | , | c | ) , if a, b, c < , φ ∈ [0 , α s in the solid cell [ x i , x i +1 ] is constructed as(( α s ) x ) n +1 i + = minmod φ ( α s ) n +1 i + − ( α s ) n +1 i + ∆ x , ( α s ) n +1 , − i +1 − ( α s ) n +1 , − i ∆ x , φ ( α s ) n +1 i + − ( α s ) n +1 i − ∆ x , (47)where ( α s ) n +1 , − i is chosen as( α s ) n +1 , − i = ( α s ) ni, − + ∆ t (( α s ) t ) ni, − , if ( u s ) ni > , ( α s ) ni, + + ∆ t (( α s ) t ) ni, + , if ( u s ) ni ≤ . The algorithm for the second-order GRP scheme is summarized in
Algorithm 2 . The staggered-projection Godunov-type scheme can be extended over structural meshes in twodimensions. In this section we adopt the dimensional splitting just to show this method works well.We write the two-dimensional homogeneous BN model as u t + f ( u ) x + g ( u ) y = h ( u ) · (( α s ) x , ( α s ) y ) ⊤ , (48)26 lgorithm 2 A second-order staggered-projection GRP scheme
Input:
Piece-wise linear data e ω i ( x ) in (34) in [ x i − , x i + ] for ω = [ ρ s , u s , Q, P, H, η g ]; ( f α s ) i + ( x ) in(35) in [ x i , x i +1 ]. Output: ω n +1 i and ( ω x ) n +1 i in [ x i − , x i + ]; ( α s ) n +1 i + and (( α s ) x ) n +1 i + in [ x i , x i +1 ]. Procedure: Solve the Riemann problem RP (cid:16) ξ ; u (cid:16) w ni + , − (cid:17) , u (cid:16) w ni + , + (cid:17)(cid:17) at x i + to obtain the Riemannsolution w ni + = w ( RP ( ξ = 0)). Determine ( w t ) ni + by (36) using the acoustic approximation. The numerical fluxes are given in(38). Evaluate the nozzling term S n + i in (39) by ( p g ) ni in (40). Utilize the GRP scheme (37) to get the updated cell average u n +1 , ∗ i . Calculate ( x s ) n +1 i is according to (42) and obtain β n +1 i − ,i = 1 − β n +1 i,i + in (20). Implement Procedure 5 in
Algorithm 1 . In this algorithm, ( α s ) n +1 , ∗ i − ,i and ( α s ) n +1 , ∗ i,i + are determinedby (43). With the projection from ( x s ) n +1 i onto x i , compute ( α s ) n +1 i + in (33), ( ρ s ) n +1 i + in (44) and ( α s ρ s ) n +1 i + in (45). Update the slopes ( ω x ) n +1 i and (( α s ) x ) n +1 i + by (46) and (47), respectively.with u = α s α s ρ s α s ρ s u s α s ρ s v s α s ρ s E s α g ρ g α g ρ g u g α g ρ g v g α g ρ g E g , f = α s ρ s u s α s ρ s u s + α s p s α s ρ s u s v s α s u s ( ρ s E s + p s ) α g ρ g u g α g ρ g u g + α g p g α g ρ g u g v g α g u g ( ρ g E g + p g ) , g = α s ρ s v s α s ρ s v s u s α s ρ s v s + α s p s α s v s ( ρ s E s + p s ) α g ρ g v g α g ρ g v g u g α g ρ g v g + α g p g α g v g ( ρ g E g + p g ) , h = − u s − u s p g p g p g u s p g v s − p g − p g − p g u s − p g v s , where ( u k , v k ) is the velocity for the phase k and E k = e k + ( u k + v k ). We split the system (48) intotwo subsystems u t + f ( u ) x = h ( u )( α s ) x , (49a) u t + g ( u ) y = h ( u )( α s ) y , (49b)with [ h , h ] = h , and denote L (∆ t ) , L (∆ t ) x and L (∆ t ) y as approximate solution operators for (48),(49a) and (49b) with a time step increment, respectively. Assume that the solution operators L (∆ t ) x L (∆ t ) y are space-time second-order accurate. Then the Strang splitting algorithm [27] provides asecond order approximation to L ∆ t , L (∆ t ) = L ( ∆ t ) x L (∆ t ) y L ( ∆ t ) x . (50) Figure 11: 2-D gas-solid staggered regular grid
The solution operators L (∆ t ) x and L (∆ t ) y can be specified as the staggered-projection GRP schemein the x - and y -directions, over a 2-D gas-solid staggered grid shown in Figure 11. The Riemanninvariants ψ and solid density ρ s are piece-wise linear in the black-line gaseous cells Ω ij , while thevolume fraction α s satisfies piece-wise linear distribution in the red-dashed solid cells Ω i + ,j + . Takethe solution operator L (∆ t ) x as an example. Like the staggered-projection Godunov-type scheme forthe 1-D case, within one time step, the Riemann invariants for the λ - field and the solid density inthe gaseous cell Ω ij are obtained. The cell average of ( ρ s ) n +1 i + ,j + and ( α s ρ s ) n +1 i + ,j + is updated bythe GRP scheme in the solid cell Ω i + ,j + . Then, the cell average of α s in Ω i + ,j + at time t n +1 isevaluated by ( α s ) n +1 i + ,j + = L (∆ t ) x ( α s ) ni + ,j + = ( α s ρ s ) n +1 i + ,j + (cid:14) ( ρ s ) n +1 i + ,j + . Other steps are identical to the 1-D staggered-projection GRP scheme. Based on these Riemanninvariants and the solid density in the gaseous cell, as well as the porosity in the solid cell, we canobtain the distribution of u on the gas-solid staggered grid u n +1 = L (∆ t ) x u n . Thus we have u n +1 = L (∆ t ) u n = L ( ∆ t ) x L (∆ t ) y L ( ∆ t ) x u n . . Numerical examples In this section, several numerical examples including porosity jumps are tested, separately for thequasi-1-D duct flow model, the one- and two-dimensional homogeneous BN models. In order to avoiddifficulties of thermodynamic modeling for the two phases, we assume that both phases meet the EOSfor polytropic gases e k = e k ( ρ k , p k ) = p k ( γ k − ρ k , where constant γ k = Γ k + 1 is the specific heat ratio for the phase k . The sound speed and entropyfor the phase k are represented as c k = r γ k p k ρ k , η k = p k ρ γ k k , respectively.For these examples, numerical simulations are based on regular grids with staggered gas-solid gridcells. The initial data are given by cell averages ( ρ s ) i and exact ψ i in gaseous cells and cell averages( α s ) i + in solid cells, and exact solutions are found by using the routine package CONSTRUCT[40]. The numerical results of staggered-projection Godunov-type schemes are displayed, for whichnumerical fluxes are evaluated by the exact Riemann solver or the acoustic GRP solver, and the timestep ∆ t is restricted by the CFL condition∆ t = CFL · ∆ x/ u {| λ i,k ( u ) | + | u s | ; l = 1 , , , k = s, g } with the number CFL = 0 . × s’ represent the numerical solution by thefirst-order scheme, the black ‘ + s’ represent the numerical solution by the second GRP scheme, andthe solid blue line represents the exact solution for each case. It turns out that each numerical so-lution presents essentially no spurious oscillation near porosity interfaces and totally accord with thecorresponding exact solution. The system of Euler equations in a duct of variable cross-section (15) embodies the situation ofthe stationary solid phase. It is regarded as a reduced system of the 1-D homogeneous BN model,which effectually reflects gas dynamics in the two-phase flow. At first, a numerical experiment withtwo cases is carried out for this simplified model. The staggered Godunov-type scheme, when appliedto this model, does not need to consider the motion of contacts and the projection, only involving anEulerian step actually. The computational domain [0 , .
06] is composed of M = 111 regular grid cells,29nd the initial position of the discontinuity is x = 0 .
02. The states on both sides of the discontinuityare denoted as u L and u R . The specific heat ratio of the gas is taken as γ = 1 . Table 1: Initial data for duct flows
Case A L ρ L u L p L A R ρ R u R p R I 1 151.13 212.31 2 . × . × II 1 169.34 0 2 . × × Example 1 (Shock-tube problems) . Two cases for the duct flow with a discontinuous cross-sectionare simulated. The initial data are given in Table 1. The first case is picked up from [9], correspondingto an isolated 0-contact. The numerical solution and the exact solution are shown in Figure 12(a),which shows that the current staggered scheme preserves the constant entropy across the 0-contact,approaching accurate simulation. In contrast, a conservative scheme based on a split-step algorithmin [9] fails, leading to spurious oscillations.The second case was ever considered in [6] with a finer grid, which is a Riemann problem containinga cross-section jump within a rarefaction wave, a shock and a contact discontinuity propagating tothe right. The coincidence of the 0-contact and the rarefaction generates a resonance phenomenon.The numerical solution by the current scheme and the exact solution are shown in Figure 12(b).Since the entropy cannot remain constant across the 0-contact, the numerical solution computed bythe conservative scheme based on operator splitting in [9] seriously deviated from the exact solution,probably due to the resonance. In contrast, the numerical solutions by the current scheme are in goodagreement with the exact solution.
In the first three examples for the 1-D homogeneous BN model below, the computational domain[0 ,
1] is divided into M = 300 regular grid cells, and the initial position of the discontinuity is locatedat x = 0 .
5. Initial data for all the examples are listed in Table 2. The specific heat ratios for the twophases are taken as γ s = γ g = 1 . Example 2 (Case I. A single solid contact) . This example is
Test 1 in [15]. The solution to thisproblem consists of a single solid contact propagating to the right with the velocity 0 .
3. This examplewas also simulated in [17] by the standard Godunov scheme with an exact Riemann solver, whereoscillations are present in the vicinity of the solid contact. The reason was explained in Subsection3.2. The numerical results by the current schemes are shown in Figure 13 in nice agreement with theexact solution at t = 0 .
1. 30a) Isolated 0-contact(b) Riemann problem with a cross-section jump
Figure 12: Numerical results of the staggered Godunov-type scheme and the exact solution at t = 6 . × − . igure 13: A single solid contact. Numerical results by the current staggered-projection scheme and the exact solutionare shown at t = 0 . able 2: Initial data for the 1-D homogeneous BN model Case Phase k α kL ρ kL u kL p kL α kR ρ kR u kR p kR I s (solid) 0.8 2 0.3 5 0.3 2 0.3 12.8567 g (gas) 0.2 1 2 1 0.7 0.1941 2.8011 0.1II s g s g s g s , g α sM ρ sM u sM p sM α gM ρ gM u gM p gM Example 3 (Case II. Coinciding shocks and rarefactions) . This example is
Test 2 taken in [15], andthe solution consists of two coinciding left-going shocks for the gas and solid phase and a right-goinggaseous shock within a right-going solid rarefaction wave. Our numerical results are presented inFigure 14. Compared with the numerical solution by a finite-volume Roe method presented in [15],our schemes perform better in capturing shocks and resolving the gaseous rarefaction waves.
Example 4 (Case III. A gas shock approaching a solid contact) . This example is
Test 4 in [15]. Inthe solution the solid phase contains a left-going rarefaction wave, a contact and a right-moving shock,the same as the gas phase. This example involves the Riemann problem demonstrated in Subsection3.3, i.e., the gaseous shock coincides with the solid contact, leading to a resonance phenomenon thatcauses difficulties to solve such a problem numerically. The Riemann invariants are no longer constantacross the resonant wave containing a solid contact, and the conservation property is a main factor toreduce noticeable conservation errors near the shock. For this problem, the numerical results by thecurrent staggered-projection scheme are shown in Figure 15. Compared with the numerical solutionin [15] with drastic spurious oscillations, errors in the current solutions are reduced significantly.
Example 5 (Case IV. A shock refraction at a porosity interface) . This example is taken from [9]. Agaseous shock propagates to the right and interacts with a stationary porosity interface. Initial data aregiven by u L , u M and u R in Case IV of Table 2, in which u M is the middle state between the gaseousshock and the right porosity interface. This example is simulated in the domain [0 , .
06] composed of M = 400 regular grid cells, and the initial position of the porosity interface is x = 0 .
03. This exampledemands the conservation of numerical schemes and compatible simulation of porosity interfaces. The33 igure 14: Coinciding shocks and rarefactions. Numerical results by the current staggered-projection scheme and theexact solution are shown at t = 0 . igure 15: A gas shock approaching a solid contact. Numerical results by the current staggered-projection scheme andthe exact solution are shown at t = 0 . igure 16: Shock refraction at a porosity interface. Numerical results by the current staggered-projection scheme andthe exact solution are shown at t = 0 .
007 after the refraction. t = 0 .
007 afterthe interaction of two waves. It is observed that the current schemes provides a better approximation.
Example 6 (Accuracy test) . This example is an initial value problem in [16]. We take the initial datawith a smooth transition of α s from 0 . . v s from 0 to 1, specifically, α s ( x,
0) = 0 . . x − , v s ( x,
0) = 0 . . x − , and ρ s ( x,
0) = ρ g ( x,
0) = 1 , p s ( x,
0) = p g ( x,
0) = 1 , v g ( x,
0) = 0The computational domain [0 ,
1] is divided into M cells for M = 100 , , , M = 12 , L errors and convergence orders of numerical results for u are displayed in Table 3. Thistable shows clearly that the staggered-projection Godunov-type scheme reaches first-order accuracy.In the absence of the limiter, the second-order accuracy of the staggered-projection GRP scheme isachieved. At the local extrema, the minmod limiter reduces the accuracy. Table 3: L errors and convergence orders of the vector u for Example 6 at t = 0 .
1. The numerical methods are thefirst-order and second-order staggered-projection schemes with M cells. Godunov solution GRP solution GRP solution (no limiter)
M L error order L error order L error order100 1 . × − . × − . × −
200 5 . × − . × − . × − . × − . × − . × − . × − . × − . × − We provide two kinds of examples to show the performance of the current scheme for 2-D problems.The first are the 2-D Riemann problem mimicking those in the context of gas dynamics [41, 42, 43];and the second is mimicking the shock-bubble interaction problem [44].
Example 7 (2-D Riemann problems) . The computational domain [ − . , . × [ − . , .
5] is composedof M × M square cells, and the initial data are composed of piece-wise constant states in the four37 able 4: Initial data for 2-D Riemann problems Case Region α s ρ s u s v s p s ρ g u g v g p g I γ s = 1 . γ g = 1 .
67 ( x > , y >
0) 0.8 2 0 0 2 1.5 0 0 2( x < , y >
0) 0.4 1 0 0 1 0.5 0 0 1( x < , y <
0) 0.8 2 0 0 2 1.5 0 0 2( x > , y <
0) 0.4 1 0 0 1 0.5 0 0 1II γ s = 1 . γ g = 1 . x > , y >
0) 0.6 0.5 0 0 0.6 3 0 0 0.3( x < , y >
0) 0.4 1.2 0 0 0.12 1.2 0 0 0.12( x < , y <
0) 0.6 0.5 0 0 0.6 3 0 0 0.3( x > , y <
0) 0.4 1.2 0 0 0.12 1.0 0 0 0.12quadrants, as shown in Table 4. Case I is taken from [45, 46] for which the 1-D initial data is pickedup in [17]. Case II is a 2-D Riemann problem in [47] containing a supersonic configuration,The numerical results of the staggered-projection GRP scheme are presented with M = 200 in Fig-ure 17, 18. The corresponding reference solutions are computed by the first-order staggered-projectionscheme in a finer regular grid with M = 1000. The staggered-projection GRP scheme can capturevarious waves well and displays high resolution for Cases I, II. Example 8 (Shock-cylinder interaction) . This example carries out the numerical simulation of aphysical experiment in [44] that is used as a benchmark test evaluating numerical methods of multiphaseflows. Some existing numerical results can be found in [47, 48]. In this example, a weak shock withthe shock Mach number M s = 1 .
22 propagates in atmospheric air and collides with a stationaryhelium cylinder. The computational domain [0 , . × [0 , .
89] consists of 560 ×
200 cells and the initialconfiguration is set in Figure 19, where the diameter of the cylinder is D = 0 .
5. The upper and lowerboundaries are solid walls, whereas the left and right boundaries are non-reflective boundaries. Twophases inside and outside the cylinder are regarded as polytropic gases, with the specific heat ratio γ k = 1 . k
1) and γ k = 1 .
67 for helium (denoted as the phase k Table 5: Initial data of the shock-cylinder interaction problem
Region α k ρ k p k u k ρ k p k u k inside cylinder 0.0001 1 1 0 0.1821 1 0outside cylinder, pre-shock 0.9999 1 1 0 0.1821 1 0outside cylinder, post-shock 0.9999 1.3764 1.5698 − . − . s ρ s ρ g Staggered-projection GRP scheme Reference solution
Figure 17: 2-D Riemann problem I: Numerical results of the staggered-projection GRP scheme ( M = 200) and thereference solution ( M = 1000) at t = 0 . s ρ s ρ g Staggered-projection GRP scheme Reference solution
Figure 18: 2-D Riemann problem II: Numerical results of the staggered-projection GRP scheme ( M = 200) and thereference solution ( M = 1000) at t = 0 . =0.5 Lx=2.5Ly=0.89 cylinder postshockpreshockporosityinterface
Figure 19: Diagram of the shock-cylinder interaction problem
The staggered-projection GRP scheme is applied to simulate this experimental example numerically.The phase inside the cylinder is regarded as the solid phase or gas phase separately to carry out thesimulation, and the differences between the two settings are compared. Numerical results of the totaldensity ρ = α k ρ k + α k ρ k in the shock-cylinder interaction problem are displayed in Figure 20. It isshown that the staggered-projection GRP scheme can capture interface clearly on a sparse grid. Theinstability arising at the interface due to the shock acceleration is known as the Richtmyer-Meshkovinstability. By comparing the numerical results of different settings for the cylindrical phase, it can befound that the interface shapes after the collision are different obviously. This numerical phenomenonillustrates that when the homogeneous BN two-phase flow model is utilized to simulate two separatephases approximatively, the gas phase and the solid phase have different features because of the nozzlingterms.
6. Conclusions
It is always a challenging problem to simulate compressible multi-material/phase flows due to theconflict of shock capturing and interface tracking. Moreover, the interaction of shocks and interfacesleads to the instability of flow fields, which raises more requirements of numerical methods, even thoughthere were already a lot of contributions in literature. To this end, we make our efforts in this aspectby taking the BN model.At first, we realize, through the detailed analysis of spurious numerical oscillations near materialinterfaces, that the Riemann invariants play an essential role in the design of numerical methods,they are taken as a kind of key ingredients in practice. In order to overcome the difficulty resultingfrom the conflict of conservative and non-conservative requirements, the staggered strategy is adopted,together with the projection based on the Riemann invariants. The accuracy is improved through theacoustic GRP solver. Several numerical experiments are carried out to demonstrate the reasonableperformance. 41 = 0 t = 0 . t = 0 . t = 0 . t = 1 . t = 3 k s , k g (cylindrical gas phase) k = g , k = s (cylindrical solid phase) Figure 20: Numerical density plot of the shock-cylinder interaction problem by the staggered-projection GRP schemefor two settings: The cylindrical phase as the gas phase ( k = g ) or the solid phase ( k = s ).
42t is worth noting that the Newton-Raphson iteration method is utilized to solve nonlinear algebraicequations that result from the transformation between Riemann invariants and primitive variablesand equal Riemann invariants across solid contacts. Possible non-uniqueness or non-existence of thesolutions of the algebraic equations raises a further investigation in the future in the sense that theinteraction of nozzling term and flux gradient (resonance) should be more seriously treated and theequations should be more reasonably fitted.
Appendix A. Representation of BN model in terms of Riemann invariants
For smooth solutions, system (1) is written in terms of these variables w t + B ( w ) w x = , w = [ α s , ρ s , u s , P, Q, H, η g ] ⊤ , (A.1)with the coefficient B , B = u s u s ρ s u s ρ s α s − Vρ s α s − r rT g g P k = s α k ρ k c k +3 α g ρ g V u s − rV (2 r + 1) V + c g α g ρ g ( r + 1) V rT g V (cid:20) − g P k = s α k ρ k + α s ρ s Γ g i α g ρ g V − r u g + rV α g ρ g ( r + 1) − α g ρ g T g ( r + 1)0 0 V + c g − Vρ s α s V α s ρ s + c g α g ρ g u g + rV − T g ( r + Γ g ) V u g , where V = u g − u s is the relative velocity, r = α g ρ g α s ρ s is the ratio of the mass fraction of the two phasesand Γ k = ρ k ∂p k ( ρ k ,e k ) ∂e k is the Gruneisen coefficient for the phase k . Let Λ be a diagonal matrix withthe eigenvalues λ i of B , Λ ( w ) = diag( λ i ) = diag( u s , u s − c s , u s , u s + c s , u g − c g , u g , u g + c g )and R be the right (column) eigenvector matrix of B , R ( w ) = c s c s − ρ s ρ s α s c s + α g ρ g Vρ s α s c s − α g ρ g Vρ s u g − c g − u s − α g ρ g T g Γ g V c g u g + c g − u s α g ρ g ρ s − α g ρ g ρ s − α g ρ g T g Γ g Vc g Vρ s − Vρ s − c g α g ρ g T g c g α g ρ g ,
43o that BR = R Λ . Linearizing the system (A.1) around a state w = w ∗ , it can be diagonalized toobtain z t + Λ ( w ∗ ) z x = , where z = R ( w ∗ ) − w . Appendix B. Fitting nonlinear algebraic systems
As the Newton-Raphson method are not applicable to solve the nonlinear system (28), the Gauss-Newton method is adopted to solve a least squares problem relevant to the system (28). From theequations (28a), we can obtain symbolic relationships ( p g ) n +1 , ∗ i − ,i = ( p g ) n +1 , ∗ i − ,i (cid:16) ( ρ g ) n +1 , ∗ i − ,i , ( ρ g ) n +1 , ∗ i,i + (cid:17) and( p g ) n +1 , ∗ i,i + = ( p g ) n +1 , ∗ i,i + (cid:16) ( ρ g ) n +1 , ∗ i − ,i , ( ρ g ) n +1 , ∗ i,i + (cid:17) . On the basis of the relationships, we come to solve thetwo equations ( η g ) n +1 , ∗ i − ,i − ( η g ) n +1 , ∗ i,i + = 0 ,H n +1 , ∗ i − ,i − H n +1 , ∗ i,i + = 0 . Then we consider a least squares problem in the form,minimize f ( x ) = 12 k g ( x ) k , g = h ( η g ) n +1 , ∗ i − ,i − ( η g ) n +1 , ∗ i,i + , H n +1 , ∗ i − ,i − H n +1 , ∗ i,i + i ⊤ , subject to x = h ( ρ g ) n +1 , ∗ i − ,i , ( ρ g ) n +1 , ∗ i,i + i ⊤ ∈ (0 , ∞ ) × (0 , ∞ ) . The Gauss-Newton method is utilized to minimize the least squares cost k g ( x ) k . To enhanceconvergence, a modified form in accordance with the Cholesky factorization scheme is chosen. As faras the approximate version of the Newton method does not work, we use the Newton-Raphson method.The detailed procedure of least squares can be found in [49, Section 1.4.4]. This least squares solutionapproximates the solution of the system (28) by averaging η g and H . This approximation is reasonablesince η g and H in essence only reflect the thermodynamic relationship of the gas phase. Hence theerror of these two quantities has relatively little influence on the numerical solution.The above approach can effectively enhance the robustness of the staggered-projection method.Another part for improving the robustness is to recover the vector u from the Riemann invariants ω ,which involves a root-finding process of ρ g in the equation, h ( ρ g ) = Q α g ρ g + γ g γ g − η g ρ γ g − g − H = 0 . (B.1)As pointed out in [9], this equation may have no root. It may occur when intermediate states aregenerated when a large initial jump resolves itself into waves. We also solve the equation (B.1)approximately based on a nonlinear programming problemminimize f ( ρ g ) = 12 h ( ρ g ) . H from h ( ρ g ) = 0 in (B.1).It is noted that all above processes satisfy constraints ρ g > p g > p s > References [1] M. R. Baer, J. W. Nunziato, A two-phase mixture theory for the deflagration-to-detonation tran-sition (DDT) in reactive granular materials, Int. J. Multiphase Flow 12 (6) (1986) 861–889.[2] A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, D. S. Stewart, Two-phase modeling ofdeflagration-to-detonation transition in granular materials: Reduced equations, Phys. Fluids13 (10) (2001) 3002–3024.[3] R. Abgrall, S. Karni, A comment on the computation of non-conservative products, J. Comput.Phys. 229 (8) (2010) 2759–2763.[4] S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of theequations of hydrodynamics, Mat. Sb. 89 (3) (1959) 271–306.[5] E. F. Toro, Riemann-problem-based techniques for computing reactive two-phased flows, in: Nu-merical Combustion, Lecture Notes in Physics, Springer Berlin Heidelberg, 1989, pp. 472–481.[6] C. A. Lowe, Two-phase shock-tube problems and numerical methods of solution, J. Comput. Phys.204 (2) (2005) 598–632.[7] L. Sainsaulieu, Finite volume approximation of two phase-fluid flows based on an approximateRoe-type Riemann solver, J. Comput. Phys. 121 (1) (1995) 1–28.[8] D. Bale, R. LeVeque, S. Mitran, J. Rossmanith, A wave propagation method for conservationlaws and balance laws with spatially varying flux functions, SIAM J. Sci. Comput. 24 (3) (2003)955–978.[9] S. Karni, G. Hern´andez-Due˜nas, A hybrid algorithm for the Baer-Nunziato model using the Rie-mann invariants, J. Sci. Comput. 45 (1-3) (2010) 382–403.[10] M. Dumbser, E. F. Toro, A simple extension of the Osher Riemann solver to non-conservativehyperbolic systems, J. Sci. Comput. 48 (1) (2011) 70–88.[11] M. J. Castro, P. G. LeFloch, M. L. Muoz-Ruiz, C. Pars, Why many theories of shock waves arenecessary: Convergence error in formally path-consistent schemes, J. Comput. Phys. 227 (17)(2008) 8107–8129. 4512] R. Abgrall, P. Bacigaluppi, S. Tokareva, A high-order nonconservative approach for hyperbolicequations in fluid dynamics, Comput. Fluids 169 (2018) 10–22.[13] E. Isaacson, B. Temple, Nonlinear resonance in systems of conservation laws, SIAM J. Appl.Math. 52 (5) (1992) 1260–1278.[14] P. Goatin, P. G. LeFloch, The Riemann problem for a class of resonant hyperbolic systems ofbalance laws, Ann. Inst. Henri Poincare-Anal. Non Lineaire 21 (6) (2004) 881–902.[15] N. Andrianov, G. Warnecke, The Riemann problem for the Baer-Nunziato two-phase flow model,J. Comput. Phys. 195 (2) (2004) 434–464.[16] D. W. Schwendeman, C. W. Wahle, A. K. Kapila, The Riemann problem and a high-resolutionGodunov method for a model of compressible two-phase flow, J. Comput. Phys. 212 (2) (2006)490–526.[17] V. Deledicque, M. V. Papalexandris, An exact Riemann solver for compressible two-phase flowmodels containing non-conservative products, J. Comput. Phys. 222 (1) (2007) 217–245.[18] S. A. Tokareva, E. F. Toro, HLLC-type Riemann solver for the Baer-Nunziato equations of com-pressible two-phase flow, J. Comput. Phys. 229 (10) (2010) 3573–3604.[19] R. Saurel, R. Abgrall, A multiphase Godunov method for compressible multifluid and multiphaseflows, J. Comput. Phys. 150 (2) (1999) 425–467.[20] N. Andrianov, R. Saurel, G. Warnecke, A simple method for compressible multiphase mixturesand interfaces, Int. J. Numer. Methods Fluids 41 (2) (2003) 109–131.[21] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: A quasiconservative approach, J. Comput. Phys. 125 (1) (1996) 150–160.[22] B. Re, R. Abgrall, Non-equilibrium model for weakly compressible multi-component flows: thehyperbolic operator, arXiv:1911.00270 [physics].[23] R. Abgrall, P. Bacigaluppi, B. Re, On the simulation of multicomponent and multiphase com-pressible flows, arXiv:2006.01630 [physics].[24] F. Coquel, J.-M. Hrard, K. Saleh, A positive and entropy-satisfying finite volume scheme for theBaerNunziato model, J. Comput. Phys. 330 (2017) 401–435.[25] M. D. Thanh, A well-balanced numerical scheme for a model of two-phase flows with treatmentof nonconservative terms, Adv. Comput. Math. 45 (5) (2019) 2701–2719.4626] M. Ben-Artzi, J. Li, Hyperbolic balance laws: Riemann invariants and the generalized Riemannproblem, Numer. Math. 106 (2007) 369–425.[27] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal.5 (3) (1968) 506–517.[28] P. Embid, M. Baer, Mathematical analysis of a two-phase continuum mixture theory, ContinuumMech. Thermodyn. 4 (4) (1992) 279–312.[29] I. Menshov, A. Serezhkin, A generalized Rusanov method for the Baer-Nunziato equations withapplication to DDT processes in condensed porous explosives, Int. J. Numer. Methods Fluids86 (5) (2018) 346–364.[30] M. Ricchiuto, An explicit residual based approach for shallow water flows, J.Comput.Phys. 280(2015) 306–344.[31] G. Warnecke, N. Andrianov, On the solution to the Riemann problem for the compressible ductflow, SIAM J. Appl. Math. 64 (3) (2004) 878–901.[32] E. Han, M. Hantke, G. Warnecke, Exact Riemann solutions to compressible Euler equations inducts with discontinuous cross-section, J. Hyperbol. Differ. Eq. 9 (03) (2012) 403–449.[33] T.-P. Liu, Nonlinear resonance for quasilinear hyperbolic equation, J. Math. Phys. 28 (11) (1987)2593–2602.[34] R. Saurel, E. Franquet, E. Daniel, O. Le Metayer, A relaxation-projection method for compressibleflows. Part I: The numerical equation of state for the Euler equations, J. Comput. Phys. 223 (2)(2007) 822–845.[35] M. Ben-Artzi, J. Falcovitz, A second-order Godunov-type scheme for compressible fluid dynamics,J. Comput. Phys. 55 (1) (1984) 1–32.[36] M. Ben-Artzi, J. Falcovitz, Generalized Riemann Problems in Computational Fluid Dynamics,Cambridge University Press, 2003.[37] J. Li, Fundamentals of Lax-Wendroff type approach to hyperbolic problems with discontinuities,Adv. Appl. Math. Mech. 11 (3) (2019) 38–49.[38] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, Berlin, Hei-delberg, 1997.[39] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel toGodunov’s method, J. Comput. Phys. 32 (1) (1979) 101–136.4740] N. Andrianov, CONSTRUCT: a collection of MATLAB routines for constructing exact solutionsto the riemann problem for several non-conservative, non-strictly hyperbolic systems of partialdifferential equations, https://github.com/nikolai-andrianov/CONSTRUCThttps://github.com/nikolai-andrianov/CONSTRUCT