Second Order Unconditional Positive Preserving Schemes for Non-equilibrium Reactive Flows with Mass and Mole Balance
SS econd O rder U nconditional P ositive P reserving S chemes for N on - equilibrium R eacting F lows with M ass and M ole B alance Jianhua Pan
William. G. Lowrie Departmentof Chemical and Biomolecular Engineering,Koffolt Labs, The Ohio State University [email protected]
Yu-yen Chen
William. G. Lowrie Departmentof Chemical and Biomolecular Engineering,Koffolt Labs, The Ohio State University [email protected]
Liang-Shih Fan
William. G. Lowrie Departmentof Chemical and Biomolecular Engineering,Koffolt Labs, The Ohio State University [email protected]
June 16, 2020
Abstract
In this study, a family of second order process based modified Patankar Runge-Kuttaschemes is proposed with both the mass and mole maintained in balance while preservingthe positivity of density and pressure with the time step determined by convection terms.The accuracy analysis is conducted to derive the sufficient and necessary conditions forthe Runge-Kutta and Patankar coefficients. Coupling with the finite volume method, theproposed schemes are extended to Euler equations with non-equilibrium reacting sourceterms. Benchmark tests are given to prove the prior order of accuracy and validate thepositive-preserving property for both density and pressure. K eywords non-equilibrium flow · positive preserving · modified Patankar scheme · mass conservation · mole balance · chemical reaction. Reactive non-equilibrium flows are ubiquitous in a broad range of disciplines of science and engineering.Typical examples of such phenomena are given by the reaction engineering problems, where the fluid flowtransports with the occurring of chemical reactions , such as the multiphase fluidized bed reactors [1, 2]and combustion engines [3]. Besides reaction engineering, further examples can be found in geochemistry[4], meteorology [5], etc. The one-dimensional model for such non-equilibrium flows can be expressed by asystem of reactive Euler equations, which consists a set of hyperbolic equations with source terms, U t + F ( U ) x = S ( U ) , (1) a r X i v : . [ phy s i c s . c o m p - ph ] J un preprint - J une
16, 2020where U = ( ρ , m , E , ρ Y , ρ Y , · · · , ρ Y N ) T , F ( U ) = (cid:16) m , ρ u + p , ( E + p ) u , ρ uY , ρ uY , · · · ρ uY N (cid:17) T , S ( U ) = (
0, 0, 0, ω , ω , · · · , ω N ) T , (2)and N ∑ i = Y i = m = ρ u , E = ρ u + (cid:104) ρ Y i e i ( T ) + ρ Y i h i (cid:105) , (3)where ρ is the density for the mixture of the gas; Y i is the mass fraction for i -th specie; p is the pressure, and T is the temperature. Additionally, h i is the energy corresponding to the chemical bonds and is constant.The internal energy e i and the pressure p are correlated by the state of equation, ρ Y i e i = ρ e = p γ − e and γ are the specific internal energy and specific heat ratio denoted for the mixture, respectively.The source terms satisfy the property of conservation, i.e., ∑ Ni = ω i =
0. Note that due to the conservationproperty in Eq. (1), one relation among the N + N species isnot independent, which is presented here just for completeness. Also note that the Einstein summation ruleis applied to ρ Y i e i in Eq. (3). In the remainder of this paper, unless otherwise noted, the Einstein notation isapplied to all subscripts which appear more than once in the same term and such an index is noted as thedummy index. When a dummy index is enclosed by a bracket, it is not applied by the Einstein summationrule, e.g., ρ Y ( i ) e ( i ) which simply refers to the internal energy per volume for specie i .One of the main difficulties in solving systems of Eq. (1) is how to preserve the density, pressure andtemperature in a physically bounded space. This difficult is especially significant with the presence ofsource terms, which greatly increase the stiffness of Eq. (1). For reactions that are constituted by a complexreaction network wherein multiple species and reaction steps are involved, the disparate timescale of thereactions could greatly increase the stiffness of the system.A general framework which is with high order of accuracy and is able to keep the positivity of hyperbolicsystems has been proposed by the series work of Shu, Zhang, Xing and their collaborators [6, 7, 8, 9, 10, 11].This technique has been successfully applied to shallow water equations [8], convection-diffusion equations[11], Navier-Stokes equations [9], etc., in the context of Discontinuous Galerkin (DG) method, finite volumemethod, and finite difference method. The basic idea of this framework is utilizing the convex combinationin both spatial and temporal domains. Firstly, one need to prove the positive-preserving property of abasic forward Euler scheme combing with a well-designed slope limiter which is applied to keep both theaccuracy and boundedness of the space distribution of variables. Secondly, the strong-stability-preserving(SSP) Runge-Kutta (RK) scheme which can be expressed by a convex combination of the forward Eulerscheme is utilized to realize the high-order of accuracy in time. In [6], Zhang and Shu extended the idea tohyperbolic systems with reactions under the framework of DG method, which is however suffered fromthe limited timestep when the source terms become extremely stiff.Several schemes have already been proposed to relax the constraints on the size of time step and preserve thepositivity of solutions when solving stiff equations, including the implicit-explicit (IMEX) [12, 13, 14, 15, 16]RK method, the exponential integration method [17, 18, 19], and the Patankar method [20, 21, 22, 23, 24, 25].In IMEX-RK method, the convection part and the stiff source term are discretized in an explicit and implicitmanner, respectively. The implicit treatment of stiff source terms eliminates the constraints on time stepand preserves the positivity of desired variables. Chertock, Cui, Kurganov and Wu [12] developed afamily of second order sign-preserving IMEX-RK methods for systems of ordinary differential equations(ODEs) with a stiff damping term. Following this approach, the IMEX-RK method was applied to shallow2 preprint - J une
16, 2020water equations [13], the Kerr-Debye relaxation system [14], the stiff BGK equations [15], the Cahn-Hilliardequation [16], etc. However, the main drawback of the IMEX-RK method is that it is difficult to be extendedto systems with multiple destructions and productions while maintaining the positivity of the desiredvariables in a conserved and consistent way. The idea of exponential integration method can be datedback to 1960s [17]. In this method, the stiff source term is integrated accurately with the aid of exponentialfunction, which also retains the positivity of desired variables. Recent work regarding the hyperbolicsystems comprises those of Huang and Shu [18] and Du and Yang [19]. In Huang and Shu [18], a familyof bound-preserving modified exponential RK-DG schemes was developed to solve scalar hyperbolicequations with stiff source terms. In Du and Yang [19], a third-order conservative positive-preserving andsteady-state-preserving exponential integration scheme was developed which was applied to multi-speciesand multi-reaction chemical reactive flows.Similar to the IMEX method, the Patankar method introduces an implicit modification to the stiff sourceterms to achieve the unconditional positivity [20]; note that the "unconditional" implying that the stiffsource terms do not impose constraints to the time step. Several modified Patankar methods have beenfurther proposed to ensure the conservation [21, 22, 23], to improve the order of accuracy [26, 27], and tointegrate with SSP-RK schemes [25, 24] which are popular in the simulation of hyperbolic systems. In Eulerequations with multiple chemical reactions, both the mass conservation of the species and the mole balanceof each involved element must be preserved to obtain physically reasonable solutions. The "mole balance ofeach involved element" denotes that the ratio of the destruction or production rate (with unit mol/ ( m · s ) )of different species in a reaction must equal to the ratio of the corresponding stoichiometry coefficients.However, for some versions of the Patankar-type schemes [25, 24], only the mass conservation is keptwhile the mole balance of the involved elements is violated, since the systems are expressed in a pair-wiseproduction-destruction (PD) way. For some versions of Patankar schemes, e.g., [22, 23], to keep both themass and mole balance, a common Patankar modification is applied to all the reactions. However, withsuch a common Patankar coefficient, a very stiff source term would cause all other reactions being reducedto almost zero even they are totally irrelevant. Additionally, in the Euler system with non-equilibriumreactions, the positivity of pressure cannot be guaranteed [24, 25] when the reaction is endothermic for allthe available Patankar-type schemes.In this study, we construct a second order Patankar RK scheme which is able to keep both the mass andmole balance, while preserving the positivity of pressure and is applicable to the hyperbolic systems. Toachieve the mass and mole balance simultaneously, a process based modified Patankar RK scheme (PMPRK)is proposed herein. Different from the schemes proposed by Huang, Zhao and Shu [25, 24] for PD systems,the Patankar modifications in this study are applied to reaction by reaction. Thus, the mass conservationand mole balance are maintained at the same time. An accuracy analysis is conducted to determine theconstraints for the coefficients of the SSP-RK schemes and the detailed form of the Patankar coefficients.The positivity of the proposed scheme is validated through a large number of numerical samples. Severaltheoretical proofs are given under simplified conditions. To keep the positivity of the pressure in the Eulersystem, a simple but an efficient technique is introduced. At the end, several benchmark tests are given tovalidate the proposed Patankar scheme in Sec. 6 and conclusions are provided in Sec. 7. In order to illustrate the idea of PMPRK scheme, consider the following ordinary differential equationswith only source terms, d c i d t = ω ( i ) M ( i ) = K ∑ k = ω k ( i ) M ( i ) = R k λ ki , λ ki M i = preprint - J une
16, 2020where c i is the mole concentration for the i -th species and c i = ρ Y ( i ) M ( i ) , ω ki is the mass variance (productionor destruction) for the i -the specie due to reaction k , R k is the reaction rate for the k -th chemical reaction, λ ki is the stoichiometry coefficient for specie i in reaction k , and M i is the mass per mole for specie i , with i ∈ { · · · , N } and k ∈ {
1, 2, · · · , K } . The mass conservation and mole balance can be expressed as N ∑ i = ω ki = λ ( k ) i M i R ( k ) = ω ( k ) i ω ( k ) j = λ ( k )( i ) M ( i ) λ ( k )( j ) M ( j ) , ∀ i , j ∈ {
1, 2, · · · , N } , (6) ∀ k ∈ {
1, 2, · · · , K } .The explicit second order SSP-RK scheme in the Shu-Osher form for Eq. (4) is c ( ) i = c ni , (7) c ( ) i = α c ( ) i + ∆ t β R ( ) k λ ki , (8) c ( ) i = α c ( ) i + α c ( ) i + ∆ t (cid:16) β R ( ) k + β R ( ) k (cid:17) λ ki , (9) c n + i = c ( ) i , (10)where α = α = α and β should be positive such that c ( k ) i , k =
1, 2 can bewritten as a convex combination of forward Euler schemes.In order to keep the positivity for concentrations c ( k ) i , k =
0, 1, 2, the Patankar modification is applied to thesource terms. Moreover, to keep the mole balance, the Patankar modification to the reaction rate should beprocess based, i.e., each R k adopts a same Patankar weight in all specie equations as follows, c ( ) i = c ni , (11) c ( ) i = c ( ) i + ∆ t β R ( ) k λ ki χ ( ) k , (12) c ( ) i = α c ( ) i + α c ( ) i + ∆ t (cid:16) β R ( ) k + β R ( ) k (cid:17) λ ki χ ( ) k , (13) c n + i = c ( ) i , (14)where χ ( ) k and χ ( ) k are functions defined by χ ( ) k = χ ( ) k (cid:16) c ( ) i , c ( ) i | i =
1, 2, · · · N (cid:17) , χ ( ) k = χ ( ) k (cid:16) c ( ) i , c ( ) i , c ( ) i | i =
1, 2, · · · N (cid:17) .The introduction of χ ( j ) k , j =
1, 2 makes Eqs. (12-13) implicit, and in turn makes it possible for the positivepreserving. Additionally, since the Patankar coefficients are process based, for the k -th reaction, thenumerical source terms satisfy N ∑ i = ω ki = λ ( k ) i M i R ( k ) χ ( s )( k ) = (cid:16) λ ( k ) i M i (cid:17) R ( k ) χ ( s )( k ) = s =
1, 2 , (15) ω ( k ) i ω ( k ) j = λ ( k )( i ) M ( i ) R ( k ) χ ( s )( k ) λ ( k )( j ) M ( j ) R ( k ) χ ( s )( k ) = λ ( k )( i ) M ( i ) λ ( k )( j ) M ( j ) , s =
1, 2 , ∀ i , j ∈ {
1, 2, · · · , N } , (16)which is in accordance with Eqs. (5-6), i.e., the mass conservation and mole balance are kept.4 preprint - J une
16, 2020Let us first discuss the sufficient and necessary conditions for the PMPRK scheme to achieve the prior orderof accuracy, and the detailed expressions for the Patankar coefficients χ ( j ) k , j =
1, 2 will be provided later.Since the proposed PMPRK is just a minor adjustment of the original scheme, the truncation error of themodification coefficient χ ( ) k and χ ( ) k is assumed to have the following form, χ ( ) k = + X k ∆ t + O ( ∆ t ) , (17) χ ( ) k = + Y k ∆ t + O ( ∆ t ) . (18)To fulfill the requirements for the order of accuracy, there should exist certain constraints on X k and Y k .Substituting Eq. (17) into Eq. (12), we obtain c ( ) i = c ( ) i + ∆ t β R ( ) k λ ki ( + X k ∆ t ) + O ( ∆ t ) . (19)Meanwhile, expand the expression of c ( ) i near time t ( ) and substitute Eqs. (18) into it, c ( ) i = α c ( ) i + α c ( ) i + ∆ t (cid:110) β R ( ) k + β (cid:104) R ( ) k + R ( ) k , j ( c ( ) j − c ( ) j ) (cid:105)(cid:111) λ ki · ( + Y k ∆ t ) + O ( ∆ t ) , (20)where the subscript after the comma in R k , j denotes the index of the specie to which the operation of partialderivative is applied to, i.e., R k , j = ∂ R k ∂ c j . (21)Substitute Eq. (19) into Eq. (20), the final expression of c ( ) i can be obtained.On the other hand, the Taylor expansion of exact solution (cid:101) c i near c ( ) i is (cid:101) c i = c ( ) i + R k λ ki ∆ t + R k , j R m λ mj λ ki ∆ t + O ( ∆ t ) . (22)Subtracting (cid:101) c i from c ( ) i , the following equation is obtained, c ( ) i − (cid:101) c i = θ c ( ) i + θ ∆ t + θ ∆ t + O ( ∆ t ) , (23)where θ = α + α − θ = [( α β + β + β ) − ] R k λ ki , θ = [ α β X k + ( β + β ) Y k ] R k λ ki + (cid:18) β β − (cid:19) R m λ m , j R k , j λ k , i . (24)In order to achieve a prior second order of accuracy, θ k , k =
1, 2, 3 should all be 0. Thus, the sufficient andnecessary conditions for the second order PMPRK scheme is α = α + α = α β + β + β = β β =
12 , (25) α β X k + ( β + β ) Y k = preprint - J une
16, 2020 χ ( ) k and χ ( ) k The key novelty of this work is that the Patankar modification in the proposed PMPRK scheme is processbased and applied to the source terms reaction by reaction, and thus should be the same for all the speciesinvolved in the same reaction.An intuitive construction of χ ( ) k and χ ( ) k is χ ( ) k = ∏ ε ∈ Ret ( k ) c ( ) ε π ε q k , (27) χ ( ) k = ∏ ε ∈ Ret ( k ) c ( ) ε τ ε q k , (28)where Ret ( k ) denotes the set of reactants (the species which are destructed) in the k -th reaction; q k is areal number locates in ( + ∞ ) . In this study, it is assumed that a specie cannot be both a reactant anda production. If a specie does show up in both side of the reaction, the side with smaller stoichiometrycoefficient is to be cancelled. π ε is a function of c ( ) ε ; τ ε is a function of c ( ) ε , π ε and c ( ) ε .The form of χ ( j ) k , j =
1, 2 only involves the reactants of reaction k ; therefore, if the destructed reactantapproaches 0, χ ( j ) k , j =
1, 2 approaches zero, thereby playing as a scaling factor for reaction k . As will beshown later, numerical investigation of both randomly generated samples and the benchmark tests bothdemonstrate the effectiveness of the proposed form of χ ( j ) k , j =
1, 2.Before discussing the positivity of the scheme, let us first determine the expression of π ε and τ ε whichcan preserve the prior order of accuracy. The most straightforward design of π ε is like that of the work inHuang and Shu [24] and Huang et al.[25], π ε = c ( ) ε . (29)Thus, we have χ ( ) k = + O ( ∆ t ) . (30)The substitution of Eq. (30) into Eq. (12) yields c ( ) i = c ( ) i + ∆ t β R k λ ki + O ( ∆ t ) . (31)Then, substituting Eq. (31) into Eq. (27), the first order expansion of χ ( ) k becomes χ ( ) k = ∏ ε ∈ Ret ( k ) c ( ) ε c ( ) ε q k = + ∆ tq k β ∑ ε ∈ Ret ( k ) R m λ m ε c ( ) ε + O ( ∆ t ) . (32)Thus, X k is expressed as X k = q k β ∑ ε ∈ Ret ( k ) R m λ m ε c ( ) ε . (33)From Eq. (26), Y k is obtained as Y k = − α β β + β q k β ∑ ε ∈ Ret ( k ) R m λ m ε c ( ) ε . (34)Since Ret ( k ) is an arbitrary set containing the indices of the reactants, we have c ( ) i τ i = − α β β + β β R k λ ki c ( ) i ∆ t + O ( ∆ t ) . (35)6 preprint - J une
16, 2020Thus, the first order expansion of τ i c ( ) i is τ i c ( ) i = c ( ) i / c i c ( ) i / τ i = + R k λ ki c ( ) i ∆ t + O ( ∆ t ) − α β β + β β R k λ ki c ( ) i ∆ t + O ( ∆ t )= + (cid:18) + α β β + β β (cid:19) R k λ ki c ( ) i ∆ t + O ( ∆ t ) . (36)Like the work in Huang and Shu [24], we can construct the τ i in the form as the following equation, τ i c ( ) i = (cid:32) c ( ) i c ( ) i (cid:33) s . (37)Then the parameter s is determined as s = (cid:18) β + α β β + β (cid:19) . (38)In fact, τ i can be constructed in any form which can recover its first order expansion and can be expressedas a convex combination of several positive terms. Equation (38) is consistent with the form in Huangand Shu [25]. But the analysis shown in this study reveals that Eq. (38) does not only apply for the linearPatankar coefficients as in the PD systems, but also suitable for the nonlinear Patankar coefficients asproposed in Eqs. (27-28). Remark 1.
The analysis shows that the real number q k in Eqs. (27) and (28) does not affect the prior order ofaccuracy and is thus a free parameter. However, according to the Taylor expansion of χ ( ) k and χ ( ) k , a larger q k usuallycorresponds to a larger truncation error. Additionally, through numerical investigations in Sec. 4.2 with randomlygenerated samples, with the form of constant q for all χ k , e.g., q = , if the number of reactants increases, the cut-offerror due to the limited bytes of floating numbers makes the calculations of Patankar coefficients difficult and makes ithard to converge for the Newton iterations, which is utilized to solve the implicit function of Eqs. (12-13). Thus, inthe current work, we choose q k = Ret ( k ) , where Ret ( k ) denotes the number of reactants for the k-th reaction. For the PMPRK scheme, a system of implicit equations needs to be solved. Given only the destructedreactants are involved in the expression of χ ( j ) k , j =
1, 2; c ( j ) i , j =
0, 1 and ∆ tR ( k ) λ ( k ) i in Eqs. (12-13), π i and τ i , i =
1, 2, · · · , N in Eqs. (27-28) are all non-negative numbers, the solving of implicit Eqs. (12-13) isequivalent to solving Eq. (40). Conjecture 1.
Given a i , a ∗ i , b ki and q k as constants satisfying a i > a ∗ i > b ki ≥ , q k ∈ ( + ∞ ) , c i asunknowns and i =
1, 2, · · · , N ; k =
1, 2, · · · , K, additionally, ∑ i ∈ Ret ( k ) b ki ≥ ∑ i (cid:54)∈ Ret ( k ) b ki . (39) If a system of equations can be expressed in the formc i − a i + ∑ k , Ret ( k ) (cid:51) i b ki ∏ j ∈ Ret ( k ) c j a ∗ j q k − ∑ k , Ret ( k ) (cid:54)(cid:51) i b ki ∏ j ∈ Ret ( k ) c j a ∗ j q k = it has at least one solution. And among these solutions, there exists a unique solution satisfying c i > i =
1, 2, · · · , N. preprint - J une
16, 2020
Remark 2.
Conjecture 1 is an extension of the linear system in [24], and the constraint of Eq. (39) is introduced tobe consistent with the linear system in [24] where a strict diagonal dominant m-matrix can be obtained. However,through the numerical tests, the constraint of Eq. (39) can be relaxed a little for the nonlinear system when thepressure is introduced as shown in Secs. 5 and 6.6.
Remark 3.
The density of the reactants can be zero, i.e., a i and a ∗ i are zero, which leads the corresponding reactionrate to be also zero and thus be decoupled from Eq. (40). Therefore, only the system with strict positive densities for allspecies are discussed in this section. The correctness of Conjecture 1 is unknown. However, if we take a glance at Eq. (40), it can be foundthat Eq. (40) can be treated as a first order backward Euler discretization with monotone source terms;thus Conjecture 1 is assumed to be correct in this study. In this section, we give two conclusions undersimplified conditions. For more general cases, the conjecture is validated with exhausted numerical samplesin the parameter space of a i , a ∗ i and b ki . In the remainder of this paper, if all the elements of a vector or amatrix is greater than zero, or greater or equal than zero, it is denoted as (cid:126) c > M > (cid:126) c ≥ M ≥ If each of the reactions has only one reactant, then the reactions can be expressed by the indices of thespecies and denoted as q k = q i , ∀ Ret ( k ) (cid:51) i. Furthermore, Eq. (40) has at least one solution, and there exists an oneand only solution satisfying (cid:126) c > among all the solutions.Proof. After a variable substitution of v i = c q i i , Eq. (40) can be rewritten in the following form, (cid:126) v q + M (cid:126) v = (cid:126) a , (41)where (cid:126) v q = (cid:26) v q , v q , · · · , v qN N (cid:27) T and M ii = ∑ k ,Ret ( k ) (cid:51) i b ki a ∗ i q i , M ij = − ∑ k ,Ret ( k ) (cid:51) i b kj a ∗ i q i , i (cid:54) = j , (cid:126) v = { v , v , · · · , v N } T , (cid:126) a = { a , a , · · · , a N } T . (42)Note that the Einstein summation rule does not apply to the subscript of M ii in Eq. (42) and all the otherexpressions in the remainder of this section. According to Eq. (39) and considering i is the only reactantin Ret ( k ) , b ki ≥ ∑ j (cid:54) = i b kj . Thus M is a column diagonally dominant m -matrix. Then we introduce a smallnumber (cid:101) > ( (cid:101) I + D ) − (cid:18) (cid:126) v q + D (cid:126) v (cid:19) = ( (cid:101) I + D ) − ( R (cid:126) v + (cid:126) a ) , (43)where I is an identity matrix, D is the diagonal part of M , M = D − R and D ≥ R ≥ ∀ (cid:126) v >
0, wehave ( (cid:101) I + D ) − ( R (cid:126) v + (cid:126) a ) > (cid:126) F ( (cid:126) v ) = ( (cid:101) I + D ) − (cid:18) (cid:126) v q + D (cid:126) v (cid:19) , (44)and since (cid:101) I + D is a diagonal matrix, it is easy to verify that, in the first phase, ∀ (cid:126) f > (cid:126) F ( (cid:126) v ) = (cid:126) f hasone and only one solution. Thus, we denote the projection (cid:126) F − : (cid:126) f → (cid:126) v , (cid:126) v > (cid:126) f = ( (cid:101) I + D ) − ( R (cid:126) v n + (cid:126) a ) , (cid:126) v n + = (cid:126) F − ( (cid:126) f ) . (45)8 preprint - J une
16, 2020Since (cid:101) > M is a column diagonally dominant m -matrix, (cid:101) I + M is a strict column diagonallydominant m -matrix, the spectral radius of ( (cid:101) I + D ) − R is smaller than 1, and there must exist a p -normwhich makes (cid:13)(cid:13)(cid:13) ( (cid:101) I + D ) − R (cid:13)(cid:13)(cid:13) p <
1. Hence, for any given (cid:126) v n > (cid:126) v n >
0, we have (cid:126) f > (cid:126) f > (cid:13)(cid:13)(cid:13) (cid:126) f − (cid:126) f (cid:13)(cid:13)(cid:13) p ≤ (cid:13)(cid:13)(cid:13) ( (cid:101) I + D ) − R (cid:13)(cid:13)(cid:13) p (cid:107) (cid:126) v n − (cid:126) v n (cid:107) p < (cid:107) (cid:126) v n − (cid:126) v n (cid:107) p .On the other hand, if we can find a (cid:101) which makes (cid:13)(cid:13)(cid:13) (cid:126) v n + − (cid:126) v n + (cid:13)(cid:13)(cid:13) p ≤ (cid:13)(cid:13)(cid:13) (cid:126) f − (cid:126) f (cid:13)(cid:13)(cid:13) p , the iteration scheme isconvergent.The solution of the second step in Eq. (45) can be written element-wisely as v qi i = − D ii v i + ( (cid:101) + D ii ) f i , (46)where the superscript n + v i is omitted without confusion. Since R ≥ f i ≥ a i (cid:101) + D ii , the x-intercept ofthe linear function on the right-hand side of Eq. (46) satisfies ( (cid:101) + D ii ) f i D ii ≥ a i (cid:101) + D ii (cid:101) + D ii D ii = a i D ii ,and the y-intercept satisfies ( (cid:101) + D ii ) f i ≥ a i (cid:101) + D ii ( (cid:101) + D ii ) = a i .Denoting the first-phase solution of v qi i = − D ii v i + a i as (cid:101) v i which is greater than 0 and its vector form as (cid:126) (cid:101) v ,it can be concluded that the solution of Eq. (46) is a bounded value locates in [ (cid:101) v i , + ∞ ) .Meanwhile, ∂ (cid:126) F ∂ (cid:126) v = ( (cid:101) I + D ) − (cid:18) Diag ( q (cid:126) v q − ) + D (cid:19) , (47)where Diag ( · ) means a diagonal matrix and1 q (cid:126) v q − = (cid:26) q v q − , 1 q v q − · · · , 1 q N v qN − N (cid:27) T .If (cid:101) = min ( q i (cid:101) v qi − i ) , i =
1, 2, · · · , N , each component of (cid:126) f − (cid:126) f and (cid:126) v n + − (cid:126) v n + satisfies | (cid:126) f i − (cid:126) f i | p | (cid:126) v n + i − (cid:126) v n + i | p = (cid:32) ∂ (cid:126) F ∂ (cid:126) v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126) v ∗ (cid:33) pii ≥ (cid:126) v ∗ ∈ [ (cid:126) v n + , (cid:126) v n + ] ⊂ [ (cid:126) (cid:101) v , + ∞ ) .Thus, we have (cid:13)(cid:13)(cid:13) (cid:126) v n + − (cid:126) v n + (cid:13)(cid:13)(cid:13) p ≤ (cid:13)(cid:13)(cid:13) (cid:126) f − (cid:126) f (cid:13)(cid:13)(cid:13) p ≤ (cid:13)(cid:13)(cid:13) ( (cid:101) I + D ) − R (cid:13)(cid:13)(cid:13) p (cid:107) (cid:126) v n − (cid:126) v n (cid:107) p < (cid:107) (cid:126) v n − (cid:126) v n (cid:107) p .In conclusion, with the selected (cid:101) , for arbitrary given first phase initial values, the iteration scheme definedas in Eq. (45) converges to the only one solution. Given (cid:126) v n > (cid:126) v n + is also greater than 0 and theconverged solution locates in the first phase. The proof is completed.If q i = ∀ i =
1, 2 · · · , N , Theorem 1 reduces to the case of the PD system as in [25, 24]. From Theorem 1,it is concluded that the power of the Patankar coefficients can be chosen as any positive value. Theorem 2.
Among all the reactions, there is only one reaction involving with more than one reactant but withoutproduction. Without losing generality, assume
Ret ( ) > and Ret ( ) = {
1, 2, · · · , k } . q = p. For the reactionwith only one reactant, q k = k =
2, 3, · · · , K. Under such conditions, Eq. (40) has at least one solution, and thereis one and the only one satisfying (cid:126) c > . preprint - J une
16, 2020
Proof.
Equation (40) can be rewritten as M (cid:126) c = (cid:126) a + (cid:126) Bt , (48)where M has the same definition as in Eq. (42) except that the summation only applies to reactions from 2to K . (cid:126) B = {− b , − b , · · · , − b k , 0, 0, · · · , 0 } T and t = (cid:16) c c ··· c k a ∗ a ∗ ··· a ∗ k (cid:17) q .Thus (cid:126) c = M − (cid:16) (cid:126) a + (cid:126) Bt (cid:17) . (49)Since M is a column diagonally dominant m -matrix, all entries of M − is non-negative, and the diagonalelements are positive and are greater than other elements in the same row, denote the right-hand side of Eq.(49) as M − (cid:16) (cid:126) a + (cid:126) Bt (cid:17) = { ξ − ζ t , ξ − ζ t , · · · , ξ k − ζ k t , 0, 0, · · · , 0 } T ,where ξ i , ζ i , i =
1, 2, · · · , k are positive values. Then, we have t = (cid:18) ( ξ − ζ t ) ( ξ − ζ t ) · · · ( ξ k − ζ k t ) a ∗ a ∗ · · · a ∗ k (cid:19) q . (50)Since the first phase solution should satisfy (cid:126) c >
0, we have t < min ( ξ i ζ i ) , i =
1, 2, · · · , k . (51)With constraint of Eq. (51), in the first phase, Eq. (50) has one and only one solution. Thus, in the firstphase, Eq. (48) has one and only one solution. The proof is completed. For general cases, we have tried to prove the correctness of Conjecture 1 but have not yet found theapproach. Thus, this problem is remained to be solved. Instead of providing a rigorous and mathematicallyelegant proof, in this study, we try to validate this conjecture through large randomly generated samples.A system with 100 species ( N = K = a i , a ∗ i and b ki are drawn out from a uniform distribution in [
0, 1 ] and satisfyEq. (39). For a real reactive system, the value of a i , a ∗ i and b ki of course can be of any value. However, itshould be noted that, if c i , a i , a ∗ i and b ki hold for Eq. (40), after a scaling, hc i , ha i , ha ∗ i and hb ki still holdfor Eq. (40), where h is an arbitrary real number. This means, for any sets of coefficients, that they canbe scaled and locate in the region of [
0, 1 ] . Thus, it is reasonable to draw the samples of a i , a ∗ i , b ki outfrom [
0, 1 ] . For the power of the Patankar coefficient, case with q k = ( k ) is considered, which is also aconsideration from the viewpoint of truncation error and cut-off error.Since obtaining all the solutions for such a big system is very time-consuming, Newton iteration isused to solve the equations in the real space. 1000000 samples are drawn out and tested. Denotes χ k = (cid:16) ∏ ε ∈ Ret ( k ) c ε a ∗ ε (cid:17) q k , v i = c q im i , f i = c i − a i + ∑ k ,Ret ( k ) (cid:51) i [ b ki χ k ] − ∑ k ,Ret ( k ) (cid:54)(cid:51) i [ b ki χ k ] , J ij = ∂ f i ∂ v j and J = (cid:8) J ij (cid:9) ,where q im = min Ret ( k ) (cid:51) i { q k } . To keep the positivity of c i during the iteration, Newton iterative algorithmwith a re-initialization strategy is briefly listed in Algorithm 1, which is shown to be effective in all thetesting cases and numerical examples in Sec. 6.For each sample, 10 different randomly generated initial values of c i are tested and converge to the samesolution. Since the data for the tested coefficients of Eq. (40) are too large for the presentation, the sourcecode for generating and testing samples are directly provided in Pan [30] as a supporting information.Interested readers can download and test it. Although being not a rigorous proof, the result of this sectioncan provide as a solid supporting to Conjecture 1, which is the basis of the proposed scheme.10 preprint - J une
16, 2020
Algorithm 1
Newton iteration for the Patankar scheme. c i ← c i v i ← c q im i while f i > Threshold do f i ← c i − a i + ∑ k ,Ret ( k ) (cid:51) i [ b ki χ k ] − ∑ k ,Ret ( k ) (cid:54)(cid:51) i [ b ki χ k ] for k = K do χ k ← (cid:16) ∏ ε ∈ Ret ( k ) c ε a ∗ ε (cid:17) q k . end for ∆ v = − J − v . for i = N doif v i + ∆ v i < then c i ← a i + ∑ k ,Ret ( k ) (cid:54)(cid:51) i [ b ki χ k ] c i ← min (cid:16) c i , ∑ Ni = a i (cid:17) v i ← c q im i else v i ← v i + ∆ v i c i ← v q im i end ifend forend while For the Euler equations coupled with convection terms, the positive-preserving technology proposed inZhang and Shu [31] is firstly reviewed and utilized to guarantee the positivity of the convection part. Thenthe PMPRK scheme is used to keep the positivity concerning the source terms. As will be shown later, thePMPRK scheme of the form Eq. (27) can be easily coupled with the energy equation and help keep thepositivity of the pressure.If we define the admissible set G as in Huang and Shu [24] and Huang et al. [25], that is G = (cid:110) U = ( ρ , Y , Y , · · · , Y N , m , E ) T | ρ > Y i > i =
1, 2, · · · , N , p > (cid:111) , (52)it was found that if the reaction is endothermic, G may not be a convex set due to the existence of h i , seeSection 3.3 and 3.4 in Huang and Shu [24] and thus the pressure could not be guaranteed to be positive asstated. However, if we make a little rearrangement of Eq. (1), this problem is readily solved. Rewrite Eq. (2)as U = ( ρ , m , E , ρ Y , ρ Y , · · · , ρ Y N ) T , F ( U ) = (cid:16) m , ρ u + p , ( E + p ) u , ρ uY , ρ uY , · · · ρ Y N (cid:17) T , S ( U ) = (
0, 0, ω E , ω , ω , · · · , ω N ) T . (53)and ω E = − ∂ρ Y i h i ∂ t − ∂ρ Y i h i u ∂ x , E = ρ u + ρ Y i e i ( T ) . (54)11 preprint - J une
16, 2020All other definitions are the same as in Eq. (2). Further, according to the mass conservation equations andconsidering that h i is constant for each specie, ω E = − K ∑ k = R k λ ki M i h i . (55)Equations (53-55) are the final forms that we solved for the Euler equations in this study. Then, we definethe admissible set of G in Eq. (52) with E defined as in Eq. (54). After rewriting, G becomes a convex set asin Zhang and Shu [31, 32] and can be preserved to be positive if only the source term of energy equation istreated appropriately.Assume that the computational domain constitutes M uniformly distributed cells and cell j is defined in [ x j − , x j + ] . Consider the following finite volume scheme with forward Euler time integration, U n + j = U nj − ∆ t ∆ x (cid:20) ˆ F (cid:18) U − j + , U + j + (cid:19) − ˆ F (cid:18) U − j − , U + j − (cid:19)(cid:21) + ∆ tS n , (56)where ∆ t and ∆ x are the time and space step, respectively. U j means the averaged value in cell j . Thesuperscripts n and n + n and n +
1, respectively. Denotethe piece-wise high order reconstructed polynomial based on the cell averaged value U n for each cell as Q nj ( x ) . U − j + = Q nj ( x j + ) and U − j − = Q nj ( x j − ) . ˆ F is the numerical flux. S n is defined as S n = ∆ x (cid:90) x j + x j − S ( Q j ( x )) . (57)Let x η j be the W -point Legendre Gauss-Lobatto quadrature points for the interval [ x j − , x j + ] and ˆ w η bethe corresponding weights int the standard cell [ − ] such that ∑ W η = ˆ w η = W − Q nj ( x ) . Theorem 3. [31] For a finite volume scheme and a positive preserving flux ˆ F, if Q nj ( x α j ) ∈ G for all j and all α , thenU nj − ∆ t ∆ x (cid:20) ˆ F (cid:18) U − j + , U + j + (cid:19) − ˆ F (cid:18) U − j − , U + j − (cid:19)(cid:21) ∈ G under the CFL condition (cid:107)| u | + acou (cid:107) ∞ ∆ t ∆ x ≤ ˆ w δ . (58) acou is the acoustic wave speed defined as acou = (cid:112) γ p / ρ , δ is a constant number depending on thepositive preserving flux, e.g., the Godunov flux [33], the Lax-Friedrichs flux [34], the Boltzmann type flux[35] , the Harten-Lax-van Leer flux [36] and the HLLC [37, 38] flux. For the Lax-Friedrichs flux, δ canreach 1 [31] and for four-point Legendre Gauss-Lobatto quadrature points, ˆ w = .Then incorporating the PMPRK scheme, the final single-step forward Euler time integration scheme of theEuler equations can be written as ρ n + j = (cid:20) ρ nj − ∆ t ∆ x ( ˆ F j + ,1 − ˆ F j − ,1 ) (cid:21) , m n + j = (cid:20) m nj − ∆ t ∆ x ( ˆ F j + ,2 − ˆ F j − ,2 ) (cid:21) , ρ e n + j = E nj − ∆ t ∆ x ( ˆ F j + ,3 − ˆ F j − ,3 ) − m n + j ρ n + j − ∆ t K ∑ k = R nk λ ki M i h i χ k , ρ Y i , jn + = (cid:20) ρ Y i , jn − ∆ t ∆ x ( ˆ F j + ,3 + i − ˆ F j − ,3 + i ) (cid:21) + ∆ tR nk λ k ( i ) M ( i ) χ k , i =
1, 2, · · · , N , (59)12 preprint - J une
16, 2020where ρ e = E − m ρ . According to Theorem 3, the values inside the square bracket in Eq. (59) are allpositive. In Eq. (59), ρ e can be treated as a special kind of specie. If λ ki M i h i < ρ e is treatedas a production; if λ ki M i h i > ρ e is treated as a reactant. However, the introduction of ρ e equation breaks down the constraints of Eq. (39), i.e., if ρ e is treated as a specie and if the reaction isexothermic, ∑ i ∈ Ret ( k ) b ki < ∑ i (cid:54)∈ Ret ( k ) b ki . . Luckily, through the numerical investigations in a same way asdescribed in Sec. 4.2, even for the general system with 100 species and 100 reactions with each of whichbeing either exothermic or endothermic, there exists a unique solution for the PMPRK scheme, whichis different from the purely linear case in Huang and Shu [24] and Huang et al. [25], where a diagonaldominant requirements of Eq. (39) must be fulfilled.To illustrate this, let us consider a very specific case with two species and two reactions. Each reactioninvolves only one reactant, and one reaction is endothermic, the other is exothermic. The system of PMPRKscheme, when written in the form of Eq. (40), can be expressed as c = a − b c c ∗ + b (cid:114) c c c ∗ c ∗ , c = a + b c c ∗ − b (cid:114) c c c ∗ c ∗ , c = a + b c c ∗ − b (cid:114) c c c ∗ c ∗ , (60)where c represents the equation for the internale energy and b ≥ b , b ≤ b . Since b ki and c ∗ i , i =
1, 2, 3; k =
1, 2 are all constants, denoting b ∗ i = b i / c ∗ , i =
1, 2, 3 and b ∗ i = b i / (cid:112) c ∗ c ∗ , Eq. (60) can besimplified as c = a − b ∗ c + b ∗ √ c c , c = a + b ∗ c − b ∗ √ c c , c = a + b ∗ c − b ∗ √ c c , (61)and b ∗ ≥ b ∗ , b ∗ ≤ b ∗ . Denoting t = √ c c , c can be expressed as c = a + b ∗ t + b ∗ , the substitution of whichinto c and c yields c = a + a b ∗ + b ∗ + ( b ∗ b ∗ + b ∗ − b ∗ ) t , c = a + a b ∗ + b ∗ + ( b ∗ b ∗ + b ∗ − b ∗ ) t . (62)Since b ∗ ≤ b ∗ , b ∗ b ∗ + b ∗ − b ∗ < b ∗ − b ∗ ≤
0. Considering that c c = t and c and c in Eq. (62) should allbe greater than 0, it can be concluded that there exists a unique t satisfying t > c > c > c > ρ e n + , ρ Y in + given that ρ e n > ρ Y in > i =
1, 2, · · · , N , that is, the positivity of both density and pressure are preserved.The extension to high order SSP Runge-Kutta (SSP-RK) scheme is straightforward and as described in theliterature work [6, 7, 8, 9, 10, 11], i.e., for each step of SSP-RK, the update scheme can be written as a convexcombination of the forward Euler scheme.To preserve the positivity, the restriction of the CFL condition is (cid:107)| u | + acou (cid:107) ∞ ∆ t ∆ x ≤ ˆ w δ min α ij β ij . (63)To maximize min α ij β ij , for current second order PMPRK scheme, the following RK coefficients are selected, α = β = α = α = β = β = preprint - J une
16, 2020with min α ij β ij = Firstly, the accuracy of the proposed PMPRK scheme is verified through the following ordinary differentialequations, d c d t = − λ c c + λ c c c ,d c d t = − λ c c − λ c c c ,d c d t = λ c c − λ c c c ,d c d t = − λ c c c , (65)where two reactions are involved, λ and λ are constants. Initial values are { c , c , c , c } = { } .Time step ∆ t = t end / N t . The results of N t = { λ , λ } = { } and { λ , λ } = (cid:8) × , 2 × (cid:9) , the latterof which represents a system with stiff source terms. The final time is t end = { λ , λ } = { } , the second prior order of accuracy is observed. For the stiff system with { λ , λ } = (cid:8) × , 2 × (cid:9) , super-convergence is observed when N t ≤ N t = N t increases further. Forboth cases, the positivity of c i is preserved unconditionally.Table 1: Accuracy test with initial values { c , c , c , c } = { } . N t { λ , λ } = { } { λ , λ } = (cid:8) × , 2 × (cid:9) Error Order
Error ∞ Order
Error Order
Error ∞ Order40 1.80E-6 2.84E-6 4.00E-3 6.19E-380 4.45E-7 2.01 7.02E07 2.02 7.67E-4 2.38 1.17E-3 2.41160 1.11E-7 2.00 1.75E-7 2.00 6.79E-5 3.50 9.78E-5 3.58320 2.77E-8 2.00 4.36E-8 2.00 6.89E-6 3.30 9.84E-6 3.31640 6.92E-9 2.00 1.09E-8 2.00 3.57E-6 0.95 6.00E-6 0.711280 1.73E-9 2.00 2.74E-9 1.99 1.27E-6 1.49 2.17E-6 1.472560 4.32E-10 2.00 6.96E-10 1.98 3.14E-7 2.02 5.33E-7 2.0314 preprint - J une
16, 2020 CH In this example, the following 1D detonation problem [42] is considered, CH + O → CO + H O .The reaction rate is defined as R = B T α H ( T − T ) ρ Y CH M CH (cid:18) ρ Y O M O (cid:19) , (66)where B = × , α = T = H ( x ) is the Heaviside function. M CH = M O = M CO =
44 and M H O = h CH =
500 and h O = h CO = h H O =
0. The specific heat ratio is set to beconstant as γ = T = p / ρ . The computational domain is [ −
25, 25 ] . Initial data are given by (cid:0) ρ , u , p , Y CH , Y O , Y CO , Y H O (cid:1) = (cid:40) (
2, 10, 40, 0, 0.2, 0.475, 0.325 ) , x ≤ − (
1, 0, 1, 0.1, 0.6, 0.2, 0.1 ) , x > − CFL ≤ .From time t n to t n + ∆ t , if negative pressure or density is found, the solved variables will be reinitializedby the solution at t n and the time step will be reduced to ∆ t . Such a time advancing technique keeps theCFL number as large as possible while guarantees the positivity when needed.The pressure, density and mass fraction of CH and O at t = N = N =
600 is enoughto resolve the location of right forward detonation wave. CH In this example, a radial symmetry 2D Detonation problem [43] is studied. The reaction model is the sameas in Section 6.2 except that h CH = [
0, 50 ] × [
0, 50 ] , where the lower andleft boundaries are inviscid walls, and the upper and right boundaries are outflow conditions. The initialvalues consist of totally burnt gas inside a circle with radius of 10 and unburnt gas outside, (cid:0) ρ , u , v , p , Y CH , Y O , Y CO , Y H O (cid:1) = (cid:16)
2, 10 xr , 10 yr , 40, 0, 0.2, 0.475, 0.325 (cid:17) , r ≤
10 , (
1, 0, 0, 1, 0.1, 0.6, 0.2, 0.1 ) , r >
10 , (68)where r = (cid:112) x + y .The computational domain consists of 600 ×
600 cells. Contours of pressure, density and radial velocityare shown in Fig. 2. Solutions along the line x = y are shown in Fig. 3. A detonation wave, a contactdiscontinuity and a shock wave can be observed with the decreasing of r . Discontinuities are capturedwithout oscillations and densities are advanced without negative values. Actually, for this problem, onlyone reaction involving two reactants is considered, the proposed Patankar scheme of Eq. (40) can be solvedanalytically without Newton iterations. The results are consistent with those of the work in [19]. H This problem involves 5 species and 2 reactions [43] which can be modeled by H + O → OH , 2 OH + H → H O . (69)15 preprint - J une
16, 2020 x Y CH
20 0 2000.040.080.12
ReferenceN=600 x Y O
20 0 200.20.40.6
ReferenceN=600
Figure 1: Solutions of 1D detonation problem of Sec. 6.2 at time t = N = N plays the role of a catalyst. The corresponding reaction rate is R = B T α H ( T − T ) ρ Y H M H ρ Y O M O , R = B T α H ( T − T ) ρ Y H M H (cid:18) ρ Y OH M OH (cid:19) , (70)where T = T = B = B = × , α = α = h H = h O = h N = h OH = − h H O = − M H = M O = M OH = M H O =
18 and M N = γ = T = p / ρ . Thecomputational domain is [ −
25, 25 ] and the initial values are (cid:0) ρ , u , p , Y H , Y O , Y OH , Y H O , Y N (cid:1) = (cid:40) (
2, 10, 40, 0, 0, 0.17, 0.63, 0.2 ) , x ≤ − (
1, 0, 1, 0.08, 0.72, 0, 0, 0.2 ) , x > − N =
600 are shown in Fig. (4) where the solid lines are reference valuesobtained with N = preprint - J une
16, 2020 (a) Pressure. (b) Density.(c) Radial velocity.
Figure 2: Contours of 2D detonation problem of CH at time t = ∆ x = H A 2D problem [43] with the same model set up as in Sec. 6.4 is presented in this section. The computationaldomain is [
0, 100 ] × [
0, 25 ] . Initial values are (cid:0) ρ , u , p , Y H , Y O , Y OH , Y H O , Y N (cid:1) = (cid:40) (
2, 10, 40, 0, 0, 0.17, 0.63, 0.2 ) , x ≤ ξ ( y ) , (
1, 0, 1, 0.08, 0.72, 0, 0, 0.2 ) , x > ξ ( y ) , (72)where ξ ( y ) = (cid:40) | y − | > − | y − | , | y − | ≤ N x × N y = × y = t = preprint - J une
16, 2020 x P r ess u r e (a) Pressure x D e n s i t y (b) Density x Y CH (c) Mass fraction of CH x Y O (d) Mass fraction of O x Y C O (e) Mass fraction of CO x Y H O (f) Mass fraction of H O Figure 3: Numerical solution of 2D detonation problem of CH at time t = x = y . CFL = 0.2.Cell size is ∆ x = t = t = To illustrate that the proposed scheme is able to preserve the positivity of pressure, a problem with a coupleof reversible reactions [9] is simulated. And one of two reactions is endothermic. Considering O ↔ O , (74)with N appears as a catalyst. The reaction rate is R f = k f (cid:18) ρ Y O M O (cid:19) (cid:18) ρ Y O M O + ρ Y O M O + ρ Y N M N (cid:19) , R b = k b (cid:18) ρ Y O M O (cid:19) (cid:18) ρ Y O M O + ρ Y O M O + ρ Y N M N (cid:19) , (75)where k f = CT − e − E / T , k b = k f /exp ( b + b log z + b z + b z + b z ) and z = T . The parametersare C = × , E = b = b = b = − b = − b = − M O = M O = M N = h O = × and18 preprint - J une
16, 2020 x P r ess u r e
20 0 20020406080100
ReferenceN=600 x D e n s i t y
20 0 2011.522.533.54
ReferenceN=600 x Y H
20 0 2000.040.08
ReferenceN=600 x Y O
20 0 2000.20.40.60.8
ReferenceN=600
Figure 4: Solutions of 1D detonation problem of H at time t = N = h O = h N =
0. In this case, the specific heat ratios for each specie are set as γ O =
53 , γ O = γ N =
75 .Pressure, temperature and internal energy are correlated according to the idea gas law, i.e., p = ρ R g T = ρ e γ − ,with the mixing rules defined as 1 M = Y O M O + Y O M O + Y N M N ,1 γ − = γ O − M O M + γ O − M O M + γ N − M N M , R g = RM ,and R = ( mol · K ) .Computational domain is [ −
1, 1 ] consisting of 4000 uniformly distributed cells. Initial values are (cid:0) u , p , ρ Y O , ρ Y O , ρ Y N (cid:1) = (cid:16)
0, 1000, ρ Y L , ρ Y L , ρ Y L (cid:17) , x ≤ (cid:16)
0, 1, ρ Y R , ρ Y R , ρ Y R (cid:17) , x > ρ Y L = × − , ρ Y L = × − , ρ Y L = × − , ρ Y R = × − , ρ Y R = × − , ρ Y R = × − .19 preprint - J une
16, 2020Figure 5: Solutions along cross-section of y = H at time t =
2. CFL = 0.4and grid number is N x × N y = × preprint - J une
16, 2020 (a) t =
1. (b) t = t =
3. (d) t = Figure 6: Density contours of 2D detonation problem of H from time t = N x × N y = × O and O in Fig.7. In Fig. 7, the reference values are solutions at the grid resolutionof ∆ x = N = N = In this study, the second order process based modified Patankar Runge-Kutta (PMPRK) schemes areproposed that maintain both the mass and mole in balance while retaining the positivity of density andpressure unconditionally for Euler equations with non-equilibrium reactive source terms. The sufficientand necessary conditions of the RK and Patankar coefficients are derived to fulfill the prior second orderof accuracy. The positivity of the proposed scheme is supported by the rigorous proof under simplifiedconditions and numerical validations with large randomly generated samples for general conditions.Coupled with the finite volume method, the PMPRK is extended to and tested on the Euler equations withmulti-reactions. Work on the schemes with third order of accuracy is ongoing and will be reported in thefuture. 21 preprint - J une
16, 2020 x P r ess u r e
1 0.5 0 0.5 110 ReferenceN=4000 x D e n s i t y
1 0.5 0 0.5 110 ReferenceN=4000 x T e m p e r a t u r e
1 0.5 0 0.5 1100002000030000
ReferenceN=4000 x Y O
1 0.5 0 0.5 10.10.150.20.25
ReferenceN=4000 x Y O
1 0.5 0 0.5 100.050.10.15
ReferenceN=4000 x Y N
1 0.5 0 0.5 10.740.760.780.80.82
ReferenceN=4000
Figure 7: Solutions of 1D reversible reactions of Oxygen at time t = N = N = preprint - J une
16, 2020
Acknowledgment
This work was supported by the computational resources allocated by the Ohio Supercomputer Center.
References [1] Qiang Zhou, Liang Zeng, and Liang-Shih Fan. Syngas chemical looping process: Dynamic modelingof a moving-bed reducer.
AIChE Journal , 59(9):3432–3443, 2013.[2] Qiang Zhou and Liang-Shih Fan. A second-order accurate immersed boundary-lattice Boltzmannmethod for particle-laden flows.
Journal of Computational Physics , 268:269–301, 2014.[3] Tim C Lieuwen.
Unsteady combustor physics . Cambridge University Press, 2012.[4] Hans Burchard, Eric Deleersnijder, and Andreas Meister. Application of modified Patankar schemesto stiff biogeochemical models for the water column.
Ocean Dynamics , 55(3-4):326–337, 2005.[5] Wilfried Kühn and Günther Radach. A one-dimensional physical-biological model study of the pelagicnitrogen cycling during the spring bloom in the northern North Sea (FLEX’76).
Journal of MarineResearch , 55(4):687–734, 1997.[6] Xiangxiong Zhang and Chi-Wang Shu. On positivity-preserving high order discontinuous Galerkinschemes for compressible Euler equations on rectangular meshes.
Journal of Computational Physics ,229(23):8918–8934, 2010.[7] Xiangxiong Zhang and Chi-Wang Shu. On maximum-principle-satisfying high order schemes forscalar conservation laws.
Journal of Computational Physics , 229(9):3091–3120, 2010.[8] Yulong Xing, Xiangxiong Zhang, and Chi-Wang Shu. Positivity-preserving high order well-balanceddiscontinuous Galerkin methods for the shallow water equations.
Advances in Water Resources ,33(12):1476–1493, 2010.[9] Xiangxiong Zhang and Chi-Wang Shu. Positivity-preserving high order finite difference WENOschemes for compressible Euler equations.
Journal of Computational Physics , 231(5):2245–2258, 2012.[10] Xiangxiong Zhang, Yinhua Xia, and Chi-Wang Shu. Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes.
Journal of Scientific Computing , 50(1):29–62, 2012.[11] Yifan Zhang, Xiangxiong Zhang, and Chi-Wang Shu. Maximum-principle-satisfying second orderdiscontinuous Galerkin schemes for convection–diffusion equations on triangular meshes.
Journal ofComputational Physics , 234:295–316, 2013.[12] Alina Chertock, Shumo Cui, Alexander Kurganov, and Tong Wu. Steady state and sign preservingsemi-implicit Runge–Kutta methods for ODEs with stiff damping term.
SIAM Journal on NumericalAnalysis , 53(4):2008–2029, 2015.[13] Alina Chertock, Shumo Cui, Alexander Kurganov, and Tong Wu. Well-balanced positivity preservingcentral-upwind scheme for the shallow water system with friction terms.
International Journal forNumerical Methods in Fluids , 78(6):355–383, 2015.[14] Juntao Huang and Chi-Wang Shu. A second-order asymptotic-preserving and positivity-preservingdiscontinuous Galerkin scheme for the Kerr–Debye model.
Mathematical Models and Methods in AppliedSciences , 27(03):549–579, 2017.[15] Jingwei Hu, Ruiwen Shu, and Xiangxiong Zhang. Asymptotic-preserving and positivity-preservingimplicit-explicit schemes for the stiff BGK equation.
SIAM Journal on Numerical Analysis , 56(2):942–973,2018. 23 preprint - J une
16, 2020[16] Hyun Geun Lee. Stability Condition of the Second-Order SSP-IMEX-RK Method for the Cahn–HilliardEquation.
Mathematics , 8(1):11, 2020.[17] David A Pope. An exponential method of numerical integration of ordinary differential equations.
Communications of the ACM , 6(8):491–493, 1963.[18] Juntao Huang and Chi-Wang Shu. Bound-preserving modified exponential Runge–Kutta discontinuousGalerkin methods for scalar hyperbolic equations with stiff source terms.
Journal of ComputationalPhysics , 361:111–135, 2018.[19] Jie Du and Yang Yang. Third-order conservative sign-preserving and steady-state-preserving time inte-grations and applications in stiff multispecies and multireaction detonations.
Journal of ComputationalPhysics , 395:489–510, 2019.[20] Suhas Patankar.
Numerical heat transfer and fluid flow . Taylor & Francis, 2018.[21] Hans Burchard, Eric Deleersnijder, and Andreas Meister. A high-order conservative Patankar-typediscretisation for stiff systems of production–destruction equations.
Applied Numerical Mathematics ,47(1):1–30, 2003.[22] Jorn Bruggeman, Hans Burchard, Bob W. Kooi, and Ben Sommeijer. A second-order, unconditionallypositive, mass-conserving integration scheme for biochemical systems.
Applied Numerical Mathematics ,57(1):36–58, 2007.[23] N. Broekhuizen, Graham J. Rickard, J. Bruggeman, and A. Meister. An improved and generalizedsecond order, unconditionally positive, mass conserving integration scheme for biochemical systems.
Applied Numerical Mathematics , 58(3):319–340, 2008.[24] Juntao Huang and Chi Wang Shu. Positivity-Preserving Time Discretizations for Production-Destruction Equations with Applications to Non-equilibrium Flows.
Journal of Scientific Computing ,78(3):1811–1839, 2019.[25] Juntao Huang, Weifeng Zhao, and Chi Wang Shu. A Third-Order Unconditionally Positivity-PreservingScheme for Production-Destruction Equations with Applications to Non-equilibrium Flows.
Journal ofScientific Computing , 79(2):1015–1056, 2019.[26] Stefan Kopecz and Andreas Meister. Unconditionally positive and conservative third order modifiedPatankar–Runge–Kutta discretizations of production–destruction systems.
BIT Numerical Mathematics ,58(3):691–728, 2018.[27] Philipp Öffner and Davide Torlo. Arbitrary high-order, conservative and positivity preservingPatankar-type deferred correction schemes.
Applied Numerical Mathematics , 153:15–34, 2020.[28] Ernst Hairer, Syvert P Nørsett, and Gerhard Wanner.
Solving ordinary differential equations. I, Nonstiffproblems . Springer-Vlg, 1991.[29] Sigal Gottlieb, David I Ketcheson, and Chi-Wang Shu.
Strong stability preserving Runge-Kutta andmultistep time discretizations . World Scientific, 2011.[30] Jianhua Pan. Patankar subroutine for PMPRKS. https://github.com/LucienPan0903/PatankarSubroutineForPMPRK , 2020.[31] Xiangxiong Zhang and Chi Wang Shu. On positivity-preserving high order discontinuous Galerkinschemes for compressible Euler equations on rectangular meshes.
Journal of Computational Physics ,229(23):8918–8934, 2010.[32] Xiangxiong Zhang and Chi Wang Shu. Positivity-preserving high order discontinuous Galerkinschemes for compressible Euler equations with source terms.
Journal of Computational Physics ,230(4):1238–1248, 2011.[33] B. Einfeldt, C. D. Munz, P. L. Roe, and B. Sjögreen. On Godunov-type methods near low densities.
Journal of Computational Physics , 92(2):273–295, 1991.24 preprint - J une
16, 2020[34] Benoit Perthame and Chi Wang Shu. On positivity preserving finite volume schemes for Eulerequations.
Numerische Mathematik , 73(1):119–130, 1996.[35] Benoıt Perthame. Second-order Boltzmann schemes for compressible Euler equations in one and twospace dimensions.
SIAM Journal on Numerical Analysis , 29(1):1–19, 1992.[36] Amiram Harten, Peter D. Lax, and Bram van Leer. On Upstream Differencing and Godunov-TypeSchemes for Hyperbolic Conservation Laws.
SIAM Review , 25(1):35–61, 1983.[37] Paul Batten, Nicholas Clarke, Claire Lambert, and Derek M Causon. On the choice of wavespeeds forthe HLLC Riemann solver.
SIAM Journal on Scientific Computing , 18(6):1553–1570, 1997.[38] Juan Cheng and Chi-Wang Shu. Positivity-preserving Lagrangian scheme for multi-material compress-ible flow.
Journal of Computational Physics , 257:143–168, 2014.[39] Zhen-Sheng Sun, Yu-Xin Ren, Cédric Larricq, Shi-ying Zhang, and Yue-cheng Yang. A class offinite difference schemes with low dispersion and controllable dissipation for DNS of compressibleturbulence.
Journal of computational physics , 230(12):4616–4635, 2011.[40] QiuJu Wang, YuXin Ren, ZhenSheng Sun, and YuTao Sun. Low dispersion finite volume schemebased on reconstruction with minimized dispersion and controllable dissipation.
Science China Physics,Mechanics and Astronomy , 56(2):423–431, 2013.[41] Qiuju Wang, Ralf Deiterding, Jianhua Pan, and Yu-Xin Ren. Consistent High Resolution Interface-Capturing Finite Volume Method for Compressible Multi-material Flows.
Computers & Fluids , page104518, 2020.[42] Weizhu Bao and Shi Jin. The random projection method for stiff multispecies detonation capturing.
Journal of Computational Physics , 178(1):37–57, 2002.[43] Wei Wang, Chi-Wang Shu, Helen C Yee, Dmitry V Kotov, and Björn Sjögreen. High order finite differ-ence methods with subcell resolution for stiff multispecies discontinuity capturing.