A Two-fluid Model for Plasma with Prandtl Number Correction
aa r X i v : . [ phy s i c s . p l a s m - ph ] M a y A Two-fluid Model for Plasma with Prandtl Number Correction
Ruo Li ∗ , Yixiao Lu † , Yanli Wang ‡ May 14, 2020
Abstract
A two-fluid model is derived from the plasma kinetic equations using the moment modelreduction method. The moment method we adopt was recently developed with a glob-ally hyperbolic regularization where the moment model attained is locally well-posed intime. Based on the hyperbolic moment model with well-posedness, the Maxwellian iterationmethod is utilized to get the closure relations for the resulted two-fluid model. By takingthe Shakhov collision operator in the Maxwellian iteration, the two-fluid model inherits thecorrect Prandtl number from the plasma kinetic equations. The new model is formally thesame as the five-moment two-fluid model except for the closure relations, where the pressuretensor is anisotropic and the heat flux is presented. This provides the model the capacity todepict problems with anisotropic pressure tensor and large heat flux.
Keywords:
Plasma kinetic equations; Maxwellian iteration; two-fluid model; moment modelreduction
The multi-species plasma model was proposed to capture dynamics in systems with multiplespecies of charges [17]. Each species is treated separately and coupled together through elec-tromagnetic and collisional processes. In the collisionless problems, the evolution may be mostaccurately [24, 31] modeled using the kinetic theory such as the Vlasov-Maxwell equations, wherethe plasma is described by the distribution functions of all species, and the electromagnetic fieldsare governed by Maxwell equations. However, it is extremely expensive to numerically solve theplasma kinetic equations due to the high dimensionality of the plasma kinetic equations[26].Several reduced models were introduced [10, 21, 11] to approximate the plasma kinetic equa-tions. When plasma behaves like an electrically conducting fluid, where the motion of theelectrons and ions is locked together by electrostatic forces, the MHD model is utilized [8]. Inthe MHD model, the plasma is treated as a single electrically conducting fluid. Several algo-rithms have been designed based on the MHD model, which are already used to simulate severalplasma phenomena successfully [16, 23]. It was pointed out that the MHD model ignores theelectron mass and the finite Larmor radius effects [5]. This may lead to the limitation to treatthe Hall effect and diamagnetic terms [14]. Hall effect can be added to the MHD equationsand we get the Hall-MHD model [12]. Though a distinction is made between the bulk plasma ∗ CAPT, LMAM & School of Mathematical Sciences, Peking University, Beijing, China, email: [email protected] . † School of Mathematical Sciences, Peking University, Beijing, China, 100871, email: [email protected] . ‡ Beijing Computational Science Research Center, email: [email protected] . In the collisionless plasma where the collective interactions dominate the plasma dynamics, theVlasov-Maxwell (VM) equations are mostly used to describe the evolution of the plasma. Whenthe distribution function of the particles is not far from Maxwellian, the movement of the plasmacan be described by the two-fluid model [11, 18]. In this section, we will introduce the two-fluidplasma model and the VM model briefly.
The two-fluid model, which captures the separate motion of the electrons and ions, is derivedby taking moments of the plasma kinetic model for each species. This process eliminates themicroscopic velocity space, and at the same time, the macroscopic variables, including thedensity, momentum and energy, are utilized to describe the evolution of the plasma. The detailedform of the two-fluid model is as belowdensity : d n β d t + n β ( ∇ x · u β ) = 0 , momentum : m β n β d u β d t = n β q β [ E + u β × B ] − ∇ x p β − ∇ x · σ β , energy : 32 n β d T β d t = − n β T β ( ∇ x · u β ) − ∇ x · q β − σ β · ∇ x u β , (2.1)where β = i or e , represents electrons and ions respectively. The expression σ · ∇ x u is definedas σ · ∇ x u = X i,j =1 σ ij ∂u j ∂x i . (2.2) m β and q β is the mass and charge of the particles. n β , u β and T β represent the density,macroscopic velocity and temperature. Moreover, p β is the isotropic pressure, σ β is the shearstress tensor, and q β is the heat flux. E and B are the electric and magnetic field respectively,which is decided by the Maxwell equations, the exact form of which is as below ∂ B ∂t + ∇ x × E = 0 , c ∂ E ∂t − ∇ x × B = − ν j , ∇ x · E = ρ c ǫ , ∇ x · B = 0 . (2.3)Here ν and ǫ are the permeability and permittivity of free space [10], and c = ( ν ǫ ) − / is thespeed of light. ρ c and j are the charged density and current density defined by ρ c = X β = i,e q β n β , j ( t, x ) = X β = i,e q β n β u β . (2.4)The two-fluid model (2.1) is not closed, and there are several kinds of closure methods, whichwill lead to different two-fluid models. For example, in the five-moment ideal two-fluid model,3he anisotropic stress σ β and the heat flux q β are set as zero. Combined with the equation ofstate for the ideal gas, we can get the closed five-moment ideal two-fluid model [11]. Similarly,we can also get the ten-moment two-fluid equations [10]. In [4], the closure of σ and q arederived in the kinetic theory of gas with the BGK collision term q = − κ ∇ x T , σ = − µ (cid:18) (cid:2) ∇ x u + ( ∇ x u ) T (cid:3) − I ( ∇ x · u ) (cid:19) , (2.5)where κ is the coefficient of heat conductivity and µ is the coefficient of viscosity. Since onlyone species is considered in the kinetic theory of gas, the footnote index β in (2.5) is eliminated.However, in (2.5) the BGK collision model will lead to wrong Prandtl number, which equals 1for the ideal gas, while the correct one equals 2 /
3. Therefore, other collision models should beutilized to get the correct Prandtl number.Based on the Shakhov collision model, we will try to deduce the two-fluid model with thecorrect Prandtl number based on the kinetic plasma equations. In the next section, the plasmakinetic equations will be introduced briefly.
In the plasma kinetic equations, the plasma is described by the distribution functions f ( t, x , v ),which depends on time, physical space and microscopic velocity space. The plasma kineticequations then have the form below ∂f β ∂t + v · ∇ x f β + q β m β ( E + v × B ) · ∇ v f β = ∂f β ∂t (cid:12)(cid:12)(cid:12) collision , ( x , v ) ∈ Ω × R , (2.6)where β also represents the ions, or electrons respectively, and Ω ∈ R . The electron-magneticfields ( E , B ) are given by the classical Maxwell system (2.3). In the VM model, where for thelow-collisional plasma, the collisions are neglected and collective interactions are assumed todominate the plasma dynamics, which means that ∂f β ∂t (cid:12)(cid:12)(cid:12) collision = 0 . (2.7)In what follows, we are focusing on the derivation of the closed two-fluid model with correctPrandtl number. Therefore, the Shakhov collision model is applied [25], whose exact form is Q Shakhov ( f β ) = ν β ( f sβ − f β ) , f sβ = P β ( t, x , v ) M β ( t, x , v ) , (2.8)with P β = 1 + (1 − Pr)( v − u β ) · q β ( D + 2) ρ β ( t, x ) T β ( t, x ) (cid:18) | v − u β | T β ( t, x ) − ( D + 2) (cid:19) , D = 3 , (2.9)where ν β is the collision frequencies, D is the dimension number of the microscopic velocityspace and M β ( t, x , v ) is Maxwellian, which has the form M β ( t, x , v ) = n β (2 π T β ) / exp (cid:18) − ( v − u β ) T β (cid:19) . (2.10)Moreover, Pr is the Prandtl number which is decided by the type of the particles. u β , T β and q β are the macroscopic variables defined in the last section, whose relationships with the4istribution function f β ( t, x , v ) are as below n β = Z R v f β d v , n β u β = Z R v f β d v , n β T β = 12 Z R | v − u β | f β d v , q β = 12 Z R | v − u | ( v − u β ) f β d v , p β,ij = Z R ( v i − u β,i )( v j − u β,j ) f β d v , i, j = 1 , , . (2.11)The pressure p β is defined as p β = 13 X i =1 p β,ii . (2.12)Since the plasma kinetic equation is seven-dimensional, it is only used to simulate the plasmato capture the essential physics [26]. Some simplified models are introduced to simulate theevolution of the plasma, just by taking moments of the plasma kinetic equations.Different collision models and closure methods will lead to different two-fluid models. In thefollowing sections, we will consider the Shakhov collision model purposely for a correct Prandtlnumber in the reduced model. In this section, we will first derive the hyperbolic moment model for the plasma kinetic equations,which is the base for us to deduce the new two-fluid model. For simplicity, we consider at firstthe single-species case of the non-relativistic electrons under the self-consistent electromagneticfields while the ions are treated as a uniform fixed background. The whole procedure can betrivially extended to the general plasma.Without loss of generality, we consider the dimensionless form of the governing equationsfor one species. Eliminating the footnote index for the distribution function f β ( t, x , v ), thedimensionless equation for the electrons is as below ∂f∂t + v · ∇ x f + ( E + v × B ) · ∇ v f = ν ( f s − f ) , ( x , v ) ∈ Ω × R , (3.1)and the dimensionless form for the classical Maxwell system is ∂ E ∂t − ∇ x × B = − j ,∂ B ∂t + ∇ x × E = 0 , ∇ x · E = ρ − h, ∇ x · B = 0 , (3.2)where the relations of the current density j and charge density ρ with the distribution function f are reduced into j ( t, x ) = Z R f ( t, x , v ) v d v , ρ ( t, x ) = Z R f ( t, x , v ) d v , (3.3)with h the initial background density satisfying Z Ω ( ρ − h ) d x = 0 . (3.4)5he relationship between the macroscopic variables and the distribution function is listed in(2.11). Besides, the Maxwellian is reduced into M ( t, x , v ) = ρ (2 π T ) / exp (cid:18) − ( v − u ) T (cid:19) . (3.5)Following the method in [3], we approximate the distribution function using Hermite basisfunctions as f ( v ) ≈ X | α | M f α H T ,α ( ξ ) , ξ = v − u √T , (3.6)where M is the expansion order, and α = ( α , α , α ) is a three-dimensional multi-index. Thebasis functions H T ,α are defined as H T ,α ( ξ ) = Y d =1 √ π T − αd +12 He α d ( ξ d ) exp (cid:18) − ξ d (cid:19) , (3.7)where He α d is the Hermite polynomial He n ( x ) = ( − n exp (cid:18) x (cid:19) d n d x n exp (cid:18) − x (cid:19) . (3.8)For convenience, He n ( x ) is taken as zero if n <
0, thus H T ,α ( ξ ) is zero when any component of α is negative. Based on the expansion (3.6), we can find that the Maxwellian M in (3.5) equals M ( t, x , v ) = f ( t, x ) H T , ( ξ ) , ξ = v − u √T , (3.9)and it holds for the coefficients f α that f = ρ ( t, x ) , f e i = 0 , X d =1 f e d = 0 , i = 1 , , . (3.10)Moreover, the relationship between heat flux q i , shear stress σ ij and the expansion coefficients f α are σ ij = (1 + δ ij ) f e i + e j , q i = 2 f e i + X d =1 f e d + e i , i, j = 1 , , , (3.11)where the relationship between the shear stress σ ij and the distribution function is σ ij = Z R (cid:18) ( v i − u i )( v j − u j ) − δ ij | v − u | (cid:19) f d v , i, j = 1 , , . (3.12)With the equation of state of the ideal gas, we can derive the relationship between the pressuretensor and the shear stress as p = ρ T , σ ij = p ij − δ ij p. (3.13)Below we briefly introduce the moment systems for (3.1), and refer [6] for the detailedprocedure. With the relationship ∂∂v j H T ,α ( ξ ) = −H T ,α + e j ( ξ ) , (3.14) v j H T ,α ( ξ ) = T H T ,α + e j ( ξ ) + u j H T ,α ( ξ ) + α j H T ,α − e j ( ξ ) , (3.15)6e can get that ∇ v · (cid:2) ( E + v × B ) f (cid:3) = (3.16) − X d =1 E d − X k,m ǫ dkm u k B m f α − e d − X d,k,m =1 ǫ dkm ( α k + 1) B m f α − e d + e k , where the Levi-Civita symbols ǫ dkm are defined as ǫ dkm = , d = k = m cyclic permutation of 1 , , , − , d = k = m anti-cyclic permutation of 1 , , , , ( d − k )( k − m )( m − d ) = 0 . (3.17)The Shakhov collision term can be expanded as ν ( f s − f ) = ν X α ∈ N Q α H T ,α ( ξ ) , ξ = v − u √T , (3.18)with Q α = , | α | < , − Pr5 q j − f α , α = 2 e i + e j , i, j = 1 , , , − f α , otherwise . (3.19)By plugging the expansion (3.6) into (3.1), the general moment equations can be obtained as ∂f α ∂t + X d =1 ∂u d ∂t + X j =1 u j ∂u d ∂x j − E d − X k,m =1 ǫ dkm u k B m f α − e d − X d,k,m =1 ǫ dkm ( α k + 1) B m f α − e d + e k + 12 ∂ T ∂t + X j =1 u j ∂ T ∂x j X d =1 f α − e d + X j,d =1 (cid:20) ∂u d ∂x j (cid:0) T f α − e d − e j + ( α j + 1) f α − e d + e j (1 − H [ | α | − M ]) (cid:1) + 12 ∂ T ∂x j (cid:0) T f α − e d − e j + ( α j + 1) f α − e d + e j (1 − H [ | α | − M ]) (cid:1)(cid:21) + X j =1 (cid:18) T ∂f α − e j ∂x j + u j ∂f α ∂x j + ( α j + 1) ∂f α + e j ∂x j (1 − H [ | α | − M ]) (cid:19) = νQ α , | α | M, (3.20)where M is the expansion order and H [ x ] is defined as H [ x ] = (cid:26) , x = 0 , , otherwise . (3.21)Similar to that in [3], we can deduce the conservation of mass, momentum and the equation for7he temperature as ∂ρ∂t + X j =1 (cid:18) u j ∂ρ∂x j + ρ ∂u j ∂x j (cid:19) = 0 ,∂u d ∂t + X j =1 u j ∂u d ∂x j + 1 ρ X j =1 ∂p jd ∂x j = E d + X k,m ǫ dkm u k B m ,ρ ∂ T ∂t + X j =1 u j ∂ T ∂x j + 23 X j =1 ∂q j ∂x j + X d =1 p jd ∂u d ∂x j ! = 0 . (3.22)Substituting p ij in (3.22) with (3.12), we can derive that ∂ρ∂t + X j =1 (cid:18) u j ∂ρ∂x j + ρ ∂u j ∂x j (cid:19) = 0 ,∂u d ∂t + X j =1 u j ∂u d ∂x j + 1 ρ X j =1 (cid:18) ∂σ jd ∂x j + δ jd ∂p∂x j (cid:19) = E d + X k,m ǫ dkm u k B m ,ρ ∂ T ∂t + X j =1 u j ∂ T ∂x j + 23 X j =1 ∂q j ∂x j + X d =1 ( σ jd + δ jd p ) ∂u d ∂x j ! = 0 . (3.23)Finally, substituting ∂u d ∂t and ∂ T ∂t in (3.20) with (3.23), we can get the quasi-linear system forthe higher order moments of the Shakhov collision term as F α = − νf α , ∀ | α | = 2 ,F α = ν (cid:18) − Pr5 q j − f α (cid:19) , α = 2 e i + e j , i, j = 1 , , ,F α = − νf α , ∀ | α | M and α = 2 e i + e j , i, j = 1 , , , (3.24)where F α = ∂f α ∂t + X j =1 (cid:18) T ∂f α − e j ∂x j + u j ∂f α ∂x j + ( α j + 1) ∂f α + e j ∂x j (1 − H [ | α | − M ]) (cid:19) + X j =1 3 X d =1 T f α − e d − e j + ( α j + 1) f α − e d + e j (1 − H [ | α | − M ]) − p jd ρ X k =1 f α − e k ! ∂u d ∂x j + X j =1 X k =1 (cid:0) T f α − e k − e j + ( α j + 1) f α − e k + e j (1 − H [ | α | − M ]) (cid:1) − T ρ ∂ρ∂x j + 16 ρ X d =1 ∂p dd ∂x j !! − X j =1 3 X d =1 f α − e d ρ ∂p jd ∂x j − ρ X k =1 f α − e k ! X j =1 ∂q j ∂x j − X d,k,m =1 ε dkm ( α k + 1) B m f α − e d + e k , (3.25)Together with (3.23) and (3.24), we derive the hyperbolic moment system for the one speciesplasma kinetic equations (3.1). 8t is known that one can deduce the Euler and Navier-Stokes equations from the Boltzmannequation based on the zeroth- and first-order expansions of the distribution function. For thesimilarity of the Boltzmann and the plasma kinetic equation, we will try to deduce the two-fluidplasma model based on the moment equations. To derive our new two-fluid model based on HME for the plasma kinetic equations, we carry outa Maxwellian iteration to give the closure relations. The Maxwellian iteration was introducedby Ikenberry and Truesdell [13, 30] as a technique to derive the NSF and Burnett equations fromthe moment equations. It is then used to analyze the order for the magnitude of each moment,which is known as the COET (Consistently Ordered Extended Thermodynamics) method.Here, Maxwellian iteration is utilized to derive models with the first order of accuracy. Tobegin with, we introduce the scaling t = t ′ /ǫ and x i = x ′ i /ǫ , and rewrite the moment equations(3.24) with time and spatial variables t ′ and x ′ i . Thus, a factor ǫ − is introduced to the right-sideof (3.24). In this section, we will work on the scaled equations, and the prime symbol t ′ and x ′ will be omitted. With some rearrangement, (3.24) is rewritten as F α = − νǫ f α , ∀ | α | = 2 ,F α = νǫ (cid:18) − Pr5 q j − f α (cid:19) , α = 2 e i + e j , i, j = 1 , , ,F α = − νǫ f α , ∀ | α | M, and α = 2 e i + e j , i, j = 1 , , . (4.1)Assuming that the asymptotic expansions for all the moments have the form f α = f (0) α + ǫf (1) α + ǫ f (2) α + · · · , (4.2)we will begin the iteration based on (4.1). In the original Maxwellian iteration method, werequire that f (0) is the local Maxwellian M ( t, x , v ). The corresponding assumption for thecoefficients is f (0) α = (cid:26) ρ, if | α | = 0 , , otherwise . (4.3)Then noting that f α = 0 if | α | = 1 based on (3.10), the iteration begins from | α | >
2, which canbe written as the iterative scheme below f ( n +1) α = −G α (cid:16) f ( n ) β | β ∈ N (cid:17) , ∀ α ∈ N and | α | > , n = 0 , , , · · · , (4.4)and the macroscopic variables ρ , u and T are treated as ρ = O (1) , u = O (1) , T = O (1) . (4.5)When | α | = 2, (4.1) is changed into f α = − ǫν F α , ∀ | α | = 2 . (4.6)In the iteration, { f α , | α | > } are treated as of high order, and the terms { f α , | α | = 2 } in F α arealso omitted, due to the ǫ on the right hand side. Thus, we can easily deduce the approximation9o f α when | α | = 2 as f (1)2 e i = − ν − X j =1 X d =1 p jd ∂u d ∂x j + ∂q j ∂x j ! + ρ T ∂u i ∂x i , i = 1 , , ,f (1) e i + e k = − ν (cid:20) T ρ (cid:18) ∂u d ∂x j + ∂u j ∂x d (cid:19)(cid:21) , i, k = 1 , , , i = k. (4.7)Similarly, when | α | = 3, (4.1) is changed into f α = − ǫν F α + 1 − Pr5 q j , α = 2 e i + e j , i, j = 1 , , ,f α = − ǫν F α , α = e i + e j + e k , i, j, k = 1 , , , i = j = k. (4.8)From (3.11), we can find that q is at the same order as { f α , | α | = 3 } . Assuming that theasymptotic expansion for heat flux q i is q i = q (0) i + ǫq (1) i + ǫ q (2) i + · · · , i = 1 , , , (4.9)then q (0) i = 0 due to (4.3). Thus it holds for { f (1) α , | α | = 3 } that f (1)2 e i + e k = − ν T − T ∂ρ∂x k + 16 X d =1 ∂p dd ∂x k ! + 1 − Pr5 q (1) k , i, k = 1 , , ,f (1) e i + e j + e k = 0 , i, j, k = 1 , , , i = j = k. (4.10)With the relationship between the pressure tensor and the temperature (2.12) and (3.13), (4.10)is reduced into f (1)2 e i + e k = − ν ρ T ∂ T ∂x k + 1 − Pr5 q (1) k , i, k = 1 , , ,f (1) e i + e j + e k = 0 , i, j, k = 1 , , , i = j = k. (4.11)Moreover, with the same method, we can derive that for the higher order of coefficients, f (1) α = 0 , | α | > . (4.12)For now, the first order of iteration is finished. From (3.11) and the expression of f (1) α , we canget the approximation to the shear stress and heat flux at order O ( ǫ ) σ ii = 2 f e i ≈ − ǫν − X j =1 X d =1 p jd ∂u d ∂x j + ∂q j ∂x j ! + ρ T ∂u i ∂x i , i = 1 , , , (4.13) σ ij = f e i + e j ≈ − ǫν ρ T (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) , i = j, i, j = 1 , , , (4.14) q i = 2 f e i + X d =1 f e d + e i ≈ − ǫν ρ T ∂ T ∂x i , i = 1 , , . (4.15)From (4.14) and (4.15), we can find that σ ij and q i are all at order O ( ǫ ). Substituting (4.14)and (4.15) into (4.13) and omitting the terms at order O ( ǫ ), we can get the simplified form of(4.13) with (3.12) σ ii ≈ − ǫν ρ T ∂u i ∂x i − X j =1 ∂u j ∂x j , i = 1 , , . (4.16)10ntroducing the parameter viscosity conductivity µ and thermal conductivity κ , we can get thefinal closure term by letting ǫ → σ ij = − µ ∂u i ∂x j + ∂u j ∂x i − δ ij X k =1 ∂u k ∂x k ! , q i = − κ ∂ T ∂x i , i = 1 , , , (4.17)where µ = ρ T ν , κ = 1Pr 52 ρ T ν . (4.18) Remark . From the definition of Prandtl numberPr = 52 µκ , (4.19)we can find that (4.18) could give the correct Prandtl number for any type of particles. However,if the BGK collision model is utilized instead of the Shakhov collision model in the deduction,the wrong Prandtl number Pr = 1 will be derived.For the moment, we have got the approximation of the macroscopic variables such as theshear stress σ ij and the heat flux q i with the density ρ , macroscopic velocity u and temperature T . Next, we will deduce the new reduced model with these approximations and the momentequations (3.22) and (3.23). Rewriting (3.23) into the vector form and adopting (4.17), we canget the equation system as ∂ρ∂t + ∇ x · ( ρ u ) = 0 ,ρ d u d t = −∇ x p − ∇ x · σ + ρ ( E + u × B ) , ρ d T d t = [( − p I − σ ) · ∇ x ] · u − ∇ x · q , (4.20)where d · d t = ∂ · ∂t + u · ∇ x · , and I is the 3 × σ = ( σ ij ) × = − µ (cid:20) ∇ x u + ( ∇ x u ) T −
23 ( ∇ x · u ) I (cid:21) ,p = ρ T , q = − κ ∇ x T , µ = ρ T ν , κ = 1Pr 52 ρ T ν . (4.21)Until now, we have deduced the reduced model for (3.1). Compared with the two-fluid modelin Section 2.1, we can find (4.21) has the same form of (2.5), but with different coefficients ofviscosity and thermal conductivity. This is due to the different collision models and closuremethods, where the Shakhov collision model is adopted here and will inherit the correct Prandtlnumber from the plasma kinetic equations.Since the shear stress and the heat flux are all expressed by density, macroscopic velocity andthe temperature, this new reduced model has the same number of degree of freedom as the five-moment two-fluid model, which greatly decreases the computational complexity compared withthe plasma kinetic equations. More importantly, the shear stress and heat flux are not simply setas zero, which is done in the five-moment model [26], but expressed by the macroscopic variablessuch as the density, macroscopic velocity and temperature. Thus, it is expected that this newreduced model could be more capable to describe the problems with anisotropic pressure tenserand large heat flux compared to the five-moment model.11hough this new model is only deduced for the single-species of plasma, it can be extendednaturally to the plasma kinetic equation (2.6). For the collisionless plasma where the collectiveinteractions dominate the plasma dynamics, we can neglect the momentum and heat transferbetween different species. Thus, we can derive the new two-fluid model for multi-species ofplasma with the same Maxwellian iteration. In this case, the new reduced model (4.20) ischanged into density : d n β d t + n β ( ∇ x · u β ) = 0 , momentum : m β n β d u β d t = n β q β [ E + u β × B ] − ∇ x p β − ∇ x · σ β , energy : 32 n β d T β d t = − n β T β ( ∇ x · u β ) − ∇ x · q β − σ β · ∇ x u β , (4.22)with the closure relations σ β = − µ β (cid:20) ∇ x u β + ( ∇ x u β ) T −
23 ( ∇ x · u β ) I (cid:21) ,p = n β T β , q β = − κ β ∇ x T β , µ β = n β T β ν β , κ β = 1Pr 52 n β T β ν β . (4.23) Remark . We only derive the new reduced model for the plasma which does not contain theinteractions between different species. For the general case, the corresponding new reducedmodel can also be derived but with some differences in the closure relations.For the different collision models, the same Maxwellian iteration method can also be utilizedto get the reduced models for the plasma kinetic equations, but the final form of the two-fluidmodel may be different.
In this section, the validation of this new model is studied, especially for the problem wherePrandtl number will bring negligible effect. To show the capability of the new model, thesteady-state and time-dependent problems are tested.
First example
In this example, we focus on the steady-state problem. The Hartmann flowproblem is often used to test the resistive, viscous MHD models near the high collision regime[16, 21]. Below, a similar problem is examined to show that the new two-fluid model can capturecorrect Prandtl number.The problem setting shown in Figure 1 is a plasma set between two infinitely large plates.Two infinite conducting plates are driving shear flow along x -direction with different velocitiesand temperatures. A constant magnetic field is imposed along z -direction, which drives a currentin y -direction. The flow generated along y is interacting with B z , which will suppress the viscousboundary layer. Since we are focusing on the effect of Prandtl number, the dimensionless problemof only one species with a fixed background is tested. To further reduce this problem, we justset the viscosity and thermal conductivity µ and κ in (4.22) as constant. Moreover, the wallsare no-slip, but no adiabatic. 12 zy B z u b u t u x Figure 1: Hartmann problem.For the steady state problem, the governing equations are reduced into X j =1 ∂∂x j ( ρu j ) = 0 , X j =1 u j ∂u d ∂x j + 1 ρ X j =1 (cid:18) ∂σ jd ∂x j + δ jd ∂p∂x j (cid:19) = E d + X k,m ǫ dkm u k B m ,ρ X j =1 u j ∂ T ∂x j + 23 X j =1 ∂q j ∂x j + X d =1 ( σ jd + δ jd p ) ∂u d ∂x j ! = 0 , (5.1)with the closure p = ρ T , σ = − µ (cid:20) ∇ x u + ( ∇ x u ) T −
23 ( ∇ x · u ) I (cid:21) , q = − κ ∇ x T , (5.2)where x , x , x are corresponding to x, y, z . Different boundary conditions will lead to differentsteady state solutions. In this test, the plasma is assumed to be incompressible and the initialdensity is set as ρ = 1. The velocity and temperature of the upper wall are set as u up = 1, and T up = 1, while those of the bottom wall are set as u bottom = 0, and T bottom = 0. The heightbetween the two plates is L = 1, with the imposed magnetic field B z = 1. For this problem isdesigned for a slab geometry, there is no gradient in the x or y direction, which means that ∂ · ∂x = ∂ · ∂y = 0 . (5.3)Then from the continuity equation, we can derive u z = 0 . (5.4)Steady state Faraday’s Law in a slab geometry [21] could give E x = E y = 0 . (5.5)Then the governing system (5.1) is reduced into µ ∂ u x ∂z = − ρu y B z ,µ ∂ u y ∂z = ρu x B z ,κ ∂ T ∂z = − µ (cid:18) ∂u x ∂z (cid:19) + (cid:18) ∂u y ∂z (cid:19) ! . (5.6)13hen we can get the solution to the steady state problem as u x ( z ) = a exp ( bz ) + a exp ( − bz ) + a exp (cid:0) ¯ bz (cid:1) + a exp (cid:0) − ¯ bz (cid:1) ,u y ( z ) = − i (cid:2) a exp ( bz ) + a exp ( − bz ) − a exp (cid:0) ¯ bz (cid:1) − a exp (cid:0) − ¯ bz (cid:1) (cid:3) , T ( z ) = − a a h exp (cid:0) ( b + ¯ b ) z (cid:1) − exp (cid:0) ( b − ¯ b ) z (cid:1) − exp (cid:0) ( − b + ¯ b ) z (cid:1) + exp (cid:0) ( − b − ¯ b ) z (cid:1) i + Cz, (5.7)where a = 12 (exp( b ) − exp( − b )) , a = −
12 (exp( b ) − exp( − b )) ,a = 12 (cid:0) exp(¯ b ) − exp( − ¯ b ) (cid:1) , a = − (cid:0) exp(¯ b ) − exp( − ¯ b ) (cid:1) , b = 1 √ µ (1 + i) ,C = 1 + 4Pr5 a a " exp − √ √ µ ! − exp − √ √ µ ! − exp √ √ µ ! + exp √ √ µ ! . (5.8)In Figure 2, the solutions of the temperature for different Prandtl numbers are plotted wherethe viscous conductivity is set as µ = 0 .
01. We can see that the appearance of the temperaturewith different Prandtl numbers varies greatly, which means Prandtl number is quite importantto get the correct model. The new reduced two-fluid model could inherit the correct Prandtlnumber from the plasma kinetic equations, and could be more capable to describe the behaviorof the plasma.
Pr=0.5Pr=0.6Pr=0.67Pr=0.75Pr=0.85Pr=1
Figure 2: Temperature for different Prandtl numbers. x zy B y T u x u z Figure 3: Time-dependent problem.14 econd example
In this test, we are focusing on the time-dependent problem, where a plasmaset between two infinitely large plates is also studied, but one plate is at z = 0, and the otheris at z = + ∞ . The periodic boundary is utilized in the x -direction. The setting of the problemis shown in Figure 3. A constant magnetic field of strength B y = 1 acts in the direction of the y axis. The density and velocity of the left wall are set as ρ left = 1 and u left = 0. The initialtemperature is changing with x as T left = 1 − . πx ) . (5.9)The initial condition of the particles is set as ρ = 1 , u = 0 , T = 1 . (5.10)The similar problem is studied in [7, 28] to show the ghost effect brought by the kinetic theory.In this problem setting, there is also no gradient in y -direction, and we can derive ∂ · ∂y = 0 , u = 0 . (5.11)The external electric field is added to balance the self-consistent electric field. Thus the system(5.1) is reduced into ∂ρ∂t + ∂ρu ∂x + ∂ρu ∂x = 0 ,∂u ∂t + u ∂u ∂x + u ∂u ∂x − µ ρ ∂ u ∂x − µρ ∂ u ∂x − µ ρ ∂ u ∂x ∂x + ∂ T ∂x + T ρ ∂ρ∂x = − Bu ,∂u ∂t + u ∂u ∂x + u ∂u ∂x − µ ρ ∂ u ∂x − µρ ∂ u ∂x − µ ρ ∂ u ∂x ∂x + ∂ T ∂x + T ρ ∂ρ∂x = Bu ,∂ T ∂t + u ∂ T ∂x + u ∂ T ∂x + X i =1 , " − µ ρ (cid:18) ∂u i ∂x i (cid:19) − κ ρ ∂ T ∂x i + 23 T ∂u i ∂x i − µ ρ (cid:18) ∂u ∂x + ∂u ∂x (cid:19) + 8 µ ρ ∂u ∂x ∂u ∂x = 0 . (5.12)The numerical solutions at time t = 0 .
01 with different Prandtl numbers are studied. Theupwind scheme is adopted to approximate the first-order derivatives with the central differencescheme to approximate the second-order derivatives in (5.12). In the numerical test, the gridsize is N = 400. Due to the similar behavior of the solutions with different Prandtl numbers,Figure 4 shows the behavior of ρ , u , u and T for Prandtl number 2 / µ = 0 .
01. Wecan find u and T are all changing periodically in the x -direction, which is consistent with theproblem setting. To show the affection of Prandtl number, the local numerical solutions withdifferent Prandtl numbers at different x positions are plotted. Figure 5 illustrates the differencesof the local solution ρ , u , u and T at position x = 0 and 0 . a) ρ (b) u (c) u (d) T Figure 4: Numerical solution with Prandtl number Pr = 2 / t = 0 . (a) ρ, x = 0 -4 (b) u , x = 0 (c) u , x = 0 (d) T , x = 0 (e) ρ, x = 0 . -4 (f) u , x = 0 . (g) u , x = 0 . (h) T , x = 0 . Figure 5: Numerical solutions with different Prandtl numbers at t = 0 .
01. The top column is at x = 0, and the bottom column is at x = 0 . We derived a new two-fluid model for the plasma using Maxwell iteration based on the HMEsystem. With the Shakhov collision model, the new two-fluid model takes the correct Prandtlnumber in its closure relations. Though it has the same degree of freedom as the five-momenttwo-fluid model, this new reduced model is expected to have improved performance for problemswith anisotropic pressure tensor and large heat flux, since anisotropic pressure tensor and heatflux are expressed by the density, macroscopic velocity and temperature instead of simply settingas zero. It is our future work to investigate by numerical simulations the performance of thenew model in practical applications.
Acknowledgements
We thank Prof. Yana Di at BNU-HKBU United International College for the useful sugges-tions. The work of Ruo Li is partially supported by the National Science Foundation of China16Grant No. 91630310) and Science Challenge Project (No. TZ2016002). Yixiao Lu is partiallysupported by the elite undergraduate training program of School of Mathematical Sciences inPeking University. Yanli Wang is supported by Science Challenge Project (No. TZ2016002) andNational Natural Science Foundation of China (Grant No. U1930402).
References [1] Z. Cai, Y. Fan, and R. Li. Globally hyperbolic regularization of Grad’s moment system.
Comm. Pure Appl. Math. , 67(3):464–518, 2014.[2] Z. Cai, Y. Fan, and R. Li. A framework on moment model reduction for kinetic equation.
SIAM J. Appl. Math. , 75(5):2001–2023, 2015.[3] Z. Cai, R. Li, and Y. Wang. Solving Vlasov equation using NR xx method. SIAM J. Sci.Comput. , 35(6):A2807–A2831, 2013.[4] D. Callen.
Fundamentals of Plasma Physics . University of Wisconsin, http://homepages.cae.wisc.edu/ callen/book.html , 2006.[5] C. Cheng and J. Johnson. A kinetic-fluid model.
J. Geophys. Res. , 104(A1):413–427, 1999.[6] Y. Di, Z. Kou, and R. Li. High order moment closure for Vlasov-Maxwell equations.
Front.Math. China , 10(5):1087–1100, 2015.[7] F. Filbet and T. Xiong. A hybrid discontinuous Galerkin scheme for multi-scale kineticequations.
J. Comput. Phy. , 372:841 – 863, 2018.[8] J. Friedberg.
Ideal Magnetohydrodynamics . Plenum Press, New York, 1987.[9] S. Gilliam.
A 13-moment two-fluid plasma physics model based on a Pearson type-IV dis-tribution function . Master’s thesis (University of Wisconsin, Madison), 2011.[10] A. Hakim and H. Ammar. Extended MHD modelling with the ten-moment equations.
J.Fusion Energy , 27(1):36–43, Jun 2008.[11] A. Hakim, J. Loverich, and U. Shumlak. A high resolution wave propagation scheme forideal two-fluid plasma equations.
J. Comput. Phys. , 219(1):418 – 442, 2006.[12] J. Huba. Hall magnetohydrodynamics - A tutorial.
Space Plasma Simulation , 615:166–192,2003.[13] E. Ikenberry and C. Truesdell. On the pressures and the flux of energy in a gas accordingto Maxwell’s kinetic theory I.
J. Rat. Mech. Anal. , 5(1):1–54, 1956.[14] N. Iwasawa, A. Ishida, and L. Steinhauer. Tilt mode stability scaling in field-reversedconfigurations with finite Larmor radius effect.
Phys. Plasma , 7(931), 2000.[15] E. Johnson.
Gaussian-moment relaxation closures for verifiable numerical simulation of fastmagnetic reconnection in plasma . Ph.D thesis (University of Wisconsin, Madison), 2011.[16] O. Jones, U. Shumlak, and D. Eberhardt. An implicit scheme ofr nonideal magnetohydro-dynamics.
J. Comput. Phys. , 130:231–242, 1997.1717] N. Krall and A. Trivelpiece.
Principles of Plasma Physics . McGraw-Hill, New York, NY,1973.[18] A. Laguna, N. Ozak, A. Lani, H. Deconinck, and S. Poedts. Fully-implicit finite volumemethod for the ideal two-fluid plasma model.
Comput. Phys. Commun. , 231:31–44, 2018.[19] R. Mason. An electromagnetic field algorithm for 2D implicit plasma simulation.
J. Comput.Phys , 71:429–473, 1987.[20] R. Mason and C. Cranfill. Hybrid two-dimensional electron transport in self-consistentelectromagnetic fields.
IEEE Trans. Plasma Sci. , 14:45–52, 1986.[21] S. Miller and U. Shumlak. A multi-species 13-moment model for moderately collisionalplasmas.
Phys. Plasmas , 23(8):082303, 2016.[22] V. Oraevskii, R. Chodura, and W. Feneberg. Hydrodynamic equations for plasmas in strongmagnetic fields - I: Collisionless approximation.
Plasma Physics , 10(9):819–828, jan 1968.[23] R. Peterkin, M. Frese, and C. Sovinec. Transport of magnetic flux in an arbitary coordinateALE code.
J. Comput. Phys. , 140(1):148–171, 1997.[24] P. Petkaki, M. Freeman, T. Kirk, C. Watt, and R. Horne. Anomalous resistivity and thenonlinear evolution of the ion-acoustic instability.
J. Geophys. Res. , 111:A01205, 2006.[25] E. Shakhov. Generalization of the Krook kinetic relaxation equation.
Fluid Dyn. , 3(5):95–96, 1968.[26] U. Shumlak and J. Loverich. Approximate Riemann solver for the two-fluid plasma model.
J. Comput. Phys. , 187(2):620 – 638, 2003.[27] P. Snyder, G. Hammett, and W. Dorland. Landau fluid models of collisionless magnetohy-drodynamics.
Phys. Plasmas , 4(11):3974–3985, 1997.[28] Y. Sone, K. Aoki, S. Takata, H. Sugimoto, and A. Bobylev. Inappropriateness of theheatconduction equation for description of a temperature field of a stationary gas in thecontinuum limit: Examination by asymptotic analysis and numerical computation of theBoltzmann equation.
Phys. Fluids , 8:628, 1996.[29] M. Tippett. Tokamak transport based on the 13-moment model.
J. Plasma Phys. ,54(1):77104, 1995.[30] C. Truesdell and R. Muncaster.
Fundamentals of Maxwell’s Kinetic Theory of a SimpleMonatomic Gas: Treated as a Branch of Rational Mechanics . Academic Press, 1980.[31] R. Vann, R. Dendy, G. Rowlands, T. Arber, and N. Ambrumenil. Fully nonlinear phe-nomenology of the Berk-Breizman augmentation of the Vlasov-Maxwell system.
Phys.Plasma , 10:632, 2003.[32] L. Wang, A. Hakim, A. Bhattacharjee, and K. Germaschewski. Comparison of multi-fluid moment models with particle-in-cell simulations of collisionless magnetic reconnection.