Moment Method for the Boltzmann Equation of Reactive Quaternary Gaseous Mixture
MMoment Method for the Boltzmann Equation of Reactive Quaternary GaseousMixture
Neeraj Sarna a , Georgii Oblapenko b , Manuel Torrilhon a a Center for Computational Engineering, Department of Mathematics, RWTH Aachen University, Schinkelstr 2, 52062, Germany b Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, 2201 E 24th St, Stop C0200,Austin, TX 78712
Abstract
We are interested in solving the Boltzmann equation of chemically reacting rarefied gas flows using the Grad’s-14moment method. We first propose a novel mathematical model that describes the collision dynamics of chemicallyreacting hard spheres. Using the collision model, we present an algorithm to compute the moments of the Boltzmanncollision operator. Our algorithm is general in the sense that it can be used to compute arbitrary order momentsof the collision operator and not just the moments included in the Grad’s-14 moment system. For a first-orderchemical kinetics, we derive reaction rates for a chemical reaction outside of equilibrium thereby, extending theArrhenius law that is valid only in equilibrium. We show that the derived reaction rates (i) are consistent in thesense that at equilibrium, we recover the Arrhenius law and (ii) have an explicit dependence on the scalar fourteenthmoment, highlighting the importance of considering a fourteen moment system rather than a thirteen one. Throughnumerical experiments we study the relaxation of the Grad’s-14 moment system to the equilibrium state.
Keywords:
Moment method, reacting flow, rarefied gas, Boltzmann equation
1. Introduction
The Boltzmann equation (BE) accurately models gas flows in all thermodynamic regimes. It governs theevolution of a probability density function that is defined on a seven dimensional space-time-velocity domain. Weconsider the BE of a gaseous mixture of four mono-atomic gases that react via a chemical reaction given as A + A (cid:28) A + A . (1)In a concise form, the BE for such a gaseous mixture can be written as d t f α ( c α , t ) = Q α ( f ( c , t ) , f ( c , t ) , f ( c , t ) , f ( c , t )) . (2)Above, α ∈ { , , , } is an index for the different gases in the mixture, f α represents the probability densityfunction of the α -th gas, and c α ∈ R is the molecular velocity. For simplicity, we restrict our study to a spatiallyhomogeneous BE. The collision operator Q α models the inter-particle interaction and contains contributions fromboth the chemical and the mechanical interactions. It is noteworthy that a collision between molecules that canundergo a chemical reaction ( A and A , for instance) does not necessarily result in a chemical reaction—the kineticenergy should be large enough to trigger a chemical reaction. Therefore, the part of Q α that models mechanicalcollisions contains interactions between all the gas molecules. The explicit form of Q α discussed later providesfurther clarification. Email addresses: [email protected] (Neeraj Sarna), [email protected] (Georgii Oblapenko), [email protected] (Manuel Torrilhon)
Preprint submitted to Elsevier September 8, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p e want to solve the above equation using a deterministic method i.e., with a method that does not requireMonte-Carlo sampling, which introduces an undesirable stochastic noise in the numerical approximation. Weconsider a Galerkin type approach where we approximate f α ( · , t ) in a span of the basis functions { φ α,i } i . Thisprovides the approximation f α ( · , t ) ≈ f α,N ( · , t ) = N (cid:88) i =0 α i ( t ) φ α,i f α . (3)We choose φ α,i as polynomials in the velocity variable c α and scale them with a Gaussian distribution function f α —the exact form of φ α,i and f α is discussed later in section 3. The above approximation was first proposed byGrad [1] and has several desirable properties, it (i) preserves the Galelian invariance of the BE, (ii) results in mass,momentum and energy conservation, (iii) at least in the linearized regime, converges to the BE [2, 3, 4], (iv) resultsin a hyperbolic moment system under appropriate regularization [5], etc. We refer to the review article [6] and thereferences therein for an exhaustive discussion.To compute the approximation f α,N , we replace f α by f α,N in the BE, multiply by the test functions { φ α,i } i and integrate over the velocity-domain R . This results in a time-dependent ordinary-differential-equation (ODE)for the expansion coefficients { α i } i , or the so-called moment equations. To well-define the moment equations, weneed integrals of the form I α ( t ) := (cid:90) R φ α,i ( c α , t ) Q α ( f ( c , t ) , f ( c , t ) , f ( c , t ) , f ( c , t )) dc α . (4)The main objective of the paper is to compute the above integral. To this end, we study the following threequestions related to chemically reactive collisions.1. Through which molecular potential do the molecules interact?2. For a given molecular potential and for binary collisions, how to relate the pre and the post collision relativevelocities using mass, momentum and energy conservation?3. For a given relation between the post and the pre collisional relative velocities, how to compute I α ( t )?Note that we only consider the contribution from binary collision in the collision operator Q α —a standard assump-tion in the kinetic gas theory [7]. Furthermore, all the above questions have been extensively studied in the contextof mechanical collisions, see [8, 9, 10, 7, 1].We use the Direction Simulation Monte Carlo (DSMC) method as a motivation to answer the first question andconsider a hard sphere interaction potential for both the chemical and the mechanical interactions [11]. Intuitively,the hard sphere potential treats molecules like billiard balls, which interact with each other only when they ”touch”each other. To answer the second question, we assume that the collisions are symmetric that is, the pre collisionalrelative velocity changes in equal proportions along the normal and the tangential direction of collision. We referto section 2 for a precise relation between the pre and the post collisional relative velocities. An answer to thethird question follows by extending the framework for a mono-atomic binary gas mixture (proposed in [10, 9]) to amixture of chemically reacting gases.To the best of our knowledge, only the work in [12] proposes an algorithm to compute the integral I α ( t ). Itdiffers from our work in the following sense. Authors in [12] accommodate the chemically reacting molecules atthe Maxwell-Boltzmann distribution function corresponding to the thermal equilibrium. This largely simplifies theexpression for I α ( t ), making the computations simpler. In contrast, we do not make any such assumption. Ourframework allows for chemical reactions that occur outside of a thermal equilibrium, which is of particular interestfor rarefied gas flows. Here, we emphasis that a thermal equilibrium is different from a chemical equilibrium. In a2hermal equilibrium, the probability density function is a Maxwell-Boltzmann distribution function but the numberdensities do not need to be related by the law of mass action, see [13, 12] or section 4 for details. Furthermore,authors in [12] restrict to a Grad’s-13 approximation i.e., to a particular choice of N in f α,N . The framework wepropose is not restricted to a particular N . However, for demonstration purposes, we focus particularly on theGrad’s-14 approximation.Let us mention that in DSMC computations, simple scattering models are used to simulate chemically reactivecollisions [11]. For instance, see the isotropic variable hard sphere model (VHS) or the anisotropic variable softsphere model (VSS) [14]. The DSMC literature focuses largely on incorporating effects of preferential dissociationfrom high-energy internal states into the chemical cross-section models [15]. Moreover, it has been argued that themicroscopic details of chemically reactive collisions do not have a significant influence on the flow dynamics [16].However, to realize/well-define a deterministic method, it is imperative to compute the integral I α ( t ), which requiresthe microscopic details of chemically reactive collisions. Our work proposes one possible way to define thesemicroscopic details.The BE of chemically reacting gases has not received much attention from a numerical approximation standpoint.Nevertheless, several works in the previous decades have contributed to its theoretical understanding. For the sakeof completeness, we briefly recall some of these works. One can show that similar to a single mono-atomic gas,the BE exhibits the H -theorem [17, 18]. Furthermore, the equation has a unique equilibrium that is stable andcan be given in terms of the Maxwell-Boltzmann distribution function [18]. Using the Chapman-Enksog expansion,one can also understand the behaviour of the transport coefficients of the gas mixture [19]. Furthermore, with ageneralized Chapman–Enksog method, one can provide viscous corrections to the non-equilibrium process rates andstudy the effect of flow compressibility on the reaction rate coefficients [20, 21].
2. Dynamics of reactive binary collisions
We start with discussing the dynamics of reactive binary collisions. The dynamics of mechanical binary collisionsis well studied in the literature [8] and is a special case of reactive collision—details included in remark 2. With m α we denote the mass of the molecule A α and with (cid:15) α we denote its energy of formation. We consider chemicalreactions of the type A + A (cid:28) A + A . With (cid:15) f and (cid:15) r we represent the activation energy of the forward andthe reverse reaction, respectively. The heat of the reactions is represented by Q and reads Q = (cid:15) f − (cid:15) r = ( (cid:15) + (cid:15) ) − ( (cid:15) + (cid:15) ) . For brevity, we only consider the forward reaction—the relations for the reverse reaction are similar. Throughoutthe article, we ignore the angular momentum of the molecules. Consider two molecules with masses m and m travelling with velocities c and c , respectively. The molecules collide to form two molecules with masses m and m , respectively, and with the post collisional velocities c and c , respectively. We define the pre and the postcollisional relative velocities as g = c − c , g = c − c . (5)We assume that A and A have enough kinetic energy to trigger a chemical reaction i.e., m g ≥ (cid:15) f , m g ≥ (cid:15) r , (6)where m αβ = m α m β m α + m β . 3ince mass, momentum and energy are conserved in a binary collision, we find [13] m + m = m + m , m c + m c = m c + m c , m | g | m | g | Q. (7)Above, |·| represents the Eucledian norm of a vector. For later convenience, we use the energy balance to derivethe following two relations | g | = ˆ Q g , | g | = ˆ Q | g | , (8)where ˆ Q = m m (cid:18) − Qm g (cid:19) , ˆ Q = m m (cid:18) Qm g (cid:19) . (9) Remark 1.
The constraint on the kinetic energies to trigger a chemical reaction given in (6) does not take intoaccount the geometry of collision—see [13, 22, 23]. Taking these geometrical effects into account is beyond the scopeof the present work and we hope to discuss them in one of our future works.2.2. Hard sphere interaction potential
With r αβ we represent the distance between the centres of the two molecules A α and A β . We assume that thesetwo molecules interact via a molecular potential Ψ αβ ( r αβ ) that depends solely on r αβ . We consider a hard sphereinteraction potential given as Ψ αβ ( r αβ ) = (cid:40) ∀ r αβ > r min αβ ∞ ∀ r αβ ≤ r min αβ , (10)where r min αβ is the minimum distance between the two molecules. When molecules undergo a mechanical collision, r min αβ is the average of the molecular diameters d α and d β i.e., r min αβ = ( d α + d β ) . For the molecules that undergoa chemical collision, we define r min αβ separately for the forward and the reverse reaction as r min12 ,f = d f , r min34 ,r = d r . (11)Here, d f and d r are the reactive diameters for the forward and the reverse reactions, respectively. Similar to [24],we assume that d f and d r are linearly related to d and d in the following way d f = s f d and d r = s r d , where s f and s r are the so called steric factors—a relation between s f and s r follows from the law of mass action and isgiven later in (70).For clarity of the following discussion, we visualize a collision between two hard spheres in Figure 1. The vector k is a unit vector (i.e., | k | = 1) that passes through the line joining the centres of the two interacting molecules and b is the so-called impact parameter. The angles θ and χ are the angle-of-attack and the scattering angle, respectively. The aim of the following discussion is to express the post-collisional velocities c and c in terms of the pre-collisional c and c ones and vice-versa. We will use these relations later to find the moments of the collisionoperators given in (2). We start with relating the relative velocity g to g .4 igure 1: Binary collision for a hard sphere interaction potential. g : relative pre-collisional velocity, g : relative post-collisionalvelocity, b : impact parameter, k : a unit vector that passes through the centres of the two colliding spheres, e n : normal direction ofcollision, and e t : tangential direction of collision. Let e n and e t be two orthonormal unit-vectors that define the normal and the tangential direction of collision,respectively. Both of these vectors are shown in Figure 1 and can be related to the vector k via e n = − k and e t = ( − k ) ⊥ . Let g n = ( g · e n ) e n and g t = ( g · e t ) e t . Here, u · v represents an Eucledian inner-product betweentwo vectors u and v . Similarly, define g n and g t . Since the span of e n and e t approximates every vector in thecollision plan, using the orthogonality of e n and e t , we find g = g n + g t , g = g n + g t . (12)Let E n and E t be the coefficients of restitution along the normal and the tangential direction, respectively, definedas g n = E n g n , g t = E t g t . (13)We want to find a relation for E n and E t in terms of the pre/post collisional velocities g /g and the heat ofthe reaction Q . One possibility is to use the Hertzian contact theory of rough inelastic hard spheres [25]. However,involvement of chemical reactions makes the application of Hertzian contact theory complicated and, as yet, it isunclear how to extend this contact theory to chemically reacting molecules. In the present work, we propose tocompute E n and E t by assuming that the pre-collisional relative velocity g changes in equal proportions along thenormal and the tangential direction of collision. Equivalently, E n = E t = E . (14)To compute E , we first note that | g | = E | g | . Then, the energy balance given in (8) provides E = (cid:113) ˆ Q . (15)5ote that we have only considered the positive root in the above expression because E n cannot be negative. Alsonote that ˆ Q ≥ Q ,and (ii) a change in mass of the interacting molecules. As per the above assumption, both these effects influence thepre collisional relative velocity g equally along the normal and the tangential direction of collision. Equivalently,there is no preferred direction of motion along which a chemical reaction changes the pre-collisional relative velocitydifferently.A collision between two rough inelastic hard spheres also results in a change in both the tangential and thenormal pre-collisional relative velocity. We discuss the differences of our collision model to the collision model of arough inelastic hard sphere. We refer to [26] for a detailed discussion on a rough inelastic hard sphere.1. For an inelastic rough hard sphere, two different mechanisms change the pre-collisional relative velocity—friction changes the tangential velocity and the inelasticity changes the normal velocity. Since both these mech-anisms are different, usually, E n (cid:54) = E t . In contrast, for a chemically reactive collision, a single mechanism—thechemical reaction—is responsible for a change in both the normal and the tangential pre-collisional velocity,which motivates our assumption E n = E t .2. For a rough inelastic hard sphere, if the friction is sufficiently large then E t can be negative, resulting in achange in both the tangential and the normal pre-collisional velocity. In our model E t is positive thus, thepre-collisional relative velocity can only change in magnitude. We now express c /c in terms of c /c . We first summarise our results, and devote the rest of the discussionin proving it. The post-collisional velocities c / are related to the pre-collisional ones c / via c = c + g (cid:18) µ (cid:113) ˆ Q − µ (cid:19) − µ (cid:113) ˆ Q ( k · g ) k, (16) c = c − g (cid:18) µ (cid:113) ˆ Q − µ (cid:19) + 2 µ (cid:113) ˆ Q ( k · g ) k. (17)Above, k is the unit vector shown in Figure 1 and ˆ Q is as given in (8). Furthermore, µ αβ is the mass fractiondefined as µ αβ = m α / ( m α + m β ) . (18)Expressions for c and c can be derived in a similar fashion and are given as c = c + g (cid:18) µ (cid:113) ˆ Q − µ (cid:19) + 2 µ (cid:113) ˆ Q ( k · g ) k, (19) c = c − g (cid:18) µ (cid:113) ˆ Q − µ (cid:19) − µ (cid:113) ˆ Q ( k · g ) k. (20)We now derive the above relations. First recall that the angles θ , φ and χ are as shown in Figure 1. Since weassume that the collision is symmetric (i.e., (14) holds) we find that (i) θ = φ , (ii) χ = π − θ , and (iii) k is a vector6hat bisects g and g and is given as k = (cid:18)(cid:113) ˆ Q g − g (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:113) ˆ Q g − g (cid:12)(cid:12)(cid:12)(cid:12) . (21)Taking a dot product of the above equation with (cid:113) ˆ Q g − g , we obtain the following relation k · (cid:18)(cid:113) ˆ Q g − g (cid:19) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:113) ˆ Q g − g (cid:12)(cid:12)(cid:12)(cid:12) (22)Since k bisects g and g , we also have k · g = −| g | cos( φ ) = − (cid:113) ˆ Q | g | cos( θ ) = − (cid:113) ˆ Q k · g . (23)Substituting (23) into (22) we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:113) ˆ Q g − g (cid:12)(cid:12)(cid:12)(cid:12) = 2 (cid:113) ˆ Q k · g (24)Substituting the above relation into the definition of k given in (21) we have g = (cid:113) ˆ Q g − (cid:113) ˆ Q ( k · g ) k. (25)The velocity of the center of mass of the system, namely h , is given as h = µ c + µ c = µ c + µ c (26)Note that the second equality is a result of the momentum conservation given in (7). Using the center of massvelocity h and the post collisional relative velocity g , the post-collisional velocities c and c can be written as c = h + µ g c = h − µ g , (27)which is equivalent to c = µ c + µ c + µ g , c = µ c + µ c − µ g . (28)Finally, using (25) in the above expression, we express g in terms of g to find the desired result. Remark 2.
For consistency, we show that the velocity transformations for a single mono-atomic gas and a binarymixture of (chemically neutral) mono-atomic gases follow from the velocity transformations of a chemically reactivegas mixture. For a single gas and for all α ∈ { , , , } , we have m α = m and Q = 0 . As a result, the velocitytransformations (16) – (20) , change to c = c − ( k · g ) k and c = c + ( k · g ) k . For a binary mixture, wehave m = m , m = m , and Q = 0 . Using these physical parameters, the velocity transformations change to c = c − µ ( k · g ) k and c = c + 2 µ ( k · g ) k . One can easily check that these relations are the same asthose given in [8, 9].
3. Moment Approximation
We present our moment approximation for the BE given in (2). We start with giving the explicit form of thecollision operator appearing in the BE (2). 7 .1. Collision operator
We split the collision operator Q α as Q α ( f , f , f , f ) = (cid:88) β =1 Q M ( f α , f β ) + Q α,R ( f , f , f , f ) . (29)Above, Q M ( f α , f β ) models the mechanical collisions between A α and A β , and Q α,R models the chemically reactivecollisions. The expression for Q M is the same as that for a single gas and reads [8] Q M ( f α , f β ) = (cid:90) R (cid:90) π (cid:90) π (cid:16) f (cid:48) α f (cid:48) β − f α f β (cid:17) g αβ σ αβ sin( χ ) dχd(cid:15)dc β , (30)where f (cid:48) α and f (cid:48) β represent the phase densities corresponding to the post collisional velocities c (cid:48) α and c (cid:48) β , respectively, (cid:15) represents the angle made by the collisional plane, and σ αβ represents the differential cross-section. For a hardsphere interaction potential, σ αβ is given as [8, 10] σ αβ = 14 d αβ , d αβ = 12 ( d α + d β ) . The operator Q α,R that accounts for chemical reactions can be given as [11, 13] Q α,R ( f , f , f , f ) dc α = ν α (cid:90) R (cid:90) R (cid:90) π (cid:90) π f f g σ f sin( χ ) dχd(cid:15)dc dc − ν α (cid:90) R (cid:90) R (cid:90) π (cid:90) π f f g σ r sin( χ ) dχd(cid:15)dc dc , (31)where ν α is a stoichiometric coefficient given as − ν = − ν = ν = ν = 1 . Furthermore, σ f and σ r are thedifferential cross-sections corresponding to the forward and the reverse reaction, respectively. For a hard sphereinteraction potential, both these cross-sections read σ f = U (cid:18) − (cid:15) f m g (cid:19) d f (cid:18) − (cid:15) f m g (cid:19) , (32) σ r = U (cid:18) − (cid:15) r m g (cid:19) d r (cid:18) − (cid:15) r m g (cid:19) . (33)Above, U ( x ) is a unit-step function that is one for x ≥ d f and d r are as givenin (11). Let i ∈ R N denote a multi-index with each entry being a natural number smaller than the velocity spacedimension d . Using i and a scalar a ∈ N , we define a polynomial over R d as φ a, (cid:104) i ...i N (cid:105) ( y ) = | y | a y (cid:104) i . . . y i N (cid:105) , y ∈ R d . (34)The tensor y (cid:104) i . . . y i N (cid:105) represents the trace free part of y i · · · y i N and its explicit form can be found in [7]. Forcompleteness, we give this explicit form for N = 2 B (cid:104) i i (cid:105) = 12 ( B i i + B i i ) − δ i i B i k i k , (35)8here δ i i is a Kronecker delta. Throughout the article, we use the Einstein’s summation convention for the tensornotation.Using φ a, (cid:104) i ...i N (cid:105) , we define a general (2 a + N )-th order moment of f α as w αa,i ...i N ( t ) = m α (cid:90) R d φ a, (cid:104) i ...i N (cid:105) ( C α ) f α ( c α , t ) dC α . (36)Above, C α is the so-called peculiar velocity given as C α = c α − v . Here, v is the macroscopic velocity of the gasmixture defined in the following way. Let the density ρ α and the macroscopic velocity v α of gas- α be defined as ρ α = w α , , ρ α v αi = w α ,i . (37)We define the mixture velocity v and the mixture density ρ as ρ = (cid:88) α =1 ρ α , ρv = (cid:88) α =1 ρ α v α . (38)The other moments of f α that have a physical relevance are given as32 kT α = 32 ρ α θ α = 12 w α , σ αi i = w α ,i i , q αi = 12 w α ,i . (39)Above, θ α = kT α /m α is the temperature in energy units with k being the Boltzmann’s constant, σ ( α ) i i is thestress-tensor and q ( α ) i is the heat flux. For convenience, we further define∆ α = w α , − ρ α θ α . We also define the temperature θ and the number density n of the gas mixture as kT = ρθ = (cid:88) i =1 ρ α θ α , n = (cid:88) i =1 n α , where n α = ρ α /m α . Remark 3.
Note that due to mass, momentum and energy conservation, the mixture velocity v is constant overtime resulting in a time-independent peculiar velocity C α .3.3. The fourteen moment system We are interested in solving for only the 14 moments contained in the set w [14] α = { w α , , w α ,i , w α , , w α ,i i , w α ,i , w α , } . Note that as per the above relations, this set contains one component for ρ α , three for v α , one for θ α , six for σ α ,three for q α and one for w α , —thus fourteen in total. Given the moment vector w [14] α , we first want to approximatethe distribution function f α . We consider the Grad’s-14 (G14) moment approximation and approximate f α by91, 26] f α ≈ f α | G = f α (cid:34) α ρ α θ α (cid:18) − | C α | θ α + 115 | C α | θ α (cid:19) + q ( α ) i C ( α ) i ρ α θ α (cid:18) | C α | θ α − (cid:19) + σ ( α ) i i ρ α θ α C ( α ) i C ( α ) i + u ( α ) i C ( α ) i θ α (cid:35) . (40)The function f α is a Gaussian-distribution function defined as f α = n α (cid:18) πθ α (cid:19) d exp (cid:18) − | C α | θ α (cid:19) . (41)To derive a governing equation for the moments w [14] α , we replace f α by f α | G in the BE (2), multiply theresulting equation by appropriate functions φ a, (cid:104) i ...i N (cid:105) ( C α ) and integrate over the velocity domain R d . For all α ∈ { , , , } , this results in the moment equations d t n α = P αR, , d t (cid:0) n α v αi (cid:1) = P αM,i + P αR,i , d t ( n α T α ) = P αM, + P αR, , d t σ ( α ) i i = P αM,i i + P αR,i i ,d t q ( α ) i = P αM, ,i + P αR, ,i , d t u α ) = P αM, + P αR, . (42)The different P α ’s are the moments of the collision operator defined as P αM,a,i ...i N = m α (cid:88) β =1 (cid:90) R d φ a, (cid:104) i ...i N (cid:105) Q M ( f α | G , f β | G ) dC α , P αR,a,i ...i N = m α (cid:90) R d φ a, (cid:104) i ...i N (cid:105) Q α,R ( f | G , f | G , f | G , f | G ) dC α . (43)Above, Q M and Q α,R are the mechanical and the chemical parts of the collision operator given in (29). Note thatthe mechanical collisions do not contribute into the mass balance equation. Remark 4.
Strictly speaking, the Grad’s-14 moment approximation is a special case of the approximation f α | G given in (40) . In f α , if we replace the ”local” gas temperature θ α by the temperature of the mixture θ , we find theGrad’s-14 moment approximation. We expect f α | G to be more accurate than a Grad’s-14 moment approximationbecause it adapts better to the changes in the ”local” temperature of the gas- α . Remark 5.
The following two observations motivate our choice of a fourteen moment system rather than a thirteenone. Firstly, out of all the fourteen moments, only the non-equilibrium moment ∆ α (and not the stress tensor andthe heat flux) influences the reaction rates, which can result in a better approximation of the reaction rates outsideof the equilibrium—see (64) and (65) derived later for explicit reaction rates. We emphasise that for a thirteenmoment system, none of the non-equilibrium moments will appear in the reaction rates. Secondly, as discussedearlier, we treat a chemically reacting mixture similar to a granular gas in the sense that we allow both the normaland the tangential pre-collisional velocities to change. Previous works have shown that including ∆ α in the set ofmoments better captures the cooling rate and the transport coefficients of a granular gas [26, 27]. We expect a similarresult to hold true for a chemically reactive mixture. Our numerical experiments corroborate our expectations bydemonstrating that non-equilibrium states where ∆ α is initially zero can trigger non-zero perturbations in ∆ α , whileall the other non-equilibrium moments remain zero.3.4. Moments of the collision operator To compute the moments of Q M we first derive moments of Q α,R and interpret the mechanical collisions as aspecial case of the chemical collisions by using the physical parameters given in remark 2. This results in the samemoments as those given in [9, 10]. For brevity, we do not discuss these moments here again.10o compute the moments of Q α,R , we restrict ourselves to α = 1, since computations for all the other gascomponents are the same. First, we decompose Q ,R into a gain Q g ,R and a loss part Q f ,R as Q ,R ( f | G , f | G , f | G , f | G ) dC = (cid:90) f | G f | G g σ r sin( χ ) dχd(cid:15)dC dC (cid:124) (cid:123)(cid:122) (cid:125) Q g dC − (cid:90) f | G f | G g σ r sin( χ ) dχd(cid:15)dC dC (cid:124) (cid:123)(cid:122) (cid:125) = Q l dC . (44)It is easy to conclude that Q g and Q l results in a gain and a loss in f , respectively. Note that because f | G isexpressed in terms of C , we have changed the integration variable from c to C . To perform the change we haveused the relation dc = dC , which immediately follows from the definition of C .To simplify the following notations, rather than computing the moments of Q α,R with respect to the polynomial φ a, (cid:104) i ...i N (cid:105) given in (34), we compute the moments with respect to the monomial φ i ...i N defined as φ i ...i N ( C α ) = C αi . . . C αi N . (45)Using the definition of C (cid:104) i . . . · C i N (cid:105) and expressing | C α | a in terms of the different components of C α , one canconclude that for any a ∈ N , the function φ a, (cid:104) i ...i N (cid:105) ( C α ) can be recovered via linear combinations of φ i ...i N ( C α ),making it straightforward to express (cid:82) R d φ a, (cid:104) i ...i N (cid:105) Q α,R dC α in terms of (cid:82) R d φ i ...i N Q α,R dC α .We define the moments of Q g and Q l with respect to φ i ...i N ( C ) as P gi ...i N = m (cid:90) φ i ...i N ( C ) Q g dC = m (cid:90) φ i ...i N ( C ) f | G f | G g σ r sin( χ ) dχd(cid:15)dC dC , P li ...i N = m (cid:90) φ i ...i N ( C ) Q l dC = − m (cid:90) φ i ...i N ( C ) f | G f | G g σ f sin( χ ) dχd(cid:15)dC dC . (46)Above, for brevity we have represented all the integrals with a single integration symbol. To compute P gi ...i N , wefollow the steps outlined in algorithm 1—the details of the algorithm are given thereafter. In the first step, weexpress φ i ...i N ( C ) in terms of the integration variables C and C . To compute P li ...i N , this transformation isnot needed. All the other steps remain the same as that for P gi ...i N therefore, we only discuss the computation of P gi ...i N . Note that the production terms resulting from algorithm 1 are included in the supplementary material. Algorithm 1
Algorithm to compute P gi ...i N Using the velocity transformations outlined in subsection 2.3, express φ i ...i N ( C ) in terms of C and C . Perform integration over the angles that determine the collision orientation i.e., the angles χ and (cid:15) appearingin (46). This results in an integral over C and C . Define new integration variables in terms of C and C such that the integral in the previous step can be writtenin terms of two separate integrals, with each integral expressed solely in terms of a single integration variable. With a method-of-choice, perform integration over the two separate integrals resulting from the previous step.We use the symbolic integration from mathematica to compute these integrals.
Remark 6.
Note that algorithm 1 can be used to compute an arbitrary order moment of the collision operator andnot just the fourteen moments. However, for the sake of demonstration, during numerical experiments and in theproduction terms outlined in the supplementary material, we stick to the Grad’s-14 moment approximation. .4.1. Step-1: Transformation of the test function Let (cid:36) = µ (cid:113) ˆ Q , where µ is the mass fraction defined in (18) and ˆ Q is as given in (8). We express thevelocity transformations given in (19) as C (1) i r = C (3) i r + ( (cid:36) − µ ) g (34) i r − (cid:36) g cos θk i r ( ∵ k · g = | g | cos( θ )) . (47)Let d = ( (cid:36) − µ ) g − k(cid:36) | g | cos θ . Then, the monomial φ i ...i N given in (45) can be expressed as φ i ...i N = m ( C (3) i + d i )( C (3) i + d i )( C (3) i + d i ) . . . ( C (3) i N + d i N ) , = m N (cid:88) β =0 (cid:18) Nβ (cid:19) d ( i d i d i . . . d i β C (3) i β . . . C (3) i N ) . (48)We simplify d i d i d i . . . d i β as d i d i . . . d i β = (( (cid:36) − µ ) g (34) i − (cid:36) g cos θk i )(( (cid:36) − µ ) g (34) i − (cid:36) g cos θk i ) . . . (( (cid:36) − µ ) g (34) w β − (cid:36) g cos θk i β ) (49)= β (cid:88) β =0 (cid:18) β β (cid:19) ( − β (2 g cos θ ) β (cid:36) β ( (cid:36) − µ ) β − β k ( i k i . . . k i β g (34) i β . . . g (34) i β ) . (50)Above and elsewhere, B ( i ...i r ) represents the symmetric part of a tensor B . For r = 2, B ( i ...i ) = ( B i i + B i i ) / r >
2, see [7]. Substituting the above simplification (50) into the expression for φ i ...i N (48), we find φ i ...i N = m N (cid:88) β =0 β (cid:88) β =0 (cid:18) Nβ , β − β (cid:19) ( − β (2 g cos θ ) β (cid:36) β × ( (cid:36) − µ ) β − β k ( i k i . . . k i β g i β . . . g (34) i β C (3) i β . . . C (3) i N ) . (51)In deriving the above expression, the following trivial property of the tensors has been used A (( i ...i r ) B i r +1 ...i N ) = A ( i ...i r B i r +1 ...i N ) . With the above expression for φ i ...i N , the expression for P gi ...i N transforms to P gi ...i N =4 m N (cid:88) β =0 β (cid:88) β =0 (cid:18) Nβ , β − β (cid:19) ( − β β × (cid:90) (cid:32)(cid:90) π (cid:90) π sin θ (cos θ ) β +1 k ( i k i . . . k i β dθd(cid:15) (cid:33) × g (34) i β . . . g (34) i β C (3) i β . . . C (3) i N ) f | G f | G σ r g β +134 (cid:36) β ( (cid:36) − µ ) β − β dC dC . (52) We first compute the underlined integral appearing above in (52). Since the computation is the same as thatdiscussed in [7], we only discuss the main results. Consider a β th order tensor F i i ...i β such that F i i ...i β = (cid:90) π (cid:90) π k i k i . . . k i β sin θ (cos θ ) β +1 dθd(cid:15) (53)12e can express F as [9] F i i ...i β = β / (cid:88) γ =0 a ( β ) γ δ ( i i . . . δ i γ − i γ g (34) i γ +1 g . . . g (34) i β ) g , (54)where, for a hard sphere interaction potential, the coefficients a β γ depend only on g . Replacing the underlinedterm in P gi ...i N (52) by the above expression, we find P gi ...i N = 4 m N (cid:88) β =0 β (cid:88) β =0 (cid:18) Nβ , β − β (cid:19) ( − β β β / (cid:88) γ =0 a ( β ) γ (cid:90) δ ( i i . . . δ i γ − i γ g (34) i γ +1 . . . g (34) i β × g (34) i β . . . g (34) i β C (3) i β . . . C (3) i N ) f | G f | G σ r g γ +134 (cid:36) β ( (cid:36) − µ ) β − β dC dC . (55)With the above expression, to complete our computation of P gi ...i N , we only need to perform integrals over C and C . The following two steps take care of this. To make the following discussion simpler, we perform non-dimensionalizationˆ σ r = U (cid:18) − (cid:15) f ˆ g (cid:19) d f (cid:18) − (cid:15) f ˆ g (cid:19) , ˆ σ r = U (cid:18) − (cid:15) r ˆ g (cid:19) d r (cid:18) − (cid:15) r ˆ g (cid:19) , ˆ f α = f α θ d α n α , ˆ C α = C α √ θ α , ˆ g αβ = g αβ (cid:112) θ αβ , ˆ (cid:15) r = (cid:15) r m θ , ˆ (cid:15) f = (cid:15) f m θ ˆ (cid:36) = µ (cid:118)(cid:117)(cid:117)(cid:116) m m (cid:32) Q ˆ g (cid:33) , ˆ Q = Qm θ , (56)where θ αβ = ( θ α + θ β ) /
2. Using the above scaling, P gi ...i N given in (55) reads P gi ...i N = 4 m n n N (cid:88) β =0 β (cid:88) β =0 (cid:18) Nβ , β − β (cid:19) ( − β β β / (cid:88) γ =0 a ( β ) γ (cid:90) δ ( i i . . . δ i γ − i γ ˆ g (34) i γ +1 . . . ˆ g (34) i β × ˆ g (34) i β . . . ˆ g (34) i β ˆ C (3) i β . . . ˆ C (3) i N ) ˆ f | G ˆ f | G ˆ σ r ˆ g γ +134 × θ β θ n − β ˆ (cid:36) β ( ˆ (cid:36) − µ ) β − β d ˆ C d ˆ C . We want to define new integration variables in terms of ˆ C and ˆ C such that above integral (57) can be separatedout into two separate integrals. The structure of the Gaussian f α in (40) motivates the transformationˆ h = 12 (cid:32)(cid:114) θ θ ˆ C + (cid:114) θ θ ˆ C (cid:33) , ˆ g = (cid:114) θ θ ˆ C − (cid:114) θ θ ˆ C . (57)Note that the above relation provides ( | ˆ C | + | ˆ C | ) = (cid:16) | ˆ g | + | ˆ h | (cid:17) . Furthermore, one can show that thedeterminant of the Jacobian of transformation given as det( J ) = (cid:12)(cid:12)(cid:12)(cid:12) det (cid:18) ∂ ( ˆ C , ˆ C ) ∂ (ˆ h , ˆ g ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) , equals one. The above trans-13ormation of variables provides P gi ...i N = (cid:88) p,q ≤ N δ ( p,q ) (cid:90) R d | ˆ h | n ˆ h (34) i . . . ˆ h (34) i p exp (cid:104) −| ˆ h | (cid:105) d ˆ h (cid:124) (cid:123)(cid:122) (cid:125) I h × (cid:90) R d | ˆ g | m (cid:36) l ˆ σ r ˆ g (34) i . . . ˆ g (34) i q exp (cid:20) − | ˆ g | (cid:21) d ˆ g (cid:124) (cid:123)(cid:122) (cid:125) I g . (58)In the above expression, δ ( p,q ) is some factor that results from the non-dimensionalization given in (56), its preciseform is not important here. From the above expression (58) it is clear that we have separated out the integrationover ˆ h and ˆ g . We now discuss the computation of I h and I g . I g and I h Expressing ˆ h in terms of spherical harmonics, I h can be given as [9] I h = 12 Γ (cid:18) N + p + 32 (cid:19) πp + 1 δ i ...i p (59)where δ i ...i p is the fully symmetric product of Kronecker deltas and Γ( · ) denotes a gamma-function. Similarly,transforming ˆ g using spherical coordinates, I g can be expressed as I g = (cid:90) π (cid:90) π ˜ n i . . . ˜ n i p sin νdνdω (cid:90) ∞ (cid:36) l ˆ σ r | ˆ g | m + q +2 exp (cid:20) − | ˆ g | (cid:21) d | ˆ g | , (60)where ˜ n i = ˆ g (34) i / ˆ g . Computing the underlined term in the above expression, we find [9] I g = 4 πp + 1 δ i ...i p (cid:90) ∞ (cid:36) l ˆ σ r | ˆ g | m + p +2 exp (cid:20) − | ˆ g | (cid:21) d | ˆ g | . (61)Inserting the expression for ˆ σ r in the above expression, we find that the underlined integral in the above expressionis of the form Ω ( l,r ) (cid:15) r = (cid:90) ∞ U (cid:18) − (cid:15) r | ˆ g | (cid:19) | ˆ g | r ˆ (cid:36) l exp (cid:20) − (cid:18) | ˆ g | (cid:19)(cid:21) d | ˆ g | . We first consider the case when l = 0. The definition of an incomplete Gamma function Γ ( a, b ) providesΩ (0 ,r ) (cid:15) r = 2 r Γ (cid:16) r , ˆ (cid:15) r (cid:17) . For l (cid:54) = 0, we simplify Ω ( l,r ) (cid:15) r in the following way and compute it explicitly using thesymbolic integration in mathematica Ω ( l,r ) (cid:15) r = (cid:90) ∞ U (cid:18) − (cid:15) r | ˆ g | (cid:19) | ˆ g | r ˆ (cid:36) l exp (cid:20) − | ˆ g | (cid:21) d | ˆ g | = (cid:90) ∞ U (cid:18) − (cid:15) r | ˆ g | (cid:19) | ˆ g | r µ l (cid:18) m m m m (cid:18) (cid:15) f − ˆ (cid:15) r ) | ˆ g | (cid:19)(cid:19) l exp (cid:20) − | ˆ g | (cid:21) d | ˆ g | = (cid:90) ∞ U (cid:18) − (cid:15) r | ˆ g | (cid:19) | ˆ g | ( r − l ) µ l (cid:18) m m m m (cid:0) | ˆ g | +2 (ˆ (cid:15) f − ˆ (cid:15) r ) (cid:1)(cid:19) l exp (cid:20) − | ˆ g | (cid:21) d | ˆ g | . (62)The second equality is a result of energy conservation given in (8). The different Ω ( l,r ) (cid:15) r integrals used in the presentwork are given in the supplementary material. 14 .5. Non-equilibrium reaction rates Assuming that the chemical reaction (1) follows first-order chemical kinetics, the forward k f and the backward k r reaction rates can be given as k f = 1 n n (cid:90) f f σ f g sin χdχdc dc , k r = 1 n n (cid:90) f f σ r g sin χdχdc dc . (63)Replacing the probability density functions by their Grad’s-14 approximation and using technique outlined in 1, wecan compute the approximations k f | G and k r | G . A concise form of both these rates can be given as k f | G = 2 (cid:112) π ( θ + θ ) ( d s f ) exp (cid:20) − (cid:18) (cid:15) f m ( θ + θ ) (cid:19)(cid:21) + ˜∆ f , (64) k r | G = 2 (cid:112) π ( θ + θ ) ( d s r ) exp (cid:20) − (cid:18) (cid:15) r m ( θ + θ ) (cid:19)(cid:21) + ˜∆ r , (65)where ˜∆ f and ˜∆ r represent the contribution from the non-equilibrium moment ∆ α and are given as˜∆ r = β r ∆ θ µ (cid:0) − θ m − θ m (cid:15) r + (cid:15) r (cid:1) m n n + β r ∆ θ µ (cid:0) − θ m − θ m (cid:15) r + (cid:15) r (cid:1) m n n , ˜∆ f = β f ∆ θ µ (cid:16) − θ m − θ m (cid:15) f + (cid:15) f (cid:17) m n n + β f ∆ θ µ (cid:16) − θ m − θ m (cid:15) f + (cid:15) f (cid:17) m n n . (66)With M = m + m = m + m , the coefficients β f and β r are given as β f = (cid:32) √ πd M n n s f θ / m (cid:33) exp (cid:20) − (cid:15) f θ m (cid:21) , β r = (cid:32) √ πd M n n s r θ / m (cid:33) exp (cid:20) − (cid:15) r θ m (cid:21) . (67)Furthermore, s f and s r represent the steric factors defined in (11). Note that, out of all the fourteen moments inthe set w [14] α , ∆ α is the only non-equilibrium moment (i.e., a moment that vanishes in equilibrium) that influencesthe reaction rates. This highlights the importance of considering a 14-moment system instead of the 13 one. Weexpect that the presence of ∆ α provides a better approximation to the true reaction rates outside of equilibrium.Our expressions for k f | G and k r | G are consistent in the sense that in equilibrium, we recover the Arrheniuslaw [24]. The details of the computations are as follows. In a chemical equilibrium, the following two conditions hold(i) all the gases have the same temperature i.e., T = T = T = T = T eq , and (ii) the non-equilibrium moment∆ α vanishes, providing ˜∆ f = ˜∆ r = 0—see [18, 12] for a detailed study of the equilibrium state. Substituting thesetwo conditions into (64)-(65), we obtain the Arrhenius law [24] k f,eq = (cid:114) πkT eq m d f exp (cid:104) − (cid:15) f kT eq (cid:105) , k r,eq = (cid:114) πkT eq m d r exp (cid:104) − (cid:15) r kT eq (cid:105) (68)In equilibrium, the reaction rates satisfy the follow law of mass-action [24] k f,eq k r,eq = (cid:18) m m m m (cid:19) exp (cid:20) − QkT eq (cid:21) = n n n n (69)Substituting equilibrium reaction rates (68) into the law of mass action (69), we find a relation between the stericfactors [13] s f √ m m d = s r √ m m d . (70)15 emark 7. Using the law of mass action, we make a distinction between a thermal and a chemical equilibrium.At a thermal equilibrium, the contribution from mechanical collisions vanishes Q M = 0 , all the mixture componentshave the same temperature T = T = T = T = T eq and the probability density function is given by f α = n α (cid:18) πθ α (cid:19) d exp (cid:18) − m | C α | kT eq (cid:19) . (71) In a chemical equilibrium, the contribution from the chemical collisions vanishes Q α,R = 0 , all the mixture compo-nents have the same temperature and the probability density function is given by the above expression. In addition,the number densities satisfy the law of mass action given in (69) . Note that a mixture in chemical equilibrium isalso in thermal equilibrium but vice-versa does not necessarily holds.
4. Numerical Experiments
Through numerical experiments we study how the Grad’s-14 moment system (given in (42)) relaxes to theequilibrium state. In particular, we study the effect of the following parameters on the relaxation behaviour (i)initial conditions, (ii) steric factors, (iii) activation energies and (iv) the heat of the reaction. To solve the momentsystem, we use the
NDSolve routine from mathematica with the default parameters. For the simplicity of exposition,we first non-dimensionalize the moment system. We use the same technique as that proposed in [10].
We non-dimensionalize the moment system (42) using the following scalesˆ x i = x i L , ˆ t = v o L t, ˆ v i = v i v o , ˆ n α = n α n oα , ˆ T α = T α T o , ˆ u ( α ) i = u ( α ) i (cid:112) θ oα ˆ σ ( α ) ij = σ ( α ) ij ρ oα θ oα , ˆ q ( α ) i = q ( α ) i ρ oα ( θ oα ) , ˆ∆ α = ∆ α ρ oα ( θ oα ) , ˆ (cid:15) f = (cid:15) f kT o , ˆ (cid:15) r = (cid:15) r kT o (72)where all the quantities with an ’o’ represent the reference scale. In the following discussion, a hat over any quantitywill denote its non-dimensionalized counterpart. Defining the total relative number density as n o = (cid:80) α =1 n oα , wedefine the mole-fraction of gas- α as F α = n oα /n o . Using the mole fraction, an averaged mass of the mixture and thevelocity scale can be defined as m = (cid:80) α =1 F α m α and v o = (cid:113) kT o m , respectively. We define the Knudsen number ofthe mixture as Kn = 516 √ πn o L ˜Ω (73)where L is the macroscopic length scale and the factor ˜Ω can be looked upon as an average cross-section for themixture. For a hard-sphere interaction potential, an explicit expression for ˜Ω can be given as ˜Ω = (cid:80) α =1 ( F α d α ).Defining the Knudsen number as given in (73) introduces a factor Ω αβ = ˜Ω d αβ in the dimensionless production terms.This factor can be looked upon as a scaled cross-section for every component in the mixture. For all the test cases, we consider the physical parameters given in Table 1. The parameters do not necessarilycorrespond to a realistic physical system and have been chosen only for demonstration purposes. Note that all themolecules are of the same size i.e., Ω αβ = 1, and that all the components of the mixture are equally diluted i.e., F α = 0 .
25. 16arameters m m m m Kn Ω αβ F α values 11.7 3.6 8 7.3 0.1 1 0.25 Table 1: Physical parameters used in all the experiments.
For all the coming experiments, we determine the equilibrium state as follows. We represent the initial numberdensity by ˆ n α, , the initial temperature by ˆ T α, , and we start from a fluid at rest. The initial temperature of themixture is denoted by T . In chemical equilibrium, all the gases have the same temperature T eq . Furthermore, thedensity of gas- α at equilibrium is given as ˆ n α, + x . The value of x and that of T eq follows by solving the equation[18] (cid:18) m m m m (cid:19) ( n + x )( n + x ) exp (cid:18) − / (cid:20) kT Q + 2 x n (cid:21)(cid:19) = ( n − x )( n − x ) , (74)where x = nk ( T ( eq ) − T ) /Q . Recall that n is the number density of the mixture and that it remains constantover time due to total mass conservation. Remark 8.
For all the test cases, we initialize the flow velocity ˆ u ( α ) i , the stress tensor ˆ σ ( α ) ij , the heat-flux ˆ q ( α ) i and ˆ∆ α with zero. For such an initialization, using the moment system (42) , it can be shown that ˆ u ( α ) i , ˆ q ( α ) i and ˆ σ ( α ) ij remain zero for all time instances therefore, we do not discuss them further. Note that ˆ∆ α is influenced by massand energy production and is therefore, perturbed by the gas mixture being outside of equilibrium. The followingexperiments study the perturbations in ˆ∆ α in detail.4.4. Influence of initial conditions4.4.1. Case-1: Exothermic Reaction The initial conditions and the physical parameters corresponding to this test case are given in Table 2 andTable 3, respectively. Furthermore, the equilibrium state computed using (74) is given in Table 2.The time variation of ˆ n α , ˆ T α and ˆ T is shown in 2(a), 2(b) and 2(c), respectively. Clearly, the Grad-14 momentsystem relaxes to the correct equilibrium state. Furthermore, due to an exothermic nature of the reaction, themixture temperature increases—see 2(c). In 2(b) and 2(d), the zoomed in plot shows the initial layer. The initiallayer is a result of the initial temperature difference between the different gases, which results in a deviation froman equilibrium. 2(d) shows the time evolution of ˆ∆ α . Even though we start with a zero value of ˆ∆ α , it is perturbedto a non-zero value. It is noteworthy that all the different ˆ∆ α ’s relax to zero at the same time scale as the numberdensities and the mixture temperatures.One can identify the two time-scales from the temperature profile shown in 2(b). At one time scale, which ismuch shorter than the other, the temperature of all the mixture components equalize to some value, say ˜ T . Atthe other time scale, the value ˜ T increases and approaches the equilibrium temperature T eq . These two distincttime-scales can also be identified in ˆ∆ α ’s profile shown in 2(d) where in the first time-scale, ˆ∆ α is pushed towardsits equilibrium value of zero and in the second time-scale, ˆ∆ α first increases and then approaches its equilibriumvalue of zero. Following is a possible explanation for the observance of these two time-scales. Initially, since thetemperature of any mixture component is not sufficiently high, mechanical interactions dominate. These interactionstry to push for a mechanical equilibrium where all the temperatures are equal and ˆ∆ α is zero. This results in thefirst time-scale. The second time-scale is triggered by chemically reactive interactions that increase the mixturetemperature and push ˆ∆ α away from its equilibrium value of zero. The system eventually reaches equilibriumwhere all the components have the same temperature, the number densities satisfy the law of mass action (69) andˆ∆ α = 0. 17nitial conditionsmoments Component-1 Component-2 Component-3 Component-4number density(ˆ n α ) 2 5 3 4Temperature( ˆ T α ) 1.5 2 3 1Equilibrium valuesmoments Component-1 Component-2 Component-3 Component-4number density(ˆ n α ) 4.8 7.8 0.19 1.19Temperature( ˆ T α ) 7.2 7.2 7.2 7.2 Table 2: Initial data and the corresponding equilibrium state for the exothermic reaction test case.
Parameters Forward Reaction Reverse Reaction s f , s r (cid:15) f , (cid:15) r
50 10
Table 3: Activation energies and steric factors for the exothermic reaction test case.(a) number density ˆ n α (b) Temperature ˆ T α (c) Mixture temperature ˆ T mix (d) The scalar moment ˆ∆ α Figure 2: Relaxation to equilibrium of the different moments.
We consider initial conditions that are in thermal equilibrium but not in chemical equilibrium. The initialconditions and the corresponding equilibrium states are given in Table 4. Since all the mixture components havethe same temperature, the mixture is in thermal equilibrium. However, the mixture is not in chemical equilibriumsince the law of mass action (69) is not satisfied. We refer to remark 7 for a distinction between a thermal and achemical equilibrium.Fig 3(a) shows the time variation of the different temperatures ˆ T α . Even though we start from the same18emperatures, due to a chemical non-equilibrium, the values of the different temperatures initially deviate fromeach other as the time progresses. Similar to the previous test case, ˆ∆ α is perturbed from its initial value and weobserve an initial layer which, due to a different initial data, is not as strong as the previous test case.Initial conditionsvariable Component-1 Component-2 Component-3 Component-4number density(ˆ n α ) 2 5 3 4Temperature( ˆ T α ) 2 2 2 2Equilibrium valuesvariable Component-1 Component-2 Component-3 Component-4number density(ˆ n α ) 4.7 7.7 0.21 1.21Temperature( ˆ T α ) 7.3 7.3 7.3 7.3 Table 4: Initial conditions in thermal equilibrium and the corresponding equilibrium state.(a) Temperature ˆ T ( α ) (b) The scalar moment ˆ∆ ( α ) Figure 3: Relaxation to equilibrium of the temperature and the scalar fourteenth moment ˆ∆ α . Table 2 shows the initial data and the corresponding equilibrium state. The activation energies and the stericfactors are given in Table 3 and Table 5, respectively. Note that we only specify s r , the value of s f follows fromthe relation in (70).As is clear from (74), the final equilibrium state is independent of the steric factors. However, the time-scaleat which the system relaxes to the equilibrium decreases as the steric factor increases—see Figure 4 that showsthe variation of ˆ T for different steric factors. Following is a possible explanation for this behaviour. From (32)and (33) we conclude that increasing s f , which also increases s r due to (70), the chemical collision cross-sectionincreases, resulting in an increased number of chemical collisions. Increase in the number of collisions results in afaster relaxation of the moment system.parameters Forward Reaction Reverse Reaction s f , s r s f , s r s f , s r s f , s r (cid:15) f , (cid:15) r
50 10
Table 5: Activation energies and the different steric factors. igure 4: Mixture temperature for different steric factors. Q ) The initial conditions and the equilibrium states are shown in Table 2 and Table 6, respectively. Note that,unlike the steric factor, the equilibrium state changes with the heat of the reaction. We fix the steric factor to s f = 0 . T for different values of the heat of the reaction. Sinceall the cases are exothermic, the equilibrium mixture temperature increases as the forward activation energy ( (cid:15) f )increases. Furthermore, increasing the heat of the reaction results in a faster increase in the mixture’s temperature,leading to a faster relaxation of the moment system to the equilibrium.Cases (cid:15) f (cid:15) r Q T ( eq ) Table 6: Activation energies and the equilibrium values for temperatures.Figure 5: Mixture temperature for different reaction heats.
5. Conclusion
We derived the Grad’s–14 moment equations for a chemically reactive quaternary gaseous mixture. The maincontribution was an algorithm to compute the moments of the Boltzmann’s collision operator. The algorithm reliedon a novel mathematical model that describes the microscopic details of a chemically reactive collision. To developthe model, we assumed that chemical reactions affect the tangential and the normal pre-collisional relative velocity20qually, leading to a symmetric collision. Using the collision model, we derived relations between the pre and thepost collisional velocities. These relations helped us extend the framework for the computation of the moments ofa single-gas collision operator to a reactive gas mixture. For first-order chemical kinetics, we derived reaction ratesfor chemical reactions outside of equilibrium and extended the Arhenius law. The reaction rates were found tohave an explicit dependence on the scalar fourteenth moment, highlighting the importance of considering a fourteenmoment system rather than a thirteen one. Our numerical experiments showcased that the 14-moment system cancorrectly reproduce the equilibrium states under different initializations and physical parameters.
References [1] Harold Grad. On the kinetic theory of rarefied gases.
Communications on Pure and Applied Mathematics ,2(4):331–407, 1949.[2] Neeraj Sarna, Jan Giesselmann, and Manuel Torrilhon. Convergence analysis of Grad’s hermite expansion forlinear kinetic equations.
SIAM Journal on Numerical Analysis , 58(2):1164–1194, 2020.[3] Christian Schmeiser and Alexander Zwirchmayr. Convergence of moment methods for linear kinetic equations.
SIAM Journal on Numerical Analysis , 36(1):74–88, 1998.[4] Manuel Torrilhon. Convergence study of moment approximations for boundary value problems of theBoltzmann-BGK equation.
Communications in Computational Physics , 18:529–557, 2015.[5] Zhenning Cai, Yuwei Fan, and Ruo Li. Globally hyperbolic regularization of Grad’s moment system.
Commu-nications on Pure and Applied Mathematics , 67(3):464–518, 2014.[6] Manuel Torrilhon. Modeling nonequilibrium gas flow based on moment equations.
Annual Review of FluidMechanics , 48(1):429–458, 2016.[7] Henning Struchtrup.
Macroscopic Transport Equations for Rarefied Gas Flows . Springer Ltd, 2010.[8] S. Chapman and T. G. Cowling.
The Mathematical Theory of Non-uniform Gases , volume 45. CambridgeUniversity Press, 1941.[9] Vinay Kumar Gupta and Manuel Torrilhon. Automated Boltzmann collision integrals for moment equations.In
AIP Conference Proceedings-American Institute of Physics , volume 1501, page 67, 2012.[10] Vinay Kumar Gupta and Manuel Torrilhon. Higher order moment equations for rarefied gas mixtures.
Pro-ceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 471:20140754, 2015.[11] Iain D Boyd and Thomas E Schwartzentruber.
Nonequilibrium Gas Dynamics and Molecular Simulation ,volume 42. Cambridge University Press, 2017.[12] Marzia Bisi, Maria Groppi, and Giampiero Spiga. Grad’s distribution functions in the kinetic equations for achemical reaction.
Continuum Mechanics and Thermodynamics , 14(2):207–222, 2002.[13] Gilberto Kremer.
An Introduction to the Boltzmann Equation and Transport Processes in Gases . Springer Ltd,2010.[14] G. A Bird.
Molecular gas dynamics and the direct simulation of gas flows . Oxford : Clarendon Press, 1995.[15] Sergey F Gimelshein and Ingrid J Wysong. Birds total collision energy model: 4 decades and going strong.
Physics of Fluids , 31(7):076101, 2019. 2116] Charles R Lilley and Michael N Macrossan. A macroscopic chemistry method for the direct simulation of gasflows.
Physics of Fluids , 16(6):2054–2066, 2004.[17] VA Rykov. Kinetic equations for chemically reacting gas mixtures.
Fluid Dynamics , 7(4):648–655, 1972.[18] A. Rossani and G. Spiga. A note on the kinetic theory of chemically reacting gases.
Physica A StatisticalMechanics and its Applications , 272:563–573, 1999.[19] Adriano W. Silva, Giselle M. Alves, and Gilberto M. Kremer. Transport phenomena in a reactive quaternarygas mixture.
Physica A: Statistical Mechanics and its Applications , 374(2):533–548, 2007.[20] Elena V Kustova and Gilberto M Kremer. Chemical reaction rates and non-equilibrium pressure of reactinggas mixtures in the state-to-state approach.
Chemical Physics , 445:82–94, 2014.[21] Elene V Kustova and Georgii P Oblapenko. Mutual effect of vibrational relaxation and chemical reactions inviscous multitemperature flows.
Physical Review E , 93(3):033127, 2016.[22] John C. Light, John Ross, and Kurt E. Shuler.
Kinetic Processes in Gases and Plasmas . Elsevier, 1969.[23] R.D. Present. Note on the simple theory of bimolecular reactions.
Proc. Natl. Acad. Sci. USA 41 , (1):415–417,1955.[24] P.W Atkins.
Physical Chemistry , volume 45. Oxford University Press, Oxford, 1997.[25] Nikolai V. Brilliantov, Frank Spahn, Jan-Martin Hertzsch, and Thorsten P¨oschel. Model for collisions ingranular gases.
Physical Review E , 53(5):5382–5392, 1996.[26] Gilberto Kremer and Wilson Marques Jr. Fourteen moment theory for granular gases.
Kinetic and RelatedModels , 4(1):317–331, 2011.[27] Dino Risso and Patricio Cordero. Dynamics of rarefied granular gases.
Phys. Rev. E , 65:021304, 2002.22 . Supplementary material: Omega Integrals Ω (3 , (cid:15) r = 4 m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17) θ M m (75)Ω (5 , (cid:15) r = 16 m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17) θ M m (76)Ω (5 , (cid:15) r = 8 m (cid:16) Q Γ (cid:16) , (cid:15) r θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17)(cid:17) θ M m (77)Ω (7 , (cid:15) r = 64 m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17) θ M m (78)Ω (9 , (cid:15) r = 256 m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17) θ M m (79)Ω (7 , (cid:15) r = 32 m (cid:16) Q Γ (cid:16) , (cid:15) r θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17)(cid:17) θ M m (80)Ω (9 , (cid:15) r = 128 m (cid:16) Q Γ (cid:16) , (cid:15) r θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17)(cid:17) θ M m (81)Ω (11 , (cid:15) r = 1024 m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17) θ M m (82)Ω (11 , (cid:15) r = 512 m (cid:16) Q Γ (cid:16) , (cid:15) r θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) r θ m (cid:17) + Q Γ (cid:16) , (cid:15) r θ m (cid:17)(cid:17)(cid:17) θ M m (83)23 .2. Omega integrals for the forward reaction Ω (3 , (cid:15) f = 4 m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17) θ M m (84)Ω (5 , (cid:15) f = 16 m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17) θ M m (85)Ω (5 , (cid:15) f = 8 m (cid:16) Q Γ (cid:16) , (cid:15) f θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17)(cid:17) θ M m (86)Ω (7 , (cid:15) f = 64 m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17) θ M m (87)Ω (9 , (cid:15) f = 256 m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17) θ M m (88)Ω (7 , (cid:15) f = 32 m (cid:16) Q Γ (cid:16) , (cid:15) f θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17)(cid:17) θ M m (89)Ω (9 , (cid:15) f = 128 m (cid:16) Q Γ (cid:16) , (cid:15) f θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17)(cid:17) θ M m (90)Ω (11 , (cid:15) f = 1024 m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17) θ M m (91)Ω (11 , (cid:15) f = 512 m (cid:16) Q Γ (cid:16) , (cid:15) f θ m (cid:17) + 4 θ m (cid:16) θ m Γ (cid:16) , (cid:15) f θ m (cid:17) − Q Γ (cid:16) , (cid:15) f θ m (cid:17)(cid:17)(cid:17) θ M m (92)where M = m + m = m + m .
7. Supplementary material: Chemical production terms • Mass Balance P α =1 R, = d r ˜ n µ exp (cid:20) − (cid:15) r θ m (cid:21) − d f ˜ n µ exp (cid:20) − (cid:15) f θ m (cid:21) + ∆ r (1)0 (93) P α =2 R, = d r ˜ n µ exp (cid:20) − (cid:15) r θ m (cid:21) − d f ˜ n µ exp (cid:20) − (cid:15) f θ m (cid:21) + ∆ r (2)0 (94) P α =3 R, = d f ˜ n µ exp (cid:20) − (cid:15) f θ m (cid:21) − d r ˜ n µ exp (cid:20) − (cid:15) r θ m (cid:21) + ∆ r (3)0 (95) P α =4 R, = d f ˜ n µ exp (cid:20) − (cid:15) f θ m (cid:21) − d r ˜ n µ exp (cid:20) − (cid:15) r θ m (cid:21) + ∆ r (4)0 (96)where ˜ n = 4 √ π √ θ M n n , ˜ n = 4 √ π √ θ M n n and ∆ r ( α )0 is the contribution from ∆ α . Explicit expressionsfor ∆ r ( α )0 have been given in subsection 8.2 and subsection 8.1.24 .2. Momentum Balance The generic form of the momentum production can be given as P αR,i = (cid:88) γ =1 (cid:16) τ αγ u ( γ ) i + ω αγ q ( γ ) i (cid:17) (97)The coefficients τ αγ and ω αγ can now be given as τ (1)1 = 16 β f (cid:112) θ µ (cid:0) θ m (cid:0) θ + 8 θ (5 θ − θ ) + 3 θ (5 θ − θ ) (cid:1) + 2 θ θ m (cid:15) f (4 θ − θ + 5 θ ) + θ (cid:15) f (cid:1) τ (1)2 = − β f θ (cid:112) θ µ (cid:0) θ (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + θ m (5 θ − θ )( θ m + (cid:15) f ) (cid:1) τ (1)3 = 16 β r (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(6 θ − θ + 5 θ ) + 56 θ θ − θ θ + 42 θ θ − θ (cid:1) + θ m (cid:15) r ( µ (2 θ )(8 θ − θ + 5 θ ) − θ (4 θ − θ + 5 θ )) + θ (cid:15) r ( θ ( µ −
1) + θ µ ) (cid:1) τ (1)4 = − β r (cid:112) θ µ (cid:0) µ (2 θ ) (cid:0) θ (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) + θ m (5 θ − θ )(4 θ m + (cid:15) r ) (cid:1) − θ θ (cid:0) θ m + 3 θ m (cid:15) r + (cid:15) r (cid:1) + θ θ ( − m )(5 θ − θ )( θ m + (cid:15) r ) (cid:1) ω (1)1 = − β f (cid:112) θ µ (cid:0) θ m (cid:0) θ + 40 θ ( θ − θ ) + 15 θ ( θ − θ ) (cid:1) + 2 θ θ m (cid:15) f (4 θ − θ + 5 θ ) + θ (cid:15) f (cid:1) ω (1)2 = 115 β f θ (cid:112) θ µ (cid:0) θ (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + 5 θ m ( θ − θ )( θ m + (cid:15) f ) (cid:1) ω (1)3 = − β r (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(6 θ + 5( θ − θ )) + 40 θ θ − θ θ + 30 θ θ − θ (cid:1) + θ m (cid:15) r ( µ (2 θ )(8 θ + 5( θ − θ )) − θ (4 θ − θ + 5 θ )) + θ (cid:15) r ( θ ( µ −
1) + θ µ ) (cid:1) ω (1)4 = 115 β r (cid:112) θ µ (cid:0) µ (2 θ ) (cid:0) θ (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) + 5 θ m ( θ − θ )(4 θ m + (cid:15) r ) (cid:1) − θ θ (cid:0) θ m + 3 θ m (cid:15) r + (cid:15) r (cid:1) − θ θ m ( θ − θ )( θ m + (cid:15) r ) (cid:1) .2.2. Component-2 τ (2)1 = − β f (cid:112) θ θ µ (cid:0) θ m (4 θ − θ + 5 θ ) + θ m (cid:15) f (3 θ − θ + 5 θ ) + θ (cid:15) f (cid:1) τ (2)2 = 16 β f (cid:112) θ µ (cid:0) θ θ m (5 θ − θ )+ θ (cid:0) θ m +8 θ m (cid:15) f + (cid:15) f (cid:1) +2 θ θ m (5 θ − θ )(4 θ m + (cid:15) f ) (cid:1) τ (2)3 = − β r (cid:112) θ µ (cid:0) θ m (4 µ (2 θ )(6 θ − θ + 5 θ ) + θ ( − θ + 14 θ − θ ))+ θ m (cid:15) r ( µ (2 θ )(8 θ − θ + 5 θ ) + θ ( − θ + 14 θ − θ )) + θ (cid:15) r ( θ µ + θ ( µ − (cid:1) τ (2)4 = 16 β r (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(5 θ − θ + 6 θ ) + 42 θ θ − θ θ + 56 θ θ − θ (cid:1) + θ m (cid:15) r ( µ (2 θ )(5 θ − θ + 8 θ ) − θ (5 θ − θ + 4 θ )) + θ (cid:15) r ( θ µ + θ ( µ − (cid:1) ω (2)1 = 115 β f (cid:112) θ θ µ (cid:0) θ m (4 θ + 5( θ − θ )) + θ m (cid:15) f (3 θ + 5( θ − θ )) + θ (cid:15) f (cid:1) ω (2)2 = − β f (cid:112) θ µ (cid:0) θ θ m ( θ − θ )+ θ (cid:0) θ m +8 θ m (cid:15) f + (cid:15) f (cid:1) +10 θ θ m ( θ − θ )(4 θ m + (cid:15) f ) (cid:1) ω (2)3 = 115 β r (cid:112) θ µ (cid:0) θ m (4 µ (2 θ )(6 θ + 5( θ − θ )) + θ ( − θ + 10 θ − θ ))+ θ m (cid:15) r ( µ (2 θ )(8 θ + 5( θ − θ )) + θ ( − θ + 10 θ − θ )) + θ (cid:15) r ( θ µ + θ ( µ − (cid:1) ω (2)4 = − β r (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(5 θ − θ + 6 θ ) + 30 θ θ − θ θ + 40 θ θ − θ (cid:1) + θ m (cid:15) r ( µ (2 θ )(5 θ − θ + 8 θ ) − θ (5 θ − θ + 4 θ )) + θ (cid:15) r ( θ µ + θ ( µ − (cid:1) τ (3)1 = 16 β f (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(6 θ − θ + 5 θ ) + 56 θ θ − θ θ + 42 θ θ − θ (cid:1) + θ m (cid:15) f ( µ (2 θ )(8 θ − θ + 5 θ ) − θ (4 θ − θ + 5 θ )) + θ (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:1) τ (3)2 = − β f (cid:112) θ µ (cid:0) µ (2 θ ) (cid:0) θ (cid:0) θ m + 8 θ m (cid:15) f + (cid:15) f (cid:1) + θ m (5 θ − θ )(4 θ m + (cid:15) f ) (cid:1) − θ θ (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + θ θ ( − m )(5 θ − θ )( θ m + (cid:15) f ) (cid:1) τ (3)3 = 16 β r (cid:112) θ µ (cid:0) θ m (cid:0) θ + 8 θ (5 θ − θ ) + 3 θ (5 θ − θ ) (cid:1) + 2 θ θ m (cid:15) r (4 θ − θ + 5 θ ) + θ (cid:15) r (cid:1) (3)4 = − β r θ (cid:112) θ µ (cid:0) θ (cid:0) θ m + 3 θ m (cid:15) r + (cid:15) r (cid:1) + θ m (5 θ − θ )( θ m + (cid:15) r ) (cid:1) ω (3)1 = − β f (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(6 θ + 5( θ − θ )) + 40 θ θ − θ θ + 30 θ θ − θ (cid:1) + θ m (cid:15) f ( µ (2 θ )(8 θ + 5( θ − θ )) − θ (4 θ − θ + 5 θ )) + θ (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:1) ω (3)2 = 115 β f (cid:112) θ µ (cid:0) µ (2 θ ) (cid:0) θ (cid:0) θ m + 8 θ m (cid:15) f + (cid:15) f (cid:1) + 5 θ m ( θ − θ )(4 θ m + (cid:15) f ) (cid:1) − θ θ (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) − θ θ m ( θ − θ )( θ m + (cid:15) f ) (cid:1) ω (3)3 = − β r (cid:112) θ µ (cid:0) θ m (cid:0) θ + 40 θ ( θ − θ ) + 15 θ ( θ − θ ) (cid:1) + 2 θ θ m (cid:15) r (4 θ − θ + 5 θ ) + θ (cid:15) r (cid:1) ω (3)4 = 115 β r θ (cid:112) θ µ (cid:0) θ (cid:0) θ m + 3 θ m (cid:15) r + (cid:15) r (cid:1) + 5 θ m ( θ − θ )( θ m + (cid:15) r ) (cid:1) τ (4)1 = − β f (cid:112) θ µ (cid:0) θ m (4 µ (2 θ )(6 θ − θ + 5 θ ) + θ ( − θ + 14 θ − θ ))+ θ m (cid:15) f ( µ (2 θ )(8 θ − θ + 5 θ ) + θ ( − θ + 14 θ − θ )) + θ (cid:15) f ( θ µ + θ ( µ − (cid:1) τ (4)2 = 16 β f (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(5 θ − θ + 6 θ ) + 42 θ θ − θ θ + 56 θ θ − θ (cid:1) + θ m (cid:15) f ( µ (2 θ )(5 θ − θ + 8 θ ) − θ (5 θ − θ + 4 θ )) + θ (cid:15) f ( θ µ + θ ( µ − (cid:1) τ (4)3 = − β r (cid:112) θ θ µ (cid:0) θ m (4 θ − θ + 5 θ ) + θ m (cid:15) r (3 θ − θ + 5 θ ) + θ (cid:15) r (cid:1) τ (4)4 = 16 β r (cid:112) θ µ (cid:0) θ θ m (5 θ − θ ) + θ (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) + 2 θ θ m (5 θ − θ )(4 θ m + (cid:15) r ) (cid:1) ω (4)1 = 115 β f (cid:112) θ µ (cid:0) θ m (4 µ (2 θ )(6 θ + 5( θ − θ )) + θ ( − θ + 10 θ − θ ))+ θ m (cid:15) f ( µ (2 θ )(8 θ + 5( θ − θ )) + θ ( − θ + 10 θ − θ )) + θ (cid:15) f ( θ µ + θ ( µ − (cid:1) ω (4)2 = − β f (cid:112) θ µ (cid:0) θ m (cid:0) − θ + 4 µ (2 θ )(5 θ − θ + 6 θ ) + 30 θ θ − θ θ + 40 θ θ − θ (cid:1) + θ m (cid:15) f ( µ (2 θ )(5 θ − θ + 8 θ ) − θ (5 θ − θ + 4 θ )) + θ (cid:15) f ( θ µ + θ ( µ − (cid:1) ω (4)3 = 115 β r (cid:112) θ θ µ (cid:0) θ m (4 θ + 5( θ − θ )) + θ m (cid:15) r (3 θ + 5( θ − θ )) + θ (cid:15) r (cid:1) (4)4 = − β r (cid:112) θ µ (cid:0) θ θ m ( θ − θ )+ θ (cid:0) θ m +8 θ m (cid:15) r + (cid:15) r (cid:1) +10 θ θ m ( θ − θ )(4 θ m + (cid:15) r ) (cid:1) β f = (cid:32) √ πd M n n s f θ / m (cid:33) exp (cid:20) − (cid:15) f θ m (cid:21) , β r = (cid:32) √ πd M n n s r θ / m (cid:33) exp (cid:20) − (cid:15) r θ m (cid:21) (98) The following are the production terms corresponding to energy balance for all the components present in themixture. P αR, = 18 β r θ µ (cid:16) θ (2 θ ) e (cid:15)r θ m ( θ m Ω (5 , (cid:15) r − (cid:15) r Ω (3 , (cid:15) r )+ 8 (cid:0) θ µ (2 θ )(4 θ m + (cid:15) r ) − θ µ (2 θ )(4 θ m + (cid:15) r ) + θ ( θ m (4 θ + 3 θ ) + θ (cid:15) r ) (cid:1)(cid:17) − β f θ θ µ ( θ m (4 θ + 3 θ ) + θ (cid:15) f ) + ∆ r (1)1 P αR, = β r θ µ (cid:16) θ µ (2 θ ) e (cid:15)r θ m ( θ m Ω (5 , (cid:15) r − (cid:15) r Ω (3 , (cid:15) r )+ 8 µ (cid:0) θ µ (2 θ )(4 θ m + (cid:15) r ) − θ µ (2 θ )(4 θ m + (cid:15) r ) + θ ( θ m (3 θ + 4 θ ) + θ (cid:15) r ) (cid:1)(cid:17) − β f θ θ µ ( θ m (3 θ + 4 θ ) + θ (cid:15) f ) + ∆ r (2)1 P αR, = 18 β f θ µ (cid:16) θ (2 θ ) e (cid:15)f θ m ( θ m Ω (5 , (cid:15) f − (cid:15) f Ω (3 , (cid:15) f )+ 8 (cid:0) θ µ (2 θ )(4 θ m + (cid:15) f ) − θ µ (2 θ )(4 θ m + (cid:15) f ) + θ ( θ m (4 θ + 3 θ ) + θ (cid:15) f ) (cid:1)(cid:17) − β r θ θ µ ( θ m (4 θ + 3 θ ) + θ (cid:15) r ) + ∆ r (3)1 P αR, = β f θ µ (cid:16) θ µ (2 θ ) e (cid:15)f θ m ( θ m Ω (5 , (cid:15) f − (cid:15) f Ω (3 , (cid:15) f )+ 8 µ (cid:0) θ µ (2 θ )(4 θ m + (cid:15) f ) − θ µ (2 θ )(4 θ m + (cid:15) f ) + θ ( θ m (3 θ + 4 θ ) + θ (cid:15) f ) (cid:1)(cid:17) − β r θ θ µ ( θ m (3 θ + 4 θ ) + θ (cid:15) r ) + ∆ r (4)1 .4. Stress Production The production term for the stress tensor can be given in a generic form as P αR,ij = (cid:88) γ =1 (cid:16) ξ ( α ) γ σ ( α ) ij (cid:17) (99)The coefficients ξ ( α ) γ can now be given as ξ (1)1 = − β f θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 15 θ (cid:1) + 2 θ θ m (cid:15) f (4 θ + 5 θ ) + θ (cid:15) f (cid:1) ξ (1)2 = 115 β f θ θ µ (cid:0) θ m + 2 θ m (cid:15) f − (cid:15) f (cid:1) ξ (1)3 = 115 β r θ µ (cid:0) θ m (cid:0) θ + 48 θ µ (2 θ ) − µ (2 θ )(6 θ + 5 θ ) + 40 θ θ + 15 θ (cid:1) + 2 θ m (cid:15) r (cid:0) θ µ (2 θ ) − µ (2 θ )(8 θ + 5 θ ) + θ (4 θ + 5 θ ) (cid:1) + (cid:15) r (cid:0) θ + 2 θ µ (2 θ ) − θ µ (2 θ ) (cid:1)(cid:1) ξ (1)4 = 115 β r θ µ (cid:0) θ (cid:0) − θ m − θ m (cid:15) r + (cid:15) r (cid:1) + 2 θ µ (2 θ ) (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) − θ µ (2 θ ) (cid:0) θ m + 3 θ m (cid:15) r + (cid:15) r (cid:1)(cid:1) ξ (2)1 = 115 β f θ θ µ (cid:0) θ m + 2 θ m (cid:15) f − (cid:15) f (cid:1) ξ (2)2 = − β f θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 24 θ (cid:1) + 2 θ θ m (cid:15) f (5 θ + 4 θ ) + θ (cid:15) f (cid:1) ξ (2)3 = 115 β r θ µ (cid:0) θ µ (2 θ ) (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) − θ µ (2 θ ) (cid:0) θ m + 3 θ m (cid:15) r + (cid:15) r (cid:1) + θ (cid:0) − θ m − θ m (cid:15) r + (cid:15) r (cid:1)(cid:1) ξ (2)4 = 115 β r θ µ (cid:0) θ m (cid:0) θ + 48 θ µ (2 θ ) − µ (2 θ )(5 θ + 6 θ ) + 40 θ θ + 24 θ (cid:1) + 2 θ m (cid:15) r (cid:0) θ µ (2 θ ) − µ (2 θ )(5 θ + 8 θ ) + θ (5 θ + 4 θ ) (cid:1) + (cid:15) r (cid:0) θ µ (2 θ ) − θ µ (2 θ ) + θ (cid:1)(cid:1) .4.3. Component-3 ξ (3)1 = 115 β f θ µ (cid:0) θ m (cid:0) θ + 48 θ µ (2 θ ) − µ (2 θ )(6 θ + 5 θ ) + 40 θ θ + 15 θ (cid:1) + 2 θ m (cid:15) f (cid:0) θ µ (2 θ ) − µ (2 θ )(8 θ + 5 θ ) + θ (4 θ + 5 θ ) (cid:1) + (cid:15) f (cid:0) θ + 2 θ µ (2 θ ) − θ µ (2 θ ) (cid:1)(cid:1) ξ (3)2 = 115 β f θ µ (cid:0) θ (cid:0) − θ m − θ m (cid:15) f + (cid:15) f (cid:1) + 2 θ µ (2 θ ) (cid:0) θ m + 8 θ m (cid:15) f + (cid:15) f (cid:1) − θ µ (2 θ ) (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1)(cid:1) ξ (3)3 = − β r θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 15 θ (cid:1) + 2 θ θ m (cid:15) r (4 θ + 5 θ ) + θ (cid:15) r (cid:1) ξ (3)4 = 115 β r θ θ µ (cid:0) θ m + 2 θ m (cid:15) r − (cid:15) r (cid:1) ξ (4)1 = 115 β f θ µ (cid:0) θ µ (2 θ ) (cid:0) θ m + 8 θ m (cid:15) f + (cid:15) f (cid:1) − θ µ (2 θ ) (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + θ (cid:0) − θ m − θ m (cid:15) f + (cid:15) f (cid:1)(cid:1) ξ (4)2 = 115 β f θ µ (cid:0) θ m (cid:0) θ + 48 θ µ (2 θ ) − µ (2 θ )(5 θ + 6 θ ) + 40 θ θ + 24 θ (cid:1) + 2 θ m (cid:15) f (cid:0) θ µ (2 θ ) − µ (2 θ )(5 θ + 8 θ ) + θ (5 θ + 4 θ ) (cid:1) + (cid:15) f (cid:0) θ µ (2 θ ) − θ µ (2 θ ) + θ (cid:1)(cid:1) ξ (4)3 = 115 β r θ θ µ (cid:0) θ m + 2 θ m (cid:15) r − (cid:15) r (cid:1) ξ (4)4 = − β r θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 24 θ (cid:1) + 2 θ θ m (cid:15) r (5 θ + 4 θ ) + θ (cid:15) r (cid:1) The production term for the heat flux can be given in a generic form as P αR, ,i = N (cid:88) γ =1 (cid:16) κ ( α ) γ u ( α ) i + ς ( α ) γ q ( α ) i (cid:17) (100)The coefficients κ ( α ) γ and ς ( α ) γ can now be given as 30 .5.1. Component-1 κ (1)1 =124 β (2) f θ / (cid:0) θ m (cid:0) θ + 168 θ (3 θ − θ ) + 140 θ θ (3 θ − θ ) + 105 θ ( θ − θ ) (cid:1) + θ θ m (cid:15) f (cid:0) θ + 56 θ (3 θ − θ ) + 35 θ (3 θ − θ ) (cid:1) + θ θ m (cid:15) f (12 θ − θ + 21 θ ) + θ (cid:15) f (cid:1) κ (1)2 = − β (2) f θ (cid:112) θ (cid:0) θ θ m (5 θ − θ ) (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + 5 θ θ m (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + θ (cid:0) θ m (43 θ − θ ) + θ m (cid:15) f (39 θ − θ ) + θ θ m (cid:15) f + θ (cid:15) f (cid:1)(cid:1) κ (1)3 = β (2) r √ θ (cid:18) θ m e (cid:15)r θ m (4Ω (5 , (cid:15) r (3 θ θ m (14 θ − θ )+ 2 θ (cid:15) r (5 θ − θ ) + 2 θ µ (cid:15) r (14 θ − θ )) − (7 , (cid:15) r (2 θ µ ( θ m (14 θ − θ ) + θ (cid:15) r ) − θ (2 θ m (7 θ − θ ) + θ (cid:15) r ))+ θ θ m Ω (9 , (cid:15) r ( θ ( µ −
1) + θ µ ) + 24 θ (cid:15) r (5 θ − θ )Ω (3 , (cid:15) r )+ 96 (cid:0) θ m (cid:0) θ µ (cid:0) − θ + 252 θ θ − θ θ + 70 θ θ − θ (cid:1) + θ (cid:0) θ (7 θ − θ ) + 210 θ ( θ − θ ) + 56 θ θ (10 θ − theta ) − θ (cid:1) + 24 θ θ µ (cid:0) θ + 84 θ ( θ − θ ) + 35 θ ( θ − θ ) (cid:1) + 192 θ µ (8 θ − θ + 5 θ ) (cid:1) + θ m (cid:15) r (cid:0) θ µ (cid:0) − θ + 336 θ θ − θ θ + 70 θ θ − θ (cid:1) + 6 θ θ µ (cid:0) θ + 112 θ ( θ − θ ) + 35 θ ( θ − θ ) (cid:1) + θ (cid:0) − θ + 56 θ (2 θ − θ ) + 35 θ (4 θ − θ ) (cid:1) + 64 θ µ (9 θ − θ + 5 θ ) (cid:1) + θ m (cid:15) r (cid:0) θ ( − θ + 14 θ − θ ) + 12 θ θ µ (6 θ − θ + 7 θ ) + 8 θ µ (12 θ − θ + 5 θ ) − θ θ µ (18 θ − θ + 13 θ ) (cid:1) + θ (cid:15) r (cid:0) − θ + 6 θ θ µ − θ θ µ + 8 θ µ (cid:1)(cid:1) (cid:19) κ (1)4 = β (2) r √ θ (cid:18) θ m (2 θ ) e (cid:15)r θ m (4Ω (5 , (cid:15) r ( θ (14 θ − θ )(3 θ m + (cid:15) r ) + µ (cid:15) r (5 θ − θ )(2 θ )+5 θ θ (cid:15) r ) − (7 , (cid:15) r ( µ (2 θ )( θ m (5 θ − θ ) − θ (cid:15) r )+ θ ( θ m ( − θ +14 θ +5 θ )+ θ (cid:15) r ))+ θ θ m Ω (9 , (cid:15) r ( θ − µ (2 θ )) + 24 θ (cid:15) r (5 θ − θ )Ω (3 , (cid:15) r )+ 96 (cid:0) θ m (cid:0) − θ µ (2 θ ) (cid:0) θ − θ θ + 67 θ θ − θ θ + 30 θ (cid:1) + θ (cid:0) θ − θ θ + 43 θ θ − θ θ + 20 θ (cid:1) + 48 θ µ (2 θ ) ( − θ + 14 θ − θ ) + 104 θ θ µ (2 θ )(5 θ − θ + 6 θ ) (cid:1) + θ m (cid:15) r (cid:0) − θ µ (2 θ ) (cid:0) θ − θ θ + 95 θ θ − θ θ + 40 θ (cid:1) + θ (cid:0) θ − θ θ + 39 θ θ − θ θ + 15 θ (cid:1) + 16 θ µ (2 θ ) ( − θ + 14 θ − θ ) + 2 θ θ µ (2 θ )(95 θ − θ + 128 θ ) (cid:1) + θ m (cid:15) r (cid:0) − θ µ (2 θ ) (cid:0) θ + 14 θ ( θ − θ ) + 5 θ (cid:1) + θ (cid:0) θ + θ ( θ − θ ) + 5 θ (cid:1) + 2 θ µ (2 θ ) ( − θ + 14 θ − θ ) + 2 θ θ µ (2 theta )(15 θ − θ + 25 θ ) (cid:1) + θ (cid:15) r (cid:0) θ − θ µ (2 θ ) − θ µ (2 θ ) + 6 θ θ µ (2 θ ) (cid:1)(cid:1) (cid:19) ς (1)1 = 31 β (2) f θ / (cid:0) θ m (cid:0) θ (4 θ − θ ) + 30 θ (14 θ − θ ) + 8 θ θ (63 θ − θ ) + 105 θ (cid:1) + θ θ m (cid:15) f (cid:0) θ + θ (168 θ − θ ) + 5 θ (21 θ − θ ) (cid:1) + θ θ m (cid:15) f (12 θ − θ + 21 θ ) + θ (cid:15) f (cid:1) ς (1)2 =160 β (2) f θ (cid:112) θ (cid:0) θ θ m ( θ − θ ) (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + 5 θ θ m (cid:0) θ m + 3 θ m (cid:15) f + (cid:15) f (cid:1) + θ (cid:0) θ m (43 θ − θ ) + θ m (cid:15) f (39 θ − θ ) + θ θ m (cid:15) f + θ (cid:15) f (cid:1)(cid:1) ς (1)3 = − β (2) r √ θ (cid:18) θ m (2 θ ) e (cid:15)r θ m (20Ω (5 , (cid:15) r (3 θ θ m (2 θ − θ ) + µ (cid:15) r (2 θ )(2 θ − θ )+ 2 θ (cid:15) r ( θ − θ )) − (7 , (cid:15) r (5 θ m (2 θ ( θ − θ ) − µ (2 θ )( θ − θ )) + θ (cid:15) r ( µ (2 θ ) − θ ))+ θ θ m Ω (9 , (cid:15) r ( θ ( µ −
1) + θ µ ) + 120 θ (cid:15) r ( θ − θ )Ω (3 , (cid:15) r )+ 96 (cid:0) θ m (cid:0) θ µ (2 θ ) (cid:0) − θ + 180 θ θ − θ θ + 50 θ θ − θ (cid:1) + 12 θ µ (2 theta ) (cid:0) θ − θ θ + 84 θ θ − θ θ + 35 θ (cid:1) + θ (cid:0) θ (5 θ − θ ) + 30 θ (5 θ − θ ) + 8 θ θ (50 θ − θ ) − θ (cid:1) + 48 θ µ (2 θ ) (8 θ + 5( θ − θ )) (cid:1) + θ m (cid:15) r (cid:0) θ µ (2 θ ) (cid:0) − θ + 240 θ θ − θ θ + 50 θ θ − θ (cid:1) + 3 θ µ (2 θ ) (cid:0) θ + 16 θ (7 θ − θ ) + 5 θ (7 θ − θ ) (cid:1) + θ (cid:0) − θ + 8 θ (10 θ − θ ) + 5 θ (20 θ − θ ) (cid:1) + 16 θ µ (2 θ ) (9 θ + 5( θ − θ )) (cid:1) + θ m (cid:15) r (cid:0) θ ( − θ + 10 θ − θ ) + 6 θ µ (2 θ )(6 θ − θ + 7 θ )+ 2 θ µ (2 θ ) (12 θ + 5( θ − θ )) − θ θ µ (2 θ )(18 θ − θ + 13 θ ) (cid:1) + θ (cid:15) r (cid:0) − θ + 3 θ µ (2 θ ) + 2 θ µ (2 θ ) − θ θ µ (2 θ ) (cid:1)(cid:1) (cid:19) ς (1)4 = β (2) r √ θ (cid:18) θ m (2 θ ) e (cid:15)r θ m ( − (5 , (cid:15) r ( − θ ( θ − θ )(3 θ m + (cid:15) r ) + µ (cid:15) r ( θ − θ )(2 θ )+ θ θ (cid:15) r ) + 2Ω (7 , (cid:15) r (5 θ m ( µ ( θ − θ )(2 θ ) + θ ( − θ + 2 θ + θ )) + θ (cid:15) r ( θ − µ (2 θ )))+ θ θ m Ω (9 , (cid:15) r ( θ ( µ −
1) + θ µ ) − θ (cid:15) r ( θ − θ )Ω (3 , (cid:15) r )+ 96 (cid:0) θ m (cid:0) θ (cid:0) − θ + 40 θ θ − θ θ + 50 θ θ − θ (cid:1) +4 θ µ (2 θ ) lef t ( θ (67 θ − θ )+40 θ ( θ − θ )+30 θ (cid:1) +48 θ µ (2 θ ) (5 θ − θ +8 θ ) − θ θ µ (2 θ )(5 θ − θ +6 θ ) (cid:1) + θ m (cid:15) r (cid:0) θ (cid:0) − θ +30 θ θ − θ θ +50 θ θ − θ (cid:1) + 5 θ µ (2 θ ) (cid:0) θ (19 θ − θ ) + 14 θ ( θ − θ ) + 8 θ (cid:1) + 16 θ µ (2 θ ) (5 θ − θ + 9 θ ) − θ θ µ (2 θ )(95( θ − θ ) + 128 θ ) (cid:1) + θ m (cid:15) r (cid:0) − θ (cid:0) θ ( θ − θ ) + θ θ + 5 θ (cid:1) + θ µ (2 θ ) (cid:0) θ ( θ − θ ) + 14 θ θ + 5 θ )+ 2 θ µ (2 θ ) (5 θ − θ + 12 θ ) − θ θ µ (2 θ )(3 θ − θ + 5 θ ) (cid:1) + θ (cid:15) r (cid:0) − θ + 3 θ µ (2 θ ) + 2 θ µ (2 θ ) − θ θ µ (2 θ ) (cid:1)(cid:1) (cid:19) .5.2. Component-2 κ (2)1 = − µ β (2) f √ θ θ µ (cid:18) θ m (cid:0) θ − θ θ + 43 θ θ − θ θ + 20 θ (cid:1) + θ m (cid:15) f (cid:0) θ − θ θ + 39 θ θ − θ θ + 15 θ (cid:1) + θ m (cid:15) f (cid:0) θ + θ θ + θ (5 θ − θ ) (cid:1) + θ θ (cid:15) f (cid:19) θ κ (2)2 = µ β (2) f θ / µ (cid:18) θ θ m ( θ − θ ) + 7 θ θ m (3 θ − θ ) (cid:0) θ m + 8 θ m (cid:15) f + (cid:15) f (cid:1) + 35 θ θ θ m (3 θ − θ )(4 θ m + (cid:15) f ) + θ (cid:0) θ m + 72 θ m (cid:15) f + 12 θ m (cid:15) f + (cid:15) f (cid:1) (cid:19) κ (2)3 = − β (2) r √ θ µ µ (cid:18) θ µ m e (cid:15)r θ m ( − (5 , (cid:15) r (3 θ θ m (14 θ − θ )+ θ (cid:15) r (5 θ + 14 θ − θ ) + 2 θ µ (cid:15) r (5 θ − θ )) − (7 , (cid:15) r (2 θ µ ( θ m (14 θ − θ ) + θ (cid:15) r ) − θ ( θ m (5 θ + 14 θ − θ ) + θ (cid:15) r ))+ θ θ m Ω (9 , (cid:15) r ( θ µ + θ ( µ − θ (cid:15) r (14 θ − θ )Ω (3 , (cid:15) r )+ 96 µ (cid:0) θ m (cid:0) θ (cid:0) − θ + 70 θ θ − θ θ + 56 θ θ − θ (cid:1) + 8 θ θ µ (cid:0) θ + θ (67 θ − θ ) + 8 θ (5 θ − θ ) (cid:1) + 192 θ µ (8 θ − θ + 5 θ ) + 208 θ θ µ ( − θ + 14 θ − θ ) (cid:1) + θ m (cid:15) r (cid:0) θ θ µ (cid:0) θ − θ θ + 95 θ θ − θ θ + 70 θ (cid:1) + θ (cid:0) − θ + 70 θ θ − θ θ + 42 θ θ − θ (cid:1) + 64 θ µ (9 θ − θ + 5 θ )+ 4 θ θ µ ( − θ + 266 θ − θ ) (cid:1) + θ m (cid:15) r (cid:0) θ θ µ (cid:0) θ + 14 θ ( θ − θ ) + 15 θ (cid:1) − θ (cid:0) θ + θ θ + θ (5 θ − θ ) (cid:1) + 8 θ µ (12 θ − θ + 5 θ ) + 4 θ θ µ ( − θ + 42 θ − θ ) (cid:1) + θ (cid:15) r (cid:0) θ µ − θ θ µ + 6 θ θ µ − θ (cid:1)(cid:1) (cid:19) κ (2)4 = 33 (2) r √ µ θ µ µ (cid:18) θ µ m e (cid:15)r θ m (4Ω (5 , (cid:15) r (3 θ θ m (14 θ − θ )+ 2 θ (cid:15) r (5 θ − θ ) − θ µ (cid:15) r (5 θ − θ ))+ 2Ω (7 , (cid:15) r ( θ m (2 θ (7 θ − θ ) + 2 θ µ (5 θ − θ )) + θ (cid:15) r ( θ − θ µ ))+ θ θ m Ω (9 , (cid:15) r ( θ µ + θ ( µ − θ (cid:15) r (5 θ − θ )Ω (3 , (cid:15) r )+ 96 µ (cid:0) θ m (cid:0) θ µ (cid:0) − θ + 2 θ (35 θ − θ ) + 36 θ (7 θ − θ ) (cid:1) + θ (cid:0) − θ + 210 θ ( θ − θ ) + 56 θ θ (10 θ − θ ) + 48 θ (7 θ − θ ) (cid:1) + 192 θ µ (5 θ − θ + 8 θ ) + 24 θ θ µ (cid:0) θ ( θ − θ ) + 35 θ ( θ − θ ) + 48 θ (cid:1)(cid:1) + θ m (cid:15) r (cid:0) θ µ (cid:0) − θ + θ (70 θ − θ ) + 24 θ (14 θ − θ ) (cid:1) + 64 θ µ (5 θ − θ + 9 θ ) + 6 θ θ µ (cid:0) θ ( θ − θ ) + 35 θ ( θ − θ ) + 72 θ (cid:1) + θ (cid:0) θ (2 θ − θ ) + 35 θ (4 θ − θ ) − θ (cid:1)(cid:1) + θ m (cid:15) r (cid:0) θ µ (5 θ − θ + 12 θ )+ 8 θ θ µ ( − θ + 21 θ − θ ) + θ ( − θ + 14 θ − θ ) + 12 θ θ µ (7 θ − θ + 6 θ ) (cid:1) + θ (cid:15) r (cid:0) θ µ − θ θ µ + 6 θ θ µ − θ (cid:1)(cid:1) (cid:19) The coefficients β (2) r and β (2) f appearing the above expressions are given as β (2) f = (cid:32) µ √ πd M n n s f θ / m (cid:33) exp (cid:20) − (cid:15) f θ m (cid:21) , β (2) r = (cid:32) µ √ πd M n n s r θ / m (cid:33) exp (cid:20) − (cid:15) r θ m (cid:21) (101)
8. Supplementary material: Contribution from ∆ α We can now look into the contribution arising from ∆ α into the production terms mentioned above. Since ∆ α is a non-equilibrium quantity and in the present work we are only dealing with linearised production terms so theproduction terms corresponding to u ( α ) i , q ( α ) i and σ ( α ) ij will not see any contribution from ∆ α . Only the reactionsrates, the mass balance and the energy balance will be influenced by ∆ α . The generic form of the contribution from ∆ α into the mass production is given as:∆ r ( α )0 = (cid:88) γ =1 (cid:16) φ ( α ) γ ∆ γ (cid:17) (102) φ (1)1 = µ β f θ (cid:18) θ m + 2 θ m (cid:15) f − (cid:15) f (cid:19) , φ (1)2 = θ µ β f (cid:18) θ m + 2 θ m (cid:15) f − (cid:15) f (cid:19) φ (1)3 = µ β r θ (cid:18) − θ m − θ m (cid:15) r + (cid:15) r (cid:19) , φ (1)4 = θ µ β r (cid:18) − θ m − θ m (cid:15) r + (cid:15) r (cid:19) φ (2)1 = µ β f θ (cid:18) θ m + 2 θ m (cid:15) f − (cid:15) f (cid:19) , φ (2)2 = θ µ β f (cid:18) θ m + 2 θ m (cid:15) f − (cid:15) f (cid:19) φ (2)3 = µ β r θ (cid:18) − θ m − θ m (cid:15) r + (cid:15) r (cid:19) , φ (2)4 = θ µ β r (cid:18) − θ m − θ m (cid:15) r + (cid:15) r (cid:19) .1.3. Component-3 φ (3)1 = 1120 β f θ µ (cid:0) − θ m − θ m (cid:15) f + (cid:15) f (cid:1) , φ (3)2 = 1120 β f θ µ (cid:0) − θ m − θ m (cid:15) f + (cid:15) f (cid:1) φ (3)3 = 1120 β r θ µ (cid:0) θ m + 2 θ m (cid:15) r − (cid:15) r (cid:1) , φ (3)4 = 1120 β r θ µ (cid:0) θ m + 2 θ m (cid:15) r − (cid:15) r (cid:1) φ (4)1 = 1120 β f θ µ (cid:0) − θ m − θ m (cid:15) f + (cid:15) f (cid:1) , φ (4)2 = 1120 β f θ µ (cid:0) − θ m − θ m (cid:15) f + (cid:15) f (cid:1) φ (4)3 = 1120 β r θ µ (cid:0) θ m + 2 θ m (cid:15) r − (cid:15) r (cid:1) , φ (4)4 = 1120 β r θ µ (cid:0) θ m + 2 θ m (cid:15) r − (cid:15) r (cid:1) The contribution from the 14 th -moment into the energy production rate, given by ∆ r ( α )1 in the above sections,can be written in a generic form as ∆ r ( α )1 = (cid:88) γ =1 (cid:16) ψ ( α ) γ ∆ γ (cid:17) (103) ψ (1)1 = − µ β f θ θ m (cid:18) θ m (cid:0) θ +29 θ θ +20 θ (cid:1) + θ m (cid:15) f (cid:0) θ +18 θ θ +20 θ (cid:1) + θ θ m (cid:15) f (2 θ +11 θ )+ θ (cid:15) f (cid:19) ψ (1)2 = µ β f θ θ θ m (cid:18) θ θ m − θ m (cid:15) f ( θ − θ ) + 3 θ m (cid:15) f (2 θ − θ ) − θ (cid:15) f (cid:19) ψ (1)3 = β r θ µ θ m (cid:18) θ m e (cid:15)r θ m (20Ω (5 , (cid:15) r (3 θ m + 2 (cid:15) r ) − (7 , (cid:15) r (10 θ m + (cid:15) r ) + θ m Ω (9 , (cid:15) r − (cid:15) r Ω (3 , (cid:15) r ) + 64 (cid:0) θ m (cid:0) θ ( µ − + θ θ (8 µ (3 µ −
7) + 29) + 4 θ ( µ − µ − (cid:1) + θ m (cid:15) r (cid:0) θ ( µ − + 2 θ θ ( µ (7 µ −
19) + 9) + θ ( µ − µ − (cid:1) + θ m (cid:15) r (cid:0) − θ µ ( θ + 2 θ ) + θ (2 θ + 11 θ ) + 8 θ µ (cid:1) + (cid:15) r ( θ ( µ −
1) + θ µ ) (cid:1) (cid:19) ψ (1)4 = β r θ µ θ m (cid:18) θ m e (cid:15)r θ m (20Ω (5 , (cid:15) r (3 θ m +2 (cid:15) r ) − (7 , (cid:15) r (10 θ m + (cid:15) r )+ θ m Ω (9 , (cid:15) r − (cid:15) r Ω (3 , (cid:15) r )+ 64 (cid:0) θ m (cid:0) θ θ µ − θ θ + 48 θ µ (cid:1) + θ m (cid:15) r (cid:0) θ θ µ + 3 θ ( θ − θ ) + 28 θ µ (cid:1) + θ m (cid:15) r (cid:0) θ θ µ + 3 θ ( θ − θ ) + 8 θ µ (cid:1) + (cid:15) r ( θ ( µ −
1) + θ µ ) (cid:1) (cid:19) .2.2. Component-2 ψ (2)1 = β f θ θ µ (cid:16) θ θ m + 3 θ m (cid:15) f (2 θ − θ ) − θ m (cid:15) f ( θ − θ ) − θ (cid:15) f (cid:17) θ m ψ (2)2 = − µ β f θ θ m (cid:18) θ m (cid:0) θ +29 θ θ +12 θ (cid:1) + θ m (cid:15) f (cid:0) θ +18 θ θ +7 θ (cid:1) + θ θ m (cid:15) f (11 θ +2 θ )+ θ (cid:15) f (cid:19) ψ (2)3 = β r θ θ µ m (cid:18) θ µ m e (cid:15)r θ m (20Ω (5 , (cid:15) r (3 θ m + 2 (cid:15) r ) − (7 , (cid:15) r (10 θ m + (cid:15) r ) + θ m Ω (9 , (cid:15) r − (cid:15) r Ω (3 , (cid:15) r )+ 64 µ (cid:0) θ m (cid:0) − θ θ + 48 θ µ + 16 θ θ µ (cid:1) + θ m (cid:15) r (cid:0) θ ( θ − θ ) + 28 θ µ + 20 θ θ µ (cid:1) + θ m (cid:15) r (cid:0) θ ( θ − θ ) + 8 θ µ + 8 θ θ µ (cid:1) + (cid:15) r ( θ µ + θ ( µ − (cid:1) (cid:19) ψ (2)4 = β r θ θ µ m (cid:18) θ µ m e (cid:15)r θ m (20Ω (5 , (cid:15) r (3 θ m + 2 (cid:15) r ) − (7 , (cid:15) r (10 θ m + (cid:15) r ) + θ m Ω (9 , (cid:15) r − (cid:15) r Ω (3 , (cid:15) r ) + 64 µ (cid:0) θ m (cid:0) θ ( µ − µ −
5) + θ θ (8 µ (3 µ −
7) + 29) + 12 θ ( µ − (cid:1) + θ m (cid:15) r (cid:0) θ ( µ − µ −
10) + 2 θ θ ( µ (7 µ −
19) + 9) + 7 θ ( µ − (cid:1) + θ m (cid:15) r (cid:0) − θ µ (2 θ + θ ) + θ (11 θ + 2 θ ) + 8 θ µ (cid:1) + (cid:15) r ( θ µ + θ ( µ − (cid:1) (cid:19) ψ (3)1 = β f θ µ θ m (cid:18) θ m e (cid:15)f θ m (20Ω (5 , (cid:15) f (3 θ m + 2 (cid:15) f ) − (7 , (cid:15) f (10 θ m + (cid:15) f ) + θ m Ω (9 , (cid:15) f − (cid:15) f Ω (3 , (cid:15) f ) + 64 (cid:0) θ m (cid:0) θ ( µ − + θ θ (8 µ (3 µ −
7) + 29) + 4 θ ( µ − µ − (cid:1) + θ m (cid:15) f (cid:0) θ ( µ − + 2 θ θ ( µ (7 µ −
19) + 9) + θ ( µ − µ − (cid:1) + θ m (cid:15) f (cid:0) − θ µ ( θ + 2 θ ) + θ (2 θ + 11 θ ) + 8 θ µ (cid:1) + (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:1) (cid:19) ψ (3)2 = β f θ µ θ m (cid:18) θ m e (cid:15)f θ m (20Ω (5 , (cid:15) f (3 θ m + 2 (cid:15) f ) − (7 , (cid:15) f (10 θ m + (cid:15) f ) + θ m Ω (9 , (cid:15) f − (cid:15) f Ω (3 , (cid:15) f )+ 64 (cid:0) θ m (cid:0) θ θ µ − θ θ + 48 θ µ (cid:1) + θ m (cid:15) f (cid:0) θ θ µ + 3 θ ( θ − θ ) + 28 θ µ (cid:1) + θ m (cid:15) f (cid:0) θ θ µ + 3 θ ( θ − θ ) + 8 θ µ (cid:1) + (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:1) (cid:19) (3)3 = − µ β r θ θ m (cid:18) θ m (cid:0) θ + 29 θ θ + 20 θ (cid:1) + θ m (cid:15) r (cid:0) θ + 18 θ θ + 20 θ (cid:1) + θ θ m (cid:15) r (2 θ + 11 θ ) + θ (cid:15) r (cid:19) ψ (3)4 = µ β r θ θ θ m (cid:18) θ θ m − θ m (cid:15) r ( θ − θ ) + 3 θ m (cid:15) r (2 θ − θ ) − θ (cid:15) r (cid:19) ψ (4)1 = β f θ θ µ m (cid:18) θ µ m e (cid:15)f θ m (20Ω (5 , (cid:15) f (3 θ m + 2 (cid:15) f ) − (7 , (cid:15) f (10 θ m + (cid:15) f ) + θ m Ω (9 , (cid:15) f − (cid:15) f Ω (3 , (cid:15) f )+ 64 µ (cid:0) θ m (cid:0) − θ θ + 48 θ µ + 16 θ θ µ (cid:1) + θ m (cid:15) f (cid:0) θ ( θ − θ ) + 28 θ µ + 20 θ θ µ (cid:1) + θ m (cid:15) f (cid:0) θ ( θ − θ ) + 8 θ µ + 8 θ θ µ (cid:1) + (cid:15) f ( θ µ + θ ( µ − (cid:1) (cid:19) ψ (4)2 = β f θ θ µ m (cid:18) θ µ m e (cid:15)f θ m (20Ω (5 , (cid:15) f (3 θ m + 2 (cid:15) f ) − (7 , (cid:15) f (10 θ m + (cid:15) f ) + θ m Ω (9 , (cid:15) f − (cid:15) f Ω (3 , (cid:15) f ) + 64 µ (cid:0) θ m (cid:0) θ ( µ − µ −
5) + θ θ (8 µ (3 µ −
7) + 29) + 12 θ ( µ − (cid:1) + θ m (cid:15) f (cid:0) θ ( µ − µ −
10) + 2 θ θ ( µ (7 µ −
19) + 9) + 7 θ ( µ − (cid:1) + θ m (cid:15) f (cid:0) − θ µ (2 θ + θ ) + θ (11 θ + 2 θ ) + 8 θ µ (cid:1) + (cid:15) f ( θ µ + θ ( µ − (cid:1) (cid:19) ψ (4)3 = µ β r θ θ θ m (cid:18) θ θ m + 3 θ m (cid:15) r (2 θ − θ ) − θ m (cid:15) r ( θ − θ ) − θ (cid:15) r (cid:19) ψ (4)4 = − µ β r θ θ m (cid:18) θ m (cid:0) θ + 29 θ θ + 12 θ (cid:1) + θ m (cid:15) r (cid:0) θ + 18 θ θ + 7 θ (cid:1) + θ θ m (cid:15) r (11 θ + 2 θ ) + θ (cid:15) r (cid:19) u α Production(Elastic)
The generic form of this production term is given as: P e ( α )1 = N (cid:88) β =1 ¯ ν αβ ( λ αβ + ς α ∆ α + ς β ∆ β ) (104)37he coefficients λ αβ and ς γ can now be given as λ αβ = − √ (cid:112) θ αβ (cid:18) θ α θ β µ β (cid:0) − µ α + 6 µ α µ β − µ β (cid:1) − θ β µ β + 2 θ α µ α (cid:0) µ α + 2 µ α µ β + 3 µ β (cid:1) + θ α θ β (cid:0) µ α − µ α µ β + 9 µ α µ β − µ β (cid:1) (cid:19) ς α = − θ α √ θ αβ ) / (cid:18) θ α µ α (cid:0) µ α + 2 µ α µ β + 3 µ β (cid:1) + θ α θ β (cid:0) µ α + 69 µ α µ β + 210 µ α µ β − µ β (cid:1) + 3 θ α θ β (cid:0) µ α + 44 µ α µ β + 81 µ α µ β − µ β (cid:1) + 4 θ β (cid:0) µ α − µ α µ β + 12 µ α µ β − µ β (cid:1) (cid:19) ς β = θ β √ θ αβ ) / (cid:18) θ α θ β µ β (cid:0) µ α − µ α µ β + 109 µ β (cid:1) + 120 θ β µ β + 2 θ α (cid:0) µ α + 6 µ α µ β − µ α µ β + 40 µ β (cid:1) + θ α θ β (cid:0) µ α + 36 µ α µ β − µ α µ β + 290 µ β (cid:1) (cid:19) u α Production(Chemical)
The generic form of the production term of the fourteenth moment can be given as P r ( α )2 = (cid:88) γ =1 (cid:16) ε ( α ) γ ∆ α (cid:17) + χ ( α ) (105)The coefficients ε ( α ) γ and χ ( α ) can now be given as χ (1) = β r µ (cid:18) θ m θ e (cid:15)r θ m (cid:16) (5 , (cid:15) r (cid:0) θ θ θ m + (cid:15) r (cid:0) − θ − θ µ θ + 2 θ µ θ (cid:1)(cid:1) + θ (cid:16) m Ω (7 , (cid:15) r (cid:0) θ + 2 θ µ θ − θ µ θ (cid:1) + 3 θ m θ Ω (7 , (cid:15) r − (cid:15) r θ Ω (5 , (cid:15) r (cid:17) − θ θ (cid:15) r Ω (3 , (cid:15) r (cid:17) + 24 (cid:0) θ (cid:0) θ m (cid:0) θ + 40 θ θ + 15 θ (cid:1) + 2 θ θ m (cid:15) r (4 θ + 5 θ ) + θ (cid:15) r (cid:1) − θ µ θ (cid:0) θ m (6 θ + 5 θ ) + θ m (cid:15) r (8 theta + 5 θ ) + θ (cid:15) r (cid:1) + 4 θ µ θ (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) − θ θ µ θ (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) + 4 θ θ µ θ (cid:0) θ (cid:0) θ m + 8 θ m (cid:15) r + (cid:15) r (cid:1) + 5 θ θ m (4 θ m + (cid:15) r ) (cid:1)(cid:1) (cid:19) − β f θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 15 θ (cid:1) + 2 θ θ m (cid:15) f (4 θ + 5 θ ) + θ (cid:15) f (cid:1) ε (1)1 = 38 β (3) f θ µ (cid:0) θ m (cid:0) θ + 424 θ θ + 539 θ θ + 280 θ θ + 40 θ (cid:1) + 6 θ θ m (cid:15) f (7 θ + 12 θ ) (cid:0) θ + 7 θ θ + 5 θ (cid:1) + 3 θ θ m (cid:15) f (cid:0) θ + 44 θ θ + 61 θ (cid:1) + 2 θ θ m (cid:15) f (3 θ + 13 θ ) + θ (cid:15) f (cid:1) ε (1)2 = − β (3) f θ θ µ (cid:0) − θ θ m + 30 θ θ m (cid:15) f ( θ − θ ) + 15 θ m (cid:15) f (cid:0) θ − θ θ + θ (cid:1) + 10 θ θ m (cid:15) f ( θ − θ ) + θ (cid:15) f (cid:1) ε (1)3 = β (3) r θ µ (cid:18) (cid:0) θ m (cid:0) θ ( µ − + 24 θ θ ( µ − (12 µ (5 µ −
9) + 53)+ θ θ (16 µ ( µ (27 µ (5 µ −
18) + 664) − θ θ ( µ − µ (36 µ (5 µ −
12) + 343) − θ ( µ − µ (9 µ (5 µ −
11) + 65) − (cid:1) + 2 θ m (cid:15) r (cid:0) θ µ (cid:0) θ + 803 θ θ + 252 θ (cid:1) + 3 θ (7 θ + 12 θ ) (cid:0) θ + 7 θ θ + 5 θ (cid:1) − θ µ (cid:0) θ + 419 θ θ + 312 θ θ + 40 θ (cid:1) − θ µ (21 θ + 16 θ ) + 1344 θ µ (cid:1) + θ m (cid:15) r (cid:0) θ µ (cid:0) theta + 178 θ θ + 44 θ (cid:1) − θ θ µ (cid:0) θ + 94 θ θ + 64 θ (cid:1) + 3 θ (cid:0) θ + 44 θ θ + 61 θ (cid:1) − θ µ (39 θ + 28 θ ) + 624 θ µ (cid:1) + 2 θ m (cid:15) r ( θ ( µ −
1) + θ µ ) (cid:0) − θ µ (3 θ + 4 θ ) + θ (3 θ + 13 θ ) + 12 θ µ (cid:1) + ( θ (cid:15) r − θ µ (cid:15) r ) (cid:1) +4 θ m e (cid:15)r θ m (cid:16) − (5 , (cid:15) r (cid:0) θ θ m (4 θ − θ )+ (cid:15) r (cid:0) θ ( µ − + 2 θ θ (cid:0) µ + µ − (cid:1) + θ ( µ + 2)(3 µ + 2) (cid:1)(cid:1) + 600 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 10 θ θ µ m Ω (11 , (cid:15) r − θ θ µ m Ω (7 , (cid:15) r + 400 θ θ µ m Ω (9 , (cid:15) r − θ θ µ m Ω (11 , (cid:15) r + 600 θ θ m Ω (7 , (cid:15) r + 180 θ θ m Ω (7 , (cid:15) r − θ θ m Ω (9 , (cid:15) r − θ θ m Ω (9 , (cid:15) r + 10 θ θ m Ω (11 , (cid:15) r + 1200 θ θ θ µ m Ω (7 , (cid:15) r − θ θ θ µ m Ω (9 , (cid:15) r + 20 θ θ θ µ m Ω (11 , (cid:15) r + 400 θ θ θ µ m Ω (7 , (cid:15) r + 240 θ θ θ µ m Ω (9 , (cid:15) r − θ θ θ µ m Ω (11 , (cid:15) r − θ θ θ m Ω (7 , (cid:15) r + 360 θ θ θ m Ω (7 , (cid:15) r + 220 θ θ θ m Ω (9 , (cid:15) r − θ θ θ m Ω (9 , (cid:15) r + 12 θ m Ω (11 , (cid:15) r + 600 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 10 θ θ µ m Ω (11 , (cid:15) r + 1600 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 800 θ θ m Ω (7 , (cid:15) r + 180 θ θ m Ω (7 , (cid:15) r − θ θ m Ω (9 , (cid:15) r + 400 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ µ (cid:15) r Ω (7 , (cid:15) r + 40 θ µ (cid:15) r Ω (9 , (cid:15) r − θ (cid:15) r Ω (5 , (cid:15) r + 400 θ (cid:15) r Ω (7 , (cid:15) r + 120 θ (cid:15) r Ω (7 , (cid:15) r − θ (cid:15) r Ω (9 , (cid:15) r − θ (cid:15) r Ω (9 , (cid:15) r + 800 θ theta µ (cid:15) r Ω (7 , (cid:15) r − θ θ µ (cid:15) r Ω (9 , (cid:15) r − θ θ µ (cid:15) r Ω (7 , (cid:15) r + 40 θ θ µ (cid:15) r Ω (9 , (cid:15) r + 2400 θ (cid:15) r (4 θ − θ )Ω (3 , (cid:15) r − θ θ (cid:15) r Ω (5 , (cid:15) r − θ θ (cid:15) r Ω (7 , (cid:15) r + 240 θ θ (cid:15) r Ω (7 , (cid:15) r − θ θ (cid:15) r Ω (9 , (cid:15) r + 400 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r + 320 θ µ (cid:15) r Ω (7 , (cid:15) r − θ (cid:15) r Ω (5 , (cid:15) r + 120 θ (cid:15) r Ω (7 , (cid:15) r − θ (cid:15) r Ω (9 , (cid:15) r (cid:17) (cid:19) ε (1)4 = 39 (3) r θ µ (cid:18) (cid:0) θ m (cid:0) θ θ µ (2 θ + 5 θ ) − θ θ − θ θ µ + 32 θ θ µ (2 θ + 15 θ ) + 5760 θ µ (cid:1) + 2 θ m (cid:15) r (cid:0) θ θ µ (4 θ + 25 θ ) + 15 θ θ ( θ − θ ) − θ θ µ − θ θ µ (12 θ − θ ) + 1344 θ µ (cid:1) + θ m (cid:15) r (cid:0) θ θ µ ( θ +10 θ )+15 θ (cid:0) θ − θ θ + θ (cid:1) − θ θ µ − θ θ µ (7 θ − θ )+624 θ µ (cid:1) + 2 θ m (cid:15) r ( θ ( µ −
1) + θ µ ) (cid:0) θ θ µ + 5 θ ( θ − θ ) + 12 θ µ (cid:1) + ( θ (cid:15) r − θ µ (cid:15) r ) (cid:1) +4 θ m e (cid:15)r θ m (cid:16) θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 10 θ θ µ m Ω (11 , (cid:15) r − θ θ µ m Ω (7 , (cid:15) r +560 θ θ µ m Ω (9 , (cid:15) r − θ θ µ m Ω (11 , (cid:15) r +3000 θ θ m Ω (7 , (cid:15) r + 180 θ θ m Ω (7 , (cid:15) r − θ θ m Ω (9 , (cid:15) r − θ θ m Ω (9 , (cid:15) r + 10 θ θ m Ω (11 , (cid:15) r − (5 , (cid:15) r (cid:0) θ (4 θ θ m − θ θ m + 5 θ (cid:15) r − θ (cid:15) r ) − θ θ µ (cid:15) r + 12 θ µ (cid:15) r (cid:1) + 1200 θ θ θ µ m Ω (7 , (cid:15) r − θ θ θ µ m Ω (9 , (cid:15) r + 20 θ θ θ µ m Ω (11 , (cid:15) r − θ θ θ µ m Ω (7 , (cid:15) r + 560 θ θ θ µ m Ω (9 , (cid:15) r − θ θ θ µ m Ω (11 , (cid:15) r − θ θ θ m Ω (7 , (cid:15) r + 360 θ θ θ m Ω (7 , (cid:15) r + 60 θ θ θ m Ω (9 , (cid:15) r − θ θ θ m Ω (9 , (cid:15) r + 12 θ m Ω (11 , (cid:15) r + 600 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 10 θ θ µ m Ω (11 , (cid:15) r + 180 θ θ m Ω (7 , (cid:15) r − θ θ m Ω (9 , (cid:15) r + 400 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ µ (cid:15) r Ω (7 , (cid:15) r + 40 θ µ (cid:15) r Ω (9 , (cid:15) r − θ (cid:15) r Ω (5 , (cid:15) r + 720 θ (cid:15) r Ω (7 , (cid:15) r + 120 θ (cid:15) r Ω (7 , (cid:15) r − θ (cid:15) r Ω (9 , (cid:15) r − θ (cid:15) r Ω (9 , (cid:15) r + 800 θ θ µ (cid:15) r Ω (7 , (cid:15) r − θ θ µ (cid:15) r Ω (9 , (cid:15) r − θ θ µ (cid:15) r Ω (7 , (cid:15) r + 40 theta θ µ (cid:15) r Ω (9 , (cid:15) r + 2400 θ (cid:15) r (4 θ − θ )Ω (3 , (cid:15) r − θ θ (cid:15) r Ω (5 , (cid:15) r − θ θ (cid:15) r Ω (7 , (cid:15) r + 240 θ θ (cid:15) r Ω (7 , (cid:15) r − θ θ (cid:15) r Ω (9 , (cid:15) r + 400 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ (cid:15) r Ω (5 , (cid:15) r + 120 θ (cid:15) r Ω (7 , (cid:15) r − θ (cid:15) r Ω (9 , (cid:15) r (cid:17) (cid:19) χ (2) = − β f θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 24 θ (cid:1) + 2 θ θ m (cid:15) f (5 θ + 4 θ ) + θ (cid:15) f (cid:1) + β r µ (cid:18) θ µ m e (cid:15)r θ m (cid:16) µ Ω (5 , (cid:15) r (cid:0) θ θ θ m − (cid:15) r ( θ − θ µ ) (cid:1) + 12 θ µ m Ω (7 , (cid:15) r + 10 θ µ m Ω (7 , (cid:15) r ( θ − θ µ ) − θ θ µ (cid:15) r Ω (3 , (cid:15) r − θ µ (cid:15) r Ω (5 , (cid:15) r (cid:17) + 96 µ (cid:0) θ m (cid:0) θ (cid:0) θ + 40 θ θ + 24 θ (cid:1) + 32 θ θ µ (5 θ + 18 θ ) − θ θ µ (5 θ + 6 θ ) + 384 θ µ − θ θ µ (cid:1) + 2 θ m (cid:15) r ( θ µ + θ ( µ − (cid:0) θ (5 theta + 4 θ ) + 16 θ µ − θ θ µ (cid:1) + (cid:15) r ( θ µ + θ ( µ − (cid:1) (cid:19) ε (2)1 = − µ β (3) f θ θ (cid:18) − θ θ m + 30 θ θ m (cid:15) f ( θ − θ )+ 15 θ m (cid:15) f (cid:0) θ − θ θ + θ (cid:1) + 10 θ θ m (cid:15) f ( θ − θ ) + θ (cid:15) f (cid:19) ε (2)2 = 40 µ β (3) f θ (cid:18) θ m (cid:0) θ + 280 θ θ + 539 θ θ + 424 θ θ + 120 θ (cid:1) + 6 θ θ m (cid:15) f (12 θ + 7 θ ) (cid:0) θ + 7 θ θ + 4 θ (cid:1) + 3 θ θ m (cid:15) f (cid:0) θ + 44 θ θ + 13 θ (cid:1) + 2 θ θ m (cid:15) f (13 θ + 3 θ ) + θ (cid:15) f (cid:19) ε (2)3 = β (3) r θ µ (cid:18) µ (cid:0) θ m (cid:0) − θ θ + 32 θ θ µ (15 θ + 2 θ )+ 32 θ θ µ (5 θ + 2 θ ) + 5760 θ µ − θ θ µ (cid:1) + 2 θ m (cid:15) r (cid:0) θ θ µ (35 θ − θ ) + 4 θ θ µ (25 θ + 4 θ ) + 15 θ θ ( θ − θ ) + 1344 θ µ − θ θ µ (cid:1) + θ m (cid:15) r (cid:0) θ (cid:0) θ − θ θ + θ (cid:1) +8 θ θ µ (10 θ − θ )+8 θ θ µ (10 θ + θ )+624 θ µ − θ θ µ (cid:1) + 2 θ m (cid:15) r ( θ − θ µ ) (cid:0) θ ( θ − θ ) + 12 θ µ + 4 θ θ µ (cid:1) + (cid:15) r ( θ − θ µ ) (cid:1) θ µ m e (cid:15)r θ m (cid:16) θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r − θ θ µ µ m Ω (9 , (cid:15) r + 10 θ θ µ µ m Ω (11 , (cid:15) r + 40 µ Ω (7 , (cid:15) r (cid:0) θ ( − θ θ m + 25 θ θ m + θ ( − (cid:15) r ) + 6 θ (cid:15) r )+ 20 θ µ (3 θ m + 2 (cid:15) r ) − θ θ µ (5 θ m + 2 (cid:15) r ) (cid:1) + 360 θ θ θ µ m Ω (7 , (cid:15) r − θ θ θ µ m Ω (9 , (cid:15) r − θ θ θ µ µ m Ω (9 , (cid:15) r + 20 θ θ theta µ µ m Ω (11 , (cid:15) r + 560 θ θ θ µ µ m Ω (9 , (cid:15) r − θ θ θ µ µ m Ω (11 , (cid:15) r + 60 θ θ θ µ m Ω (9 , (cid:15) r + 12 θ µ m Ω (11 , (cid:15) r + 180 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r − θ θ µ µ m Ω (9 , (cid:15) r + 10 θ θ µ µ m Ω (11 , (cid:15) r + 560 θ θ µ µ m Ω (9 , (cid:15) r − θ θ µ µ m Ω (11 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 10 θ θ µ m Ω (11 , (cid:15) r + 120 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ µ µ (cid:15) r Ω (9 , (cid:15) r + 240 θ θ µ (cid:15) r Ω (7 , (cid:15) r − θ θ µ (cid:15) r Ω (9 , (cid:15) r − θ θ µ µ (cid:15) r Ω (9 , (cid:15) r + 40 θ θ µ µ (cid:15) r Ω (9 , (cid:15) r − θ µ (cid:15) r Ω (5 , (cid:15) r + 120 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ µ µ (cid:15) r Ω (9 , (cid:15) r + 40 θ µ µ (cid:15) r Ω (9 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r (cid:17) − θ θ µ µ m (cid:15) r (3 θ − θ )Ω (3 , (cid:15) r e (cid:15)r θ m + 1600 θ µ µ m Ω (5 , (cid:15) r e (cid:15)r θ m (cid:0) θ (3 θ θ m − θ θ m + 2 θ (cid:15) r − θ (cid:15) r ) − θ µ (cid:15) r + 28 θ θ µ (cid:15) r (cid:1) (cid:19) ε (2)4 = β (3) r θ µ (cid:18) µ (cid:0) θ m (cid:0) θ ( µ − µ (9 µ (5 µ −
11) + 65) − θ θ ( µ − µ (36 µ (5 µ −
12) + 343) − θ θ (16 µ ( µ (27 µ (5 µ −
18) + 664) − θ θ ( µ − (12 µ (5 µ −
9) + 53)+360 θ ( µ − (cid:1) +2 θ m (cid:15) r (cid:0) θ µ (cid:0) θ +803 θ θ +504 θ (cid:1) +3 θ (12 θ +7 θ ) (cid:0) θ +7 θ θ +4 θ (cid:1) − θ µ (cid:0) θ + 312 θ θ + 419 θ θ + 168 θ (cid:1) − θ µ (16 θ + 21 θ ) + 1344 θ µ (cid:1) + θ m (cid:15) r (cid:0) θ µ (cid:0) θ + 178 θ θ + 117 θ (cid:1) − θ θ µ (cid:0) θ + 94 θ θ + 39 θ (cid:1) + 3 θ (cid:0) θ + 44 θ θ + 13 θ (cid:1) − θ µ (28 θ + 39 θ ) + 624 θ µ (cid:1) + 2 θ m (cid:15) r ( θ − θ µ ) (cid:0) − θ µ (4 θ + 3 θ ) + θ (13 θ + 3 θ ) + 12 θ µ (cid:1) + (cid:15) r ( θ − θ µ ) (cid:1) θ µ µ m Ω (5 , (cid:15) r e (cid:15)r θ m (cid:0) θ θ m (4 θ − θ )+ (cid:15) r (cid:0) θ ( µ + 2)(3 µ + 2) + 2 θ θ (cid:0) µ + µ − (cid:1) + 3 θ ( µ − (cid:1)(cid:1) + 4 θ µ m e (cid:15)r θ m (cid:16) µ Ω (7 , (cid:15) r (cid:0) θ m (cid:0) θ ( µ + 2)(3 µ + 2) + 2 θ θ (cid:0) µ + µ − (cid:1) + 3 θ ( µ − (cid:1) + (cid:15) r (cid:0) θ µ (2 θ − θ ) + θ (10 θ − θ ) + 40 θ µ (cid:1)(cid:1) + 180 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r − θ θ µ µ m Ω (9 , (cid:15) r + 10 θ θ µ µ m Ω (11 , (cid:15) r − θ θ µ µ m Ω (9 , (cid:15) r + 360 θ θ θ µ m Ω (7 , (cid:15) r − θ θ θ µ m Ω (9 , (cid:15) r − θ θ θ µ µ m Ω (9 , (cid:15) r + 20 θ θ θ µ µ m Ω (11 , (cid:15) r + 240 θ θ θ µ µ m Ω (9 , (cid:15) r − θ θ θ µ µ m Ω (11 , (cid:15) r + 220 θ θ θ µ m Ω (9 , (cid:15) r + 12 θ µ m Ω (11 , (cid:15) r + 180 θ θ µ m Ω (7 , (cid:15) r − θ θ µ m Ω (9 , epsilon r − θ θ µ µ m Ω (9 , (cid:15) r + 10 θ θ µ µ m Ω (11 , (cid:15) r + 400 θ θ µ µ m Ω (9 , (cid:15) r − θ θ µ µ m Ω (11 , (cid:15) r − θ θ µ m Ω (9 , (cid:15) r + 10 θ θ µ m Ω (11 , (cid:15) r + 120 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ µ µ (cid:15) r Ω (9 , (cid:15) r + 240 θ θ µ (cid:15) r Ω (7 , (cid:15) r − θ θ µ (cid:15) r Ω (9 , (cid:15) r − θ θ µ µ (cid:15) r Ω (9 , (cid:15) r + 40 θ θ µ µ (cid:15) r Ω (9 , epsilon r − θ µ (cid:15) r Ω (5 , (cid:15) r + 120 θ µ (cid:15) r Ω (7 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r − θ µ µ (cid:15) r Ω (9 , (cid:15) r + 40 θ µ µ (cid:15) r Ω (9 , (cid:15) r − θ µ (cid:15) r Ω (9 , (cid:15) r (cid:17) + +9600 θ θ µ µ m (cid:15) r (4 θ − θ )Ω (3 , (cid:15) r e (cid:15)r θ m (cid:19) χ (3) = − β r θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 15 θ (cid:1) + 2 θ θ m (cid:15) r (4 θ + 5 θ ) + θ (cid:15) r (cid:1) + β f µ (cid:18) θ m e (cid:15)f θ m (cid:16) θ m Ω (7 , (cid:15) f ( θ ( µ −
1) + θ µ ) − (5 , (cid:15) f (cid:0) (cid:15) f ( θ − θ µ ) − θ θ θ m (cid:1) + 12 θ m Ω (7 , (cid:15) f − θ θ (cid:15) f Ω (3 , (cid:15) f − θ (cid:15) f Ω (5 , (cid:15) f (cid:17) + 96 (cid:0) θ m (cid:0) − θ θ µ (6 θ + 5 θ ) + θ (cid:0) θ + 40 θ θ + 15 θ (cid:1) − θ θ µ + 32 θ θ µ (18 θ + 5 θ )+ 384 θ µ (cid:1) + 2 θ m (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:0) − θ θ µ + θ (4 θ + 5 θ ) + 16 theta µ (cid:1) + (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:1) (cid:19) ε (3)1 = β (3) f θ µ (cid:18) (cid:0) θ m (cid:0) θ ( µ − + 24 θ θ ( µ − (12 µ (5 µ −
9) + 53)+ θ θ (16 µ ( µ (27 µ (5 µ −
18) + 664) − θ θ ( µ − µ (36 µ (5 µ −
12) + 343) − θ ( µ − µ (9 µ (5 µ −
11) + 65) − (cid:1) + 2 θ m (cid:15) f (cid:0) θ µ (cid:0) θ + 803 θ θ + 252 θ (cid:1) + 3 θ (7 θ + 12 θ ) (cid:0) θ + 7 θ θ + 5 θ (cid:1) − θ µ (cid:0) θ + 419 θ θ + 312 θ θ + 40 θ (cid:1) − θ µ (21 θ + 16 θ ) + 1344 θ µ (cid:1) + θ m (cid:15) f (cid:0) θ µ (cid:0) theta + 178 θ θ + 44 θ (cid:1) − θ θ µ (cid:0) θ + 94 θ θ + 64 θ (cid:1) + 3 θ (cid:0) θ + 44 θ θ + 61 θ (cid:1) − θ µ (39 θ + 28 θ ) + 624 θ µ (cid:1) + 2 θ m (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:0) − θ µ (3 θ + 4 θ ) + θ (3 θ + 13 θ ) + 12 θ µ (cid:1) + ( θ (cid:15) f − θ µ (cid:15) f ) (cid:1) θ m e (cid:15)f θ m (cid:16) − (5 , (cid:15) f (cid:0) θ θ m (4 θ − θ )+ (cid:15) f (cid:0) θ ( µ − + 2 θ θ (cid:0) µ + µ − (cid:1) + θ ( µ + 2)(3 µ + 2) (cid:1)(cid:1) + 600 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 10 θ θ µ m Ω (11 , (cid:15) f − θ θ µ m Ω (7 , (cid:15) f + 400 θ θ µ m Ω (9 , (cid:15) f − θ θ µ m Ω (11 , (cid:15) f + 600 θ θ m Ω (7 , (cid:15) f + 180 θ θ m Ω (7 , (cid:15) f − θ θ m Ω (9 , (cid:15) f − θ θ m Ω (9 , (cid:15) f + 10 θ θ m Ω (11 , (cid:15) f + 1200 θ θ θ µ m Ω (7 , (cid:15) f − θ θ θ µ m Ω (9 , (cid:15) f + 20 θ θ θ µ m Ω (11 , (cid:15) f + 400 θ θ θ µ m Ω (7 , (cid:15) f + 240 θ θ θ µ m Ω (9 , (cid:15) f − θ θ θ µ m Ω (11 , (cid:15) f − θ θ θ m Ω (7 , (cid:15) f + 360 θ θ θ m Ω (7 , (cid:15) f + 220 θ θ θ m Ω (9 , (cid:15) f − θ θ θ m Ω (9 , (cid:15) f + 12 θ m Ω (11 , (cid:15) f + 600 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 10 θ θ µ m Ω (11 , (cid:15) f + 1600 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 800 θ θ m Ω (7 , (cid:15) f + 180 θ θ m Ω (7 , (cid:15) f − θ θ m Ω (9 , (cid:15) f + 400 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ µ (cid:15) f Ω (7 , (cid:15) f + 40 θ µ (cid:15) f Ω (9 , (cid:15) f − θ (cid:15) f Ω (5 , (cid:15) f + 400 θ (cid:15) f Ω (7 , (cid:15) f + 120 θ (cid:15) f Ω (7 , (cid:15) f − θ (cid:15) f Ω (9 , (cid:15) f − θ (cid:15) f Ω (9 , (cid:15) f + 800 θ theta µ (cid:15) f Ω (7 , (cid:15) f − θ θ µ (cid:15) f Ω (9 , (cid:15) f − θ θ µ (cid:15) f Ω (7 , (cid:15) f + 40 θ θ µ (cid:15) f Ω (9 , (cid:15) f + 2400 θ (cid:15) f (4 θ − θ )Ω (3 , (cid:15) f − θ θ (cid:15) f Ω (5 , (cid:15) f − θ θ (cid:15) f Ω (7 , (cid:15) f + 240 θ θ (cid:15) f Ω (7 , (cid:15) f − θ θ (cid:15) f Ω (9 , (cid:15) f + 400 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f + 320 θ µ (cid:15) f Ω (7 , (cid:15) f − θ (cid:15) f Ω (5 , (cid:15) f + 120 θ (cid:15) f Ω (7 , (cid:15) f − θ (cid:15) f Ω (9 , (cid:15) f (cid:17) (cid:19) ε (3)2 = β (3) f θ µ (cid:18) (cid:0) θ m (cid:0) θ θ µ (2 θ + 5 θ ) − θ θ − θ θ µ + 32 θ θ µ (2 θ + 15 θ ) + 5760 θ µ (cid:1) + 2 θ m (cid:15) f (cid:0) θ θ µ (4 θ + 25 θ ) + 15 θ θ ( θ − θ ) − θ θ µ − θ θ µ (12 θ − θ ) + 1344 θ µ (cid:1) + θ m (cid:15) f (cid:0) θ θ µ ( θ +10 θ )+15 θ (cid:0) θ − θ θ + θ (cid:1) − θ θ µ − θ θ µ (7 θ − θ )+624 θ µ (cid:1) + 2 θ m (cid:15) f ( θ ( µ −
1) + θ µ ) (cid:0) θ θ µ + 5 θ ( θ − θ ) + 12 θ µ (cid:1) + ( θ (cid:15) f − θ µ (cid:15) f ) (cid:1) +4 θ m e (cid:15)f θ m (cid:16) θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 10 θ θ µ m Ω (11 , (cid:15) f − θ θ µ m Ω (7 , (cid:15) f +560 θ θ µ m Ω (9 , (cid:15) f − θ θ µ m Ω (11 , (cid:15) f +3000 θ θ m Ω (7 , (cid:15) f + 180 θ θ m Ω (7 , (cid:15) f − θ θ m Ω (9 , (cid:15) f − θ θ m Ω (9 , (cid:15) f + 10 θ θ m Ω (11 , (cid:15) f − (5 , (cid:15) f (cid:0) θ (4 θ θ m − θ θ m + 5 θ (cid:15) f − θ (cid:15) f ) − θ theta µ (cid:15) f + 12 θ µ (cid:15) f (cid:1) + 1200 θ θ θ µ m Ω (7 , (cid:15) f − θ θ θ µ m Ω (9 , (cid:15) f + 20 θ θ θ µ m Ω (11 , (cid:15) f − θ θ θ µ m Ω (7 , (cid:15) f + 560 θ θ θ µ m Ω (9 , (cid:15) f − θ θ θ µ m Ω (11 , (cid:15) f − θ θ θ m Ω (7 , (cid:15) f + 360 θ θ θ m Ω (7 , (cid:15) f + 60 θ θ θ m Ω (9 , (cid:15) f − θ θ θ m Ω (9 , (cid:15) f + 12 θ m Ω (11 , (cid:15) f + 600 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 10 θ θ µ m Ω (11 , (cid:15) f + 180 θ θ m Ω (7 , (cid:15) f − θ θ m Ω (9 , (cid:15) f + 400 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ µ (cid:15) f Ω (7 , (cid:15) f + 40 θ µ (cid:15) f Ω (9 , (cid:15) f − θ (cid:15) f Ω (5 , (cid:15) f + 720 θ (cid:15) f Ω (7 , (cid:15) f + 120 θ (cid:15) f Ω (7 , (cid:15) f − θ (cid:15) f Ω (9 , (cid:15) f − θ (cid:15) f Ω (9 , (cid:15) f + 800 θ θ µ (cid:15) f Ω (7 , (cid:15) f − θ θ µ (cid:15) f Ω (9 , (cid:15) f − θ θ µ (cid:15) f Ω (7 , (cid:15) f + 40 theta θ µ (cid:15) f Ω (9 , (cid:15) f + 2400 θ (cid:15) f (4 θ − θ )Ω (3 , (cid:15) f − θ θ (cid:15) f Ω (5 , (cid:15) f − θ θ (cid:15) f Ω (7 , (cid:15) f + 240 θ θ (cid:15) f Ω (7 , (cid:15) f − θ θ (cid:15) f Ω (9 , (cid:15) f + 400 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ (cid:15) f Ω (5 , (cid:15) f + 120 θ (cid:15) f Ω (7 , (cid:15) f − θ (cid:15) f Ω (9 , (cid:15) f (cid:17) (cid:19) (3)3 = − µ β (3) r θ (cid:18) θ m (cid:0) θ + 424 θ θ + 539 θ θ + 280 θ θ + 40 θ (cid:1) + 6 θ θ m (cid:15) r (7 θ + 12 θ ) (cid:0) θ + 7 θ θ + 5 θ (cid:1) + 3 θ θ m (cid:15) r (cid:0) θ + 44 θ θ + 61 θ (cid:1) + 2 θ θ m (cid:15) r (3 θ + 13 θ ) + θ (cid:15) r (cid:19) ε (3)4 = − µ β (3) r θ θ (cid:18) − θ θ m + 30 θ θ m (cid:15) r ( θ − θ ) + 15 θ m (cid:15) r (cid:0) θ − θ θ + θ (cid:1) + 10 θ θ m (cid:15) r ( θ − θ ) + θ (cid:15) r (cid:19) χ (4) = β f µ (cid:18) θ µ m e (cid:15)f θ m (cid:16) µ Ω (5 , (cid:15) f (cid:0) θ θ θ m − (cid:15) f ( θ − θ µ ) (cid:1) + 12 θ µ m Ω (7 , (cid:15) f + 10 θ µ m Ω (7 , (cid:15) f ( θ − θ µ ) − θ θ µ (cid:15) f Ω (3 , (cid:15) f − θ µ (cid:15) f Ω (5 , (cid:15) f (cid:17) + 96 µ (cid:0) θ m (cid:0) θ (cid:0) θ + 40 θ θ + 24 θ (cid:1) + 32 θ θ µ (5 θ + 18 θ ) − θ θ µ (5 θ + 6 θ ) + 384 θ µ − θ θ µ (cid:1) + 2 θ m (cid:15) f ( θ µ + θ ( µ − (cid:0) θ (5 theta + 4 θ ) + 16 θ µ − θ θ µ (cid:1) + (cid:15) f ( θ µ + θ ( µ − (cid:1) (cid:19) − β r θ µ (cid:0) θ m (cid:0) θ + 40 θ θ + 24 θ (cid:1) + 2 θ θ m (cid:15) r (5 θ + 4 θ ) + θ (cid:15) r (cid:1) ε (4)1 = β (3) f θ µ (cid:18) µ (cid:0) θ m (cid:0) − θ θ + 32 θ θ µ (15 θ + 2 θ )+ 32 θ θ µ (5 θ + 2 θ ) + 5760 θ µ − θ θ µ (cid:1) + 2 θ m (cid:15) f (cid:0) θ θ µ (35 θ − θ ) + 4 θ θ µ (25 θ + 4 θ ) + 15 θ θ ( θ − θ ) + 1344 θ µ − θ θ µ (cid:1) + θ m (cid:15) f (cid:0) θ (cid:0) θ − θ θ + θ (cid:1) +8 θ θ µ (10 θ − θ )+8 θ θ µ (10 θ + θ )+624 θ µ − θ θ µ (cid:1) + 2 θ m (cid:15) f ( theta − θ µ ) (cid:0) θ ( θ − θ ) + 12 θ µ + 4 θ θ µ (cid:1) + (cid:15) f ( θ − θ µ ) (cid:1) θ µ m e (cid:15)f θ m (cid:16) − θ θ µ µ m Ω (9 , (cid:15) f + 10 θ θ µ µ m Ω (11 , (cid:15) f + 180 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 40 µ Ω (7 , (cid:15) f (cid:0) θ ( − θ θ m + 25 θ θ m + θ ( − (cid:15) f ) + 6 θ (cid:15) f ) + 20 θ µ (3 θ m + 2 (cid:15) f ) − θ θ µ (5 θ m + 2 (cid:15) f ) (cid:1) − θ θ θ µ µ m Ω (9 , (cid:15) f + 20 θ θ θ µ µ m Ω (11 , (cid:15) f + 560 θ θ θ µ µ m Ω (9 , (cid:15) f − theta θ θ µ µ m Ω (11 , (cid:15) f + 360 θ θ θ µ m Ω (7 , (cid:15) f − θ θ θ µ m Ω (9 , (cid:15) f + 60 θ θ θ µ m Ω (9 , (cid:15) f + 12 θ µ m Ω (11 , (cid:15) f − θ θ µ µ m Ω (9 , (cid:15) f + 10 θ θ µ µ m Ω (11 , (cid:15) f + 560 θ θ µ µ m Ω (9 , (cid:15) f − θ θ µ µ m Ω (11 , (cid:15) f + 180 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 10 θ θ µ m Ω (11 , (cid:15) f − θ µ µ (cid:15) f Ω (9 , (cid:15) f + 120 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ θ µ µ (cid:15) f Ω (9 , (cid:15) f + 40 θ θ µ µ (cid:15) f Ω (9 , (cid:15) f + 240 θ θ µ (cid:15) f Ω (7 , (cid:15) f − θ θ µ (cid:15) f Ω (9 , (cid:15) f − θ µ (cid:15) f Ω (5 , (cid:15) f − θ µ µ (cid:15) f Ω (9 , (cid:15) f + 40 θ µ µ (cid:15) f Ω (9 , (cid:15) f + 120 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f (cid:17) − θ θ µ µ m (cid:15) f (3 θ − θ )Ω (3 , (cid:15) f e (cid:15)f θ m + 1600 θ µ µ m Ω (5 , (cid:15) f e (cid:15)f θ m (cid:0) θ (3 θ θ m − θ θ m + 2 θ (cid:15) f − θ (cid:15) f ) − θ µ (cid:15) f + 28 θ θ µ (cid:15) f (cid:1) (cid:19) ε (4)2 = β (3) f θ µ (cid:18) µ (cid:0) θ m (cid:0) θ ( µ − µ (9 µ (5 µ −
11) + 65) − θ θ ( µ − µ (36 µ (5 µ −
12) + 343) − θ θ (16 µ ( µ (27 µ (5 µ −
18) + 664) − θ θ ( µ − (12 µ (5 µ −
9) + 53)+360 θ ( µ − (cid:1) +2 θ m (cid:15) f (cid:0) θ µ (cid:0) θ +803 θ θ +504 θ (cid:1) +3 θ (12 θ +7 θ ) (cid:0) θ +7 θ θ +4 θ (cid:1) − θ µ (cid:0) θ + 312 θ θ + 419 θ θ + 168 θ (cid:1) − θ µ (16 θ + 21 θ ) + 1344 θ µ (cid:1) + θ m (cid:15) f (cid:0) θ µ (cid:0) θ + 178 θ θ + 117 θ (cid:1) − θ θ µ (cid:0) θ + 94 θ θ + 39 θ (cid:1) + 3 θ (cid:0) θ + 44 θ θ + 13 θ (cid:1) − θ µ (28 θ + 39 θ ) + 624 θ µ (cid:1) + 2 θ m (cid:15) f ( θ − θ µ ) (cid:0) − θ µ (4 θ + 3 θ ) + θ (13 θ + 3 θ ) + 12 θ µ (cid:1) + (cid:15) f ( θ − θ µ ) (cid:1) θ µ µ m Ω (5 , (cid:15) f e (cid:15)f θ m (cid:0) θ θ m (4 θ − θ )+ (cid:15) f (cid:0) θ ( µ + 2)(3 µ + 2) + 2 θ θ (cid:0) µ + µ − (cid:1) + 3 θ ( µ − (cid:1)(cid:1) + 4 θ µ m e (cid:15)f θ m (cid:16) µ Ω (7 , (cid:15) f (cid:0) θ m (cid:0) θ ( µ + 2)(3 µ + 2) + 2 θ θ (cid:0) µ + µ − (cid:1) + 3 θ ( µ − (cid:1) + (cid:15) f (cid:0) θ µ (2 θ − θ ) + θ (10 θ − θ ) + 40 θ µ (cid:1)(cid:1) − θ θ µ µ m Ω (9 , (cid:15) f + 10 θ θ µ µ m Ω (11 , (cid:15) f − θ θ µ µ m Ω (9 , (cid:15) f + 180 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f − θ θ θ µ µ m Ω (9 , (cid:15) f + 20 θ θ θ µ µ m Ω (11 , (cid:15) f + 240 θ θ θ µ µ m Ω (9 , (cid:15) f − θ θ θ µ µ m Ω (11 , (cid:15) f + 360 θ θ θ µ m Ω (7 , (cid:15) f − θ θ θ µ m Ω (9 , (cid:15) f + 220 θ θ θ µ m Ω (9 , (cid:15) f + 12 θ µ m Ω (11 , (cid:15) f − θ θ µ µ m Ω (9 , (cid:15) f + 10 θ θ µ µ m Omega (11 , (cid:15) f + 400 θ θ µ µ m Ω (9 , (cid:15) f − θ θ µ µ m Ω (11 , (cid:15) f + 180 θ θ µ m Ω (7 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f − θ θ µ m Ω (9 , (cid:15) f + 10 θ θ µ m Ω (11 , (cid:15) f − θ µ µ (cid:15) f Ω (9 , (cid:15) f + 120 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ θ µ µ (cid:15) f Ω (9 , (cid:15) f + 40 θ θ µ µ (cid:15) f Ω (9 , (cid:15) f + 240 θ θ µ (cid:15) f Ω (7 , (cid:15) f − θ θ µ (cid:15) f Ω (9 , (cid:15) f − θ µ (cid:15) f Ω (5 , (cid:15) f − θ µ µ (cid:15) f Ω (9 , (cid:15) f + 40 θ µ µ (cid:15) f Ω (9 , (cid:15) f + 120 θ µ (cid:15) f Ω (7 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f − θ µ (cid:15) f Ω (9 , (cid:15) f (cid:17) + +9600 θ θ µ µ m (cid:15) f (4 θ − θ )Ω (3 , (cid:15) f e (cid:15)f θ m (cid:19) ε (4)3 = − µ β (3) r θ θ ( − θ θ m + 30 θ θ m (cid:15) r ( θ − θ )+ 15 θ m (cid:15) r (cid:0) θ − θ θ + θ (cid:1) + 10 θ θ m (cid:15) r ( θ − θ ) + θ (cid:15) r (cid:19) ε (4)4 = − µ β (3) r θ (cid:18) θ m (cid:0) θ + 280 θ θ + 539 θ θ + 424 θ θ + 120 θ (cid:1) + 6 θ θ m (cid:15) r (12 θ + 7 θ ) (cid:0) θ + 7 θ θ + 4 θ (cid:1) + 3 θ θ m (cid:15) r (cid:0) θ + 44 θ θ + 13 θ (cid:1) + 2 θ θ m (cid:15) r (13 θ + 3 θ ) + θ (cid:15) r (cid:19) The coefficients β (3) f and β (3) r are given as: β (3) f = (cid:32) √ πd M n n sf θ / m (cid:33) exp (cid:20) − (cid:15) f θ m (cid:21) , β (3) r = (cid:32) √ πd M n n sr θ / m (cid:33) exp (cid:20) − (cid:15) r θ m (cid:21)(cid:21)