Formulation of a new set of Simplified Conventional Burnett equations for computation of rarefied hypersonic flows

Aerospace Science and Technology 38 (2014) 64–75

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Formulation of a new set of Simplified Conventional Burnett equations for computation of rarefied hypersonic flows Wenwen Zhao a , Weifang Chen a , Ramesh K. Agarwal b,∗,1 a b

Zhejiang University, Hangzhou 310027, China Washington University in St. Louis, Saint Louis 63130, USA

a r t i c l e

i n f o

Article history: Received 3 March 2014 Received in revised form 8 July 2014 Accepted 29 July 2014 Available online 5 August 2014 Keywords: Burnett equations Hypersonic Rarefied gas dynamics Transitional flow

a b s t r a c t For computation of rarefied flows in continuum-transition regime with Knudsen number Kn of O(1), Burnett equations have been proposed about a century ago as a set of extended hydrodynamics equations (EHE) that represent the second-order departure from thermodynamic equilibrium in the Chapman– Enskog expansion of Boltzmann equation; the first order terms in the expansion result in the Navier– Stokes equations. Over the years, a number of variations of original Burnett equations have been proposed in the literature known as the Conventional Burnett equations, the Augmented Burnett equations and the BGK–Burnett equations. In this paper, another simpler set of Burnett equations is proposed by order of magnitude analysis in the limit of high Mach numbers for hypersonic flow applications. These equations, designated as ‘Simplified Conventional Burnett (SCB)’ equations are stable under small perturbations and do not violate the second law of thermodynamics. An implicit numerical solver is developed for the solution of SCB equations. The SCB equations are applied to compute the hypersonic flow past 2D and 3D blunt bodies for Kn in continuum and continuum-transition regime. The SCB solutions are compared with the Navier–Stokes and DSMC solutions. It is shown that the SCB equations can be employed to compute the hypersonic flow past bodies in continuum-transition regime with much less computational effort because of their simplicity compared to Conventional and Augmented Burnett equations. © 2014 Elsevier Masson SAS. All rights reserved.

1. Introduction In high altitude hypersonic flows past space vehicles, especially in low earth orbit, part of the flow field is in continuum regime and part of it is in rarefied continuum-transition regime where Knudsen number Kn ∼ O(1). This is also the case for gaseous flow in many micro-electro-mechanical (MEM) devices where the low density and small length scale result in Kn ∼ O(1). In continuumtransition regime, Navier–Stokes (NS) equations are not accurate because of continuum breakdown due to rarefaction in flow regions near the bow shock, in Knudsen layer and in the wake of the vehicle [8]. Direct Simulation Monte Carlo (DSMC) [5,6] has been widely used for calculation of rarefied flows; however it requires a large number of particles which make the simulations computationally very costly. As an alternative, higher-order extended hydrodynamics equations (EHE) beyond Navier–Stokes have been proposed for computation of flows in continuum-transition regime.

*

Corresponding author. Tel.: +1 314 935 6091. E-mail address: [email protected] (R.K. Agarwal). 1 William Palm Professor in Department of Mechanical Engineering & Materials Science. http://dx.doi.org/10.1016/j.ast.2014.07.014 1270-9638/© 2014 Elsevier Masson SAS. All rights reserved.

These equations are derived from the classical Boltzmann equation by applying the Chapman–Enskog expansion; they were first derived by Burnett and therefore known as the Burnett equations [9]. They represent the second-order departure from thermodynamic equilibrium; the first-order departure leads to the Navier–Stokes equations. These equations include higher-order stress tensor and heat flux terms in the constitutive relations. It is instructive to recall that the constitutive relations for Navier–Stokes equations are linear Stokes law for the stress tensor and Fourier’s law for the heat flux. Over the years, a number of variants to the original Burnett equations [9] have been proposed; these are summarized in Refs. [1,13] and are known as the ‘Conventional Burnett equations (CBE)’, ‘Augmented Burnett (AB) equations’ and the ‘BGK–Burnett equations’. One of the difficulties in using Burnett equations has been that they have been found to be unstable under small wavelength perturbations; this was concluded by Bobylev [7] by conducting the linearized stability analysis of 1D Burnett equations. This problem of stability was also noted by Fiscko and Chapman [12] and Fiscko [11] who found that the Burnett equations became unstable when the mesh was refined. For a long time, the linearized instability problem impeded the use of Burnett equations

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

65

Nomenclature a r k Ki Kn L L ref Ma p p∞ Pr qi R Re T T∞ T ref

sound velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s nose radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m non-dimensional circular frequency coefficients of stress terms in Burnett equations Knudsen number characteristic length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m reference length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Mach number pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . N/m2 free stream pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . N/m2 Prandtl number heat flux terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . J/(m2 s) gas constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . J/(kg K) Reynolds number temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K free stream temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K reference temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K

for computation of rarefied flows. In 1993, Zhong et al. [25] derived the third-order super Burnett equations and then introduced part of third-order terms into the conventional Burnett equations to stabilize them using the Bobylev’s method; they called this new set of equations the Augmented Burnett (AB) equations. However, the addition of third-order terms added to the complexity and in some cases it has been difficult to obtain convergent stable solutions [1]. Comeaux et al. [10] and Jin and Slemrod [14] studied the compatibility of the conventional Burnett equations with the second law of thermodynamics and concluded that the addition of some third-order terms can result in violation of this law. Balakrishnan et al. [3] and Balakrishnan and Agarwal [2] derived a new set of Burnett equations by using the Bhatnagar–Gross–Krook (BGK) model for the collision integral in the classical Boltzmann equation; they called this new set as the BGK–Burnett equations. Balakrishnan and Agarwal [2] showed that the BGK–Burnett equations are stable under small wave length disturbances and satisfy the Boltzmann H-theorem (they do not violate the second law of thermodynamics). It is also important to note the observation of Welder et al. [18] that it is not sufficient to carry out the linearized stability analysis by ignoring the high-order non-linear terms in Burnett equations to ensure stable solutions especially at larger Knudsen numbers; it does not guarantee the stability of non-linear equations on fine computational grids. Because of their complexity, there has been limited development of 3D computational codes in generalized coordinate system for the solution of Burnett equations for rarefied flows past complex geometries. Majority of the papers consider computation of one-dimensional shock structure using Burnett equations. Among 2D and 3D applications, Zhong et al. [22,24] studied the hypersonic flow past a 2D cylindrical leading edge and axisymmetric flow past blunt bodies using the AB equations while Yun [20] investigated more complex 3D cases in his Ph.D. dissertation using the AB equations. Yun and Agarwal [21] also compared the rarefied hypersonic flow simulations using the BGK Burnett equations and the Augmented Burnett equations. Recently, Bao et al. [4] have performed the 2D linearized stability analysis of different variants of Burnett equations using the Bobylev’s method and have shed additional light on the stability issue in higher dimensions. In this paper we derive a new set of Burnett equations which are simpler than the conventional Burnett equations by performing an order of magnitude analysis of higher-order stress and heat-flux terms and neglecting the terms which become negligibly small in the limit of large Mach number. We call this new set as the

Tw Ts U ui us

wall temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K jump-temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K free stream velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s velocity tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s slip-velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3 free stream density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3 reference density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3 stress tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . N/m2 viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (N·s)/m2 reference viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . (N s)/m2 thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . J/(m s K) coefficients of heat flux terms in Burnett equations mean free path of gas molecules . . . . . . . . . . . . . . . . . . . . m accommodation coefficient of momentum accommodation coefficient of temperature

ρ ρ∞ ρref τi j μ μref κ θi λ

σu σT

Simplified Conventional Burnett (SCB) equations, which is unconditionally stable according to one-dimensional linearized stability analysis. These equations are then solved by an implicit LU symmetric Gauss–Seidel scheme. The Developed SCB code is applied to calculate 2D and 3D blunt body flows and the results are compared with NS solutions and Augmented Burnett solutions as well the DSMC solutions. For hypersonic flow past a 2D cylinder, it is shown that the SCB equations are computationally 33% more efficient than the AB equations and 25% less efficient than the NS equations. 2. Derivation of a new set of Simplified Conventional Burnett (SCB) equations The original Burnett equations were derived from the classical Boltzmann equation for monoatomic gases employing the Chapman–Enskog expansion to the velocity distribution function to second order in Knudsen number, that is to O( K n2 ) [9]. The conventional Burnett equations [1] were derived by replacing all the material derivatives in the original Burnett equations by the expressions for the material derivatives in the Euler equations. In the conventional Burnett equations, the second-order terms for the stress tensor and heat flux can be written as:

τi(j2) = K 1

μ2 ∂ uk ∂ u i p ∂ xk ∂ x j

+ K2 + K3

μ2 p



−

∂ 1 ∂p ∂ uk ∂ u j ∂ u i ∂ uk − −2 ∂ xi ρ ∂ x j ∂ xi ∂ xk ∂ xk ∂ x j



μ2 ∂ 2 T μ2 ∂ p ∂ T + K4 2 2 ρ T ∂ xi ∂ x j ρ R T ∂ xi ∂ x j

μ2 ∂ T ∂ T μ2 ∂ u i ∂ uk + K 6 p ∂ xk ∂ x j ρ T 2 ∂ xi ∂ x j    2 2 μ ∂ uk ∂ T μ 2 ∂ ∂ uk ∂ uk ∂ T (2 ) qi = θ1 +2 + θ2 T ρ T ∂ xk ∂ xi ρ T 3 ∂ xi ∂ xk ∂ xi ∂ xk   2 2 2 μ ∂p μ ∂ μ ∂ T ∂ uk + θ4 + θ5 + θ3 ρ p ∂ xk ρ ∂ xk ρ T ∂ xk ∂ xi + K5

(1)

(2)

where a bar over the derivative in Eqs. (1) and (2) represents a non-divergent symmetrical tensor defined below:

fij =

1 2

1

( f i j + f ji ) − δi j f kk 3

(3)

66

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75



Table 1 Burnett coefficients for a hard sphere and Maxwell gas.

K1 K2 K3 K4 K5 K6

θ1 θ2 θ3 θ4 θ5

Hard sphere gas

Maxwell gas

4.056 2.028 2.418 0.681 1.428 7.424 11.644 −5.822 −3.090 2.418 25.158

3.333 2 3 0 3 8 9.375 −5.625 −3 3 29.25

In Eqs. (1) and (2) the coefficients K i , θi for a hard sphere and Maxwell gas are given in Table 1. In derivation of Simplified Conventional Burnett (SBC) equations, we examine the magnitude of each term in Eqs. (1) and (2) by comparing it with the first-order Navier–Stokes stress term and heat-flux term in the energy equation respectively. In terms of order of magnitude, the first-order viscous stress tensor and heat flux can be expressed as:

τi(j1) ≈ μ (1 )

qi

≈

U L

a2 μ L

aρ λ U

≈

L a3 ρ λ

≈

L

= ρ a2 KnMa

(4)

= ρ a3 Kn

(5)



K1

μ ∂ uk ∂ u i 2

K2

τi(j1) ≈

p ∂ xk ∂ x j

μ

−

p

∂ 1 ∂p ∂ xi ρ ∂ x j



μU pL

μ

≈ KnMa

τi(j1) ≈

2

μ ∂u j ∂ T ρ T ∂ x j ∂ xi 

 (1 )

qi



≈ 

λU μU ≈ = KnMa aL a2 ρ L

∂u j ∂u j ∂ T μ2 2 ∂ +2 θ2 T ρ T 3 ∂ xi ∂xj ∂ xi ∂ x j 2

θ3

μ ∂ p ∂u j ρ p ∂ x j ∂ xi



(1 )

qi

τxx =

μ2 p



α¯ 1 u 2x + α¯ 2 u 2y + α¯ 3 u 2z + α¯ 4 v 2x + α¯ 5 v 2y + α¯ 6 v 2z

+ α¯ 7 w 2x + α¯ 8 w 2y + α¯ 9 w 2z + α¯ 10 u x v y + α¯ 11 v y w z + α¯ 12 w z u x

≈

μU a2 ρ L



= KnMa

μU (1 ) qi ≈ 2 = KnMa a ρL

+ α¯ 16 R T xx + α¯ 17 R T y y + α¯ 18 R T zz + α¯ 19

RT

+ α¯ 22

RT

+ α¯ 25

ρ ρ2 R T R

ρxx + α¯ 20

RT

ρx2 + α¯ 23

RT

¯ 26 T x2 + α

R T

ρ ρ2

RT

ρ y y + α¯ 21 ρ 2y + α¯ 24

¯ 27 T 2y + α

R T

ρ

RT

ρ2

ρzz

ρz2

T z2

R



R

(8)

+ α¯ 6 w 2x + α¯ 7 u 2y + α¯ 8 u 2z + α¯ 9 u 2x

−

+ α¯ 10 v y w z + α¯ 11 w z u x + α¯ 12 u x v y + α¯ 13 v z w y + α¯ 14 w x u z + α¯ 15 u y v x + α¯ 16 R T y y + α¯ 17 R T zz + α¯ 18 R T xx

(6)

The order of magnitude ratio of the second-order heat flux terms in Eq. (2) to first order heat flux term given by Eq. (6) can be expressed as:

θ1

(2 )

p



2

It can be noticed from Eq. (6) that out of seven terms, three terms are of O(KnMa) and four terms are of O( K n/ Ma). Tsien [16] was the first to investigate the ratio of various terms in Burnett stress tensor to the NS stress tensor and pointed out that the terms of O(Kn/Ma) can be neglected at high Mach numbers. There are four such terms in Conventional Burnett equations stress tensor and they have derivatives of temperature in their expressions. These terms become negligible compared to the NS stress term and other high-order stress terms at high Mach number encountered in hypersonic flows. The second-order Burnett stress terms can be simplified using this approximation and thus the threedimensional stress components of SCB equations are obtained as given below:

+ α¯ 28 T x ρx + α¯ 29 T y ρ y + α¯ 30 T z ρz ρ ρ ρ 2 μ τ y(2y) = α¯ 1 v 2y + α¯ 2 v 2z + α¯ 3 v 2x + α¯ 4 w 2y + α¯ 5 w 2z

μ aρ λ ≈ = Kn/Ma ρ LU ρ LU

∂ uk ∂ u j ∂ u i ∂ uk μU −2 τi(j1) ≈ ≈ KnMa p ∂ xi ∂ xk ∂ xk ∂ x j pL  μ2 ∂ 2 T μ τi(j1) ≈ ≈ Kn/Ma K3 ρ T ∂ xi ∂ x j ρ LU  μ2 ∂ p ∂ T μ τi(j1) ≈ ≈ Kn/Ma K4 2 ρ LU ρ R T 2 ∂ xi ∂ x j  μ2 ∂ T ∂ T μ τi(j1) ≈ ≈ Kn/Ma K5 ρ LU ρ T 2 ∂ xi ∂ x j  μ2 ∂ u i ∂ uk μU K6 τi(j1) ≈ ≈ KnMa p ∂ xk ∂ x j pL K2

(7)

+ α¯ 13 u y v x + α¯ 14 v z w y + α¯ 15 w x u z

In Eqs. (4) and (5), U and L represent the free stream velocity and a characteristic length respectively and a is the free stream sound velocity. The order of magnitude ratio of the second-order stress terms in Eq. (1) to first order stress term given by Eq. (4) can be expressed as: 2

μ2 ∂ ∂ u j μU (1 ) qi ≈ 2 = KnMa ρ ∂ x j ∂ xi a ρL  μ2 ∂ T ∂ u j μU (1 ) θ5 qi ≈ 2 = KnMa ρ T ∂ x j ∂ xi a ρL θ4

+ α¯ 19

RT

+ α¯ 22

RT

+ α¯ 25

ρ ρ R T R

2

ρ y y + α¯ 20 ρ 2y + α¯ 23

¯ 26 T 2y + α

R T

RT

ρ

RT

ρ

2

RT

ρzz + α¯ 21

ρz2 + α¯ 24

¯ 27 T z2 + α

R T

ρ

RT

ρ2

ρxx

ρx2

T x2

 R R ¯ ¯ ¯ + α28 T y ρ y + α29 T z ρz + α30 T x ρx ρ ρ ρ 2 μ τzz(2) = α¯ 1 w 2z + α¯ 2 w 2x + α¯ 3 w 2y + α¯ 4 u 2z + α¯ 5 u 2x p

+ α¯ 6 u 2y + α¯ 7 v 2z + α¯ 8 v 2x + α¯ 9 v 2y + α¯ 10 w z u x + α¯ 11 u x v y + α¯ 12 v y w z + α¯ 13 w x u z + α¯ 14 u y v x + α¯ 15 v z w y + α¯ 16 R T zz + α¯ 17 R T xx + α¯ 18 R T y y + α¯ 19

RT

ρ

ρzz + α¯ 20

RT

ρ

ρxx + α¯ 21

RT

ρ

ρyy

(9)

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

+ α¯ 22

RT

ρ

ρ 2 + α¯ 23 2 z

+ α¯ 25

R

R

T R

+ α¯ 28

ρ

¯ 26 T z2 + α

RT

ρ

ρ 2 + α¯ 24 2 x R

¯ 27 T x2 + α

T

¯ 29 T z ρz + α

R

ρ

T



RT

ρ

ρ2 2 y

(2 )

qy =

T 2y

¯ 30 T x ρx + α

R

ρ

μ

 T yρy

(2 )

μ2

(2 )

τxz = τzx =

R

RT

T

ρ2

ρxy + β¯13 T x T y + β¯14

R

R

ρ

ρ

+ β¯15 ρx T y + β¯16

 T xρy

(2 )

qz = (11)

β¯1 w z w x + β¯2 u z u x + β¯3 v z v x + β¯4 w z u z

1

1

ρ

ρ

ρ

1

1

ρ

(15)

ρ

μ2 1 1 1 γ¯1 T z w z + γ¯2 T z u x + γ¯3 T z v y ρ T T T 1

1

T 1

T 1

T

T

+ γ¯4 T x u z + γ¯5 T x w x

+ β¯5 w x u x + β¯6 w y u y + β¯7 w x v y + β¯8 u z v y + β¯9 w y v x

+ γ¯13 ρz w z + γ¯14 ρz u x + γ¯15 ρz v y

RT

ρ

ρzx + β¯13

R T

R

R

ρ

ρ

+ β¯15 ρz T x + β¯16

1

T z T x + β¯14



RT

ρ2

ρ z ρx

p

T z ρx

(12)

R

RT

ρ2

R

R

ρ

ρ

+ β¯15 ρ y T z + β¯16

 T y ρz

ρ y ρz (13)

Examining Eq. (7), it is clear that no term can be discarded at high Mach number since they are all of O(KnMa). The second-order heat-flux terms for SCB equations can therefore be written as:



μ2 1 1 1 γ¯1 T x u x + γ¯2 T x v y + γ¯3 T x w z ρ T T T 1

1

T 1

T 1

T

T

+ γ¯4 T y v x + γ¯5 T y u y + γ¯6 T z w x + γ¯7 T z u z + γ¯8 u xx + γ¯9 u y y + γ¯10 u zz + γ¯11 v xy + γ¯12 w xz 1

1

1

ρ

ρ

ρ

+ γ¯13 ρx u x + γ¯14 ρx v y + γ¯15 ρx w z 1

1

ρ

ρ

1

1

ρ

ρ

+ γ¯16 ρ y v x + γ¯17 ρ y u y + γ¯18 ρz w x + γ¯19 ρz u z

ρ

ρ

1

1

ρ

ρ

14 9 1

 (16)

α¯ i , β¯i , γ¯i are defined by

 (14)

K2 +

2 9

K6,

1

1

3

12

α¯ 2 = K 2 +

K6

2 1 α¯ 4 = − K 2 + K 6 , K6, 12 3 12 1 7 1 1 1 α¯ 5 = − K 1 + K 2 − K 6 , α¯ 6 = K 2 − K 6 3 9 9 3 6 2 1 1 1 α¯ 7 = − K 2 + K 6 , α¯ 8 = K 2 − K 6 , 3 12 3 6 1 7 1 1 2 2 α¯ 9 = − K 1 + K 2 − K 6 , α¯ 10 = K 1 + K 2 − K 6 3 9 9 3 9 9 2 4 4 1 2 2 α¯ 11 = − K 1 − K 2 + K 6 , α¯ 12 = K 1 + K 2 − K 6 3 9 9 3 9 9 2 1 4 1 α¯ 13 = − K 2 + K 6 , α¯ 14 = K 2 − K 6 3 6 3 3 2 1 2 2 α¯ 15 = − K 2 + K 6 , α¯ 16 = − K 2 + K 3 3 6 3 3 1 1 1 1 α¯ 17 = K 2 − K 3 , α¯ 18 = K 2 − K 3 3 3 3 3 2 1 1 α¯ 19 = − K 2 , α¯ 20 = K 2 , α¯ 21 = K 2 3 3 3 2 1 1 α¯ 22 = K 2 , α¯ 23 = − K 2 , α¯ 24 = − K 2 3 3 3 2 2 1 1 α¯ 25 = K 4 + K 5 , α¯ 26 = − K 4 − K 5 3 3 3 3 1 1 2 2 α¯ 27 = − K 4 − K 5 , α¯ 28 = − K 2 + K 4 3 3 3 3 1 1 1 1 α¯ 29 = K 2 − K 4 , α¯ 30 = K 2 − K 4 3 3 3 3 3

T

ρ yz + β¯13 T y T z + β¯14

ρ

1

3 1

α¯ 3 = K 2 +

+ β¯10 w x u y + β¯11 R T yz

ρ

ρ

1

2

α¯ 1 = K 1 −

+ β¯7 v z u x + β¯8 w y u x + β¯9 v x u z RT

ρ

In Eqs. (8)–(16), the coefficients Eq. (17).

β¯1 v y v z + β¯2 w y w z + β¯3 u y u z + β¯4 v y w y

+ β¯12

1

+ γ¯18 ρ y v z + γ¯19 ρ y w y

2

μ

1

+ γ¯16 ρx u z + γ¯17 ρx w x

+ β¯5 v z w z + β¯6 v x w x

(2 )

1

+ γ¯10 w y y + γ¯11 u xz + γ¯12 v yz

+ β¯12

qx =

T

+ γ¯6 T y v z + γ¯7 T y w y + γ¯8 w zz + γ¯9 w xx

+ β¯10 u y v z + β¯11 R T zx

(2 ) (2 ) τ yz = τzy =

T



ρx ρ y



p

T 1

+ γ¯16 ρz w y + γ¯17 ρz v z ρ ρ  1 1 + γ¯18 ρx u y + γ¯19 ρx v x

+ β¯10 v z w x + β¯11 R T xy

ρ

T 1

+ γ¯13 ρ y v y + γ¯14 ρ y w z + γ¯15 ρ y u x

+ β¯7 u y w z + β¯8 v x w z + β¯9 u z w y

+ β¯12

1

+ γ¯10 v xx + γ¯11 w yz + γ¯12 u xy

+ β¯5 u y v y + β¯6 u z v z

RT

1

+ γ¯6 T x u y + γ¯7 T x v x + γ¯8 v y y + γ¯9 v zz

(10)

β¯1 u x u y + β¯2 v x v y + β¯3 w x w y + β¯4 u x v x

p

μ2 1 1 1 γ¯1 T y v y + γ¯2 T y w z + γ¯3 T y u x ρ T T T + γ¯4 T z w y + γ¯5 T z v z

2

(2 ) (2 ) τxy = τ yx =

67

68

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

β¯1 =

1 2

K1 −

β¯3 = − K 2 + β¯5 = β¯7 =

1 2 1 2

K1 − K1 +

β¯9 = − K 2 +

5

K2 +

3 1

1 6

β¯4 =

K6,

4 2

K2 +

3 1

K2 −

3 1

K6,

4

β¯2 =

K6,

1 6 1 3

1 2

β¯6 =

K6,

1

1

2

2

3 1 4 1

K1 − K2 +

5 3 1 6

K2 +

1

K6

6

K6

K6

1 1 K1 + K2 − K6 2 3 3 1 β¯10 = − K 2 + K 6 4

β¯8 =

K6,

1

1

2

2

β¯15 = − K 2 +

β¯16 = − K 2 +

2 2

β¯12 = − K 2 ,

β¯11 = − K 2 + K 3 , β¯14 = K 2 ,

K1 −

1

β¯13 = K 4 + K 5

K4

K4

8

2

2

3 2

3 1

3 1

3 2

3 1

3 1

3

3

3

γ¯1 = θ1 + θ2 + θ3 + θ5 γ¯2 = θ1 + θ2 − θ3 − θ5

1

1

2 1

2 1

2 1

2

1

1

2 2

2

γ¯5 = θ3 + θ5 ,

γ¯6 = 2θ2 + θ3 + θ5

γ¯7 = θ3 + θ5 ,

γ¯8 = θ2 + θ4 ,

2 3

3

2

1

3

6

γ¯11 = θ2 + θ4 ,

γ¯10 = θ4 , 2 2

1

γ¯13 = θ3 γ¯14 = − θ3 , 3 1

3

1

2

1

2

2

2

Fig. 1. The trajectory curves (β vs. α ) for a Maxwellian gas and variation of attenuation coefficient α with wave frequency k for one-dimensional CB equations.

1

γ¯9 = θ4

where

2

2

1

3

6 1

L ref =

γ¯12 = θ2 + θ4 1

γ¯15 = − θ3 ,

γ¯18 = θ3 ,

γ¯17 = θ3 ,

1

γ¯4 = 2θ2 + θ3 + θ5

γ¯3 = θ1 + θ2 − θ3 − θ5 ,

3

γ¯16 = θ3

1 2

(17)

It should be emphasized that the SCB equations do not include the third-order terms from super Burnett equations which is the case with AB equations and are even simpler than the original CB equations. As a result, SCB equations are computationally more efficient than the AB equations; the original CB equations of course cannot be used for calculations since they suffer from the Bobylev instability.

Bobylev [7] has analyzed the linearized stability of the Burnett equations and found that the solution of conventional Burnett (CB) equations was unstable under small wavelength perturbations. In this section, the linearized stability analysis of one-dimensional CB and SCB equations is conducted using Bobylev’s method. For small wavelength perturbation, the flow variables ρ , u and T in the linearized 1D CB and SCB equations can be assumed to be of the form:

ρ = ρ1 e(λt +ikx) u = u 1 e (λt +ikx) (18)

where λ is the amplitude of perturbation and k is the nondimensional circular frequency defined by the wavelength and the mean free path (MFP) as:

k=

2π L ref L

= 4.92

λ L

= 4.92Kn

ρref R T ref

λ=

,

16μref 5ρref



(20)

2π R T ref

λρ1 + iku 1 = 0 4

2

3

3

λu 1 + ikρ1 + ikT 1 + k2 u 1 + i K 2 k3 ρ1 2

2

+ i K 2k3 T 1 − i K 3k3 ρ1 = 0 3

3

2

5

4

4

3

2

9

9

λ T 1 + iku 1 + k2 T 1 − i θ2k3 u 1 − i θ4k3 u 1 = 0

(21)

In order to guarantee the existence of non-trivial solution, the determinant of Eq. (21) is equated to zero and the characteristic equation for λ is obtained as:

3. Linearized stability analysis of 1D CB and SCB equations

T = T 1 e (λt +ikx)



By substituting Eq. (18) in the non-dimensional linearized 1D conventional Burnett equations, one obtains:

2

γ¯19 = θ3

μref

(19)









18λ3 + 69k2 λ2 + λk2 30 + 97k2 − 14k4 + 15k4 3 + 4k2 = 0 (22) The solutions of the characteristic equation (22) is assumed to be of the form λ = f (k) with λ = α + i β . The trajectory curves (β vs. α ) and plots showing the variation of attenuation coefficient α with wave frequency k are shown in Fig. 1. It can be seen from Fig. 1 that the conventional Burnett equations are not unconditionally stable since α > 0 when the attenuation coefficient k is greater than 2.455. The same analysis procedure is applied to SCB equations and the characteristic equation is obtained as:





18λ3 + 69λ2 k2 + λ 28k6 + 121k4 + 30k2 + 60k6 + 45k4 = 0 (23) where all Burnett coefficients are considered for a Maxwellian gas as in the case of Eq. (22)

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

69

Fig. 3. Typical 40 × 80 (grid B) around the cylindrical leading edge.

5. Results and discussion 5.1. Grid convergence study Fig. 2. The trajectory curves (β vs. α ) for a Maxwellian gas and variation of attenuation coefficient α with wave frequency k for one-dimensional SCB equations.

It can be seen from Fig. 2 that the simplified conventional Burnett (SCB) equations are unconditionally stable for all Mach numbers and Knudsen numbers since α < 0 for all values of circular frequency k. 4. Numerical method for SCB equations The implicit Lower-Upper Symmetric Gauss–Seidel method [19] is employed to solve the SCB equations. The time-dependent method allows the Simplified Conventional Burnett (SCB) equations to relax to a steady state from arbitrary initial conditions. The inviscid terms in the governing equations are solved by employing the flux difference splitting scheme and the viscous terms are discretized by the second-order accurate central difference scheme. The no-slip or first-order Maxwell/Smoluchowski slip boundary conditions for velocity and temperature at a wall [15] are respectively employed if the flow is in the continuum or slip regime (the slightly rarefied flow with Kn ∼ 0.01). However, there is need to develop the higher-order boundary conditions for the numerical solution of the Burnett equations. But to date, there is no accurate physical description and definitive formulation of the boundary conditions for the Burnett equations. Nevertheless the application of the Maxwell/Smoluchowski slip boundary conditions can result in reasonably accurate results if the non-equilibrium effects near the solid wall can be considered negligible [23]. In this paper, the first-order Maxwell/Smoluchowski slip boundary conditions for velocity and temperature given by Eqs. (24) and (25) are employed in all the calculations presented in Section 5.

us =

2 − σu

σu

Ts − T w =



∂u λ ∂y



3

y =0

15 2 − σ T 8

σT

+ λ



λ



4

∂T ∂y



2R

πT



∂T ∂x

 (24) y =0

(25) y =0

In Eqs. (24) and (25), σu and σ T denote the accommodation coefficient of momentum and temperature respectively; they are assumed to be σu = 1, σ T = 1.

For 2D blunt body, three grids with 20 × 30 (Grid A), 40 × 80 (Grid B), 80 × 120 (Grid C) mesh are employed to determine a suitable mesh for obtaining the grid independent solution. For 2D blunt body, a cylindrical leading edge is considered with geometric parameters and free-stream conditions given in Ref. [20]. The following geometric and flow conditions are used:

r = 0.02 m

Ma∞ = 10

Kn∞ = 0.1

Re ∞ = 167.9

p ∞ = 2.3881 N/m2 T w = 1000.0 K Pr = 0.72

T ∞ = 208.4 K

γ = 1.4

R = 287.04 m2





s2 K

(26)

In Eq. (26), r is the leading edge radius of the blunt body and R is the gas constant for air. Sutherland’s equation is used to calculate the laminar viscosity:



μ T = μref T ref

 1 .5 

T ref + T su



T + T su

(27)

where T su = 124 K, μref = 1.7161 × 10−5 Pa s, T ref = 273.16 K. Fig. 3 shows a typical 40 × 80 mesh (Grid B) around the cylindrical leading edge. Figs. 4 and 5 show the density and temperature distribution along the stagnation streamline using grids A, B, and C. It is clear that the computations using grid B and C are quite close. Therefore, grid B is considered adequate for maintaining both accuracy and computational efficiency. 5.2. Validation of 2D and 3D SCB codes 2D cylindrical leading edge is investigated in this section to validate the 2D Simplified Conventional Burnett (SCB) equations code. This is the same case studied in Section 5.1 for grid convergence study. Several researchers have already investigated this case by using different variants of Burnett equations as well as the NS equations [20,17]. Here, we compare the present computational results using NS equations and the SCB equations with those in the literature [20] obtained with NS and Augmented Burnett (AB) equations. The non-dimensional density, temperature and velocity distributions along the stagnation streamline of the cylindrical

70

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

Fig. 4. Density distribution along the stagnation streamline.

Fig. 6. Density distribution along the stagnation streamline of the 2D cylindrical leading edge.

Fig. 5. Temperature distribution along the stagnation streamline.

leading edge are shown in Figs. 6, 7 and 8 respectively. For density distribution, there is very little difference between the NS and Burnett solutions as expected. However, for the temperature and velocity distributions, Burnett equations predict a thicker shock. The most important result of this investigation is that the SCB equations predict results very close to those predicted by the AB equations. Fig. 9 shows the comparison of Mach contours between NS and SCB equations. Fig. 10 shows the comparison of Mach contours between SCB and AB equations. Fig. 11 shows the comparison of temperature contours between NS and SCB equations. Fig. 12 shows the comparison of temperature contours between SCB and AB equations. From Figs. 9–12, it can be observed that both SCB and AB equations produce identical results. For validation of 3D SCB code, the hypersonic flow of argon past a 3D hemispherical nose is computed. This case of 3D hemispherical nose with a radius r = 0.02 m has already been computed in the literature using the Augmented Burnett (AB) equations and the DSMC method [22]. The geometric parameters and free stream flow conditions for this case are as follows:

r = 0.02 m

Ma∞ = 10

Kn∞ = 0.1

p ∞ = 2.3881 N/m2

T ∞ = 208.4 K

T w = 208.4 K

Fig. 7. Temperature distributions along the stagnation streamline of 2D blunt body.

Fig. 8. Velocity distributions along the stagnation streamline of 2D blunt body.

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

Fig. 9. Comparison of Mach contours between NS and SCB equations for 2D blunt body.

71

Fig. 12. Comparison of temperature contours between SCB and AB equations for 2D blunt body.

Fig. 10. Comparison of Mach contours between SCB and AB equations for 2D blunt body. Fig. 13. Density distribution along the stagnation streamline of a hemispherical nose.

Figs. 13 and 14 show the density and temperature distributions along the stagnation streamline. These figures show that the results from SCB equations yield the density and temperature distributions in much closer agreement with the DSMC results of Vogenitz and Takara [17] than the results from NS equations and even the Augmented Burnett (AB) equations. The results of the Augmented Burnett equations were computed by Zhong and Furumoto [23] by using an axisymmetric AB solver. Thus, the SCB equations predict the shock wave in a monoatomic gas closer to that predicted by the DSMC method than that predicted by the NS and even AB equations for both the 2D and 3D computed cases. 5.3. Effect of variation in Knudsen number and Mach number on 3D flow of Argon past a hemispherical nose Fig. 11. Comparison of temperature contours between NS and SCB equations for 2D blunt body.

γ = 1.67

Pr = 0.67

R = 208.13 m2



s2 K



(28)

The viscosity coefficient is calculated using the power law:



μ = μref

T

0.5

(29)

T ref

where μref = 2.2695 × 10−5 kg/s m, T ref

= 300.0 K.

In this section, the effect of variation in Knudsen number and Mach number on the 3D flow of Argon past a hemispherical nose is studied using both the NS and SCB equations. Three Knudsen numbers (Kn = 0.01, 0.1 and 0.2) and three Mach numbers (Ma = 5, 10 and 15) are considered. First the flow fields at Ma = 10 for Kn = 0.01, 0.1 and 0.2 are computed. The following geometric and far field flow conditions are considered:

r = 0.02 m T ∞ = 208.4 K

Ma∞ = 10 T w = 208.4 K

72

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

Fig. 14. Temperature distribution along the stagnation streamline of a hemispherical nose.

Fig. 17. Density distribution along the stagnation streamline of a hemispherical nose; Ma = 10, Kn = 0.1.

Fig. 18. Temperature distribution along the stagnation streamline of a hemispherical nose; Ma = 10, Kn = 0.1.

Fig. 15. Density distribution along the stagnation streamline of a hemispherical nose; Ma = 10, Kn = 0.01.

Fig. 16. Temperature distribution along the stagnation streamline of a hemispherical nose; Ma = 10, Kn = 0.01.

γ = 1.67

Pr = 0.67

R = 208.13 m2





s2 K

(30)

Figs. 15 and 16 respectively show the density and temperature distribution along the stagnation streamline at Kn = 0.01. At this low value of the Knudsen number, the flow is in continuum regime and therefore as expected, the NS and SCB computations are in close agreement. Figs. 17 and 18 respectively show the density and temperature distribution along the stagnation streamline at Kn = 0.1. At this value of the Knudsen number, the flow is in continuum-transition regime and therefore there is significant difference in the results between the NS and SCB computations especially for the temperature distribution. Figs. 19 and 20 respectively show the density and temperature distribution along the stagnation streamline at Kn = 0.2. At this value of the Knudsen number, the flow is again in the continuum-transition regime and therefore there is again significant difference in the results between the NS and SCB computations especially for the temperature distribution. These results clearly demonstrate the inability of NS equations to compute the rarefied flow in continuum-transition regime; extended hydrodynamics equations (EHD) such as the SCB equations are needed when Kn ∼ O(1). Next the flow fields at Kn = 0.1 for Ma = 5 and 15 are computed. The following geometric and far field flow conditions are considered:

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

Fig. 19. Density distribution along the stagnation streamline of a hemispherical nose; Ma = 10, Kn = 0.2.

73

Fig. 22. Temperature distribution along the stagnation streamline of a hemispherical nose; Ma = 5, Kn = 0.1.

Fig. 23. Density distribution along the stagnation streamline of a hemispherical nose; Ma = 15, Kn = 0.1.

Fig. 20. Temperature distribution along the stagnation streamline of a hemispherical nose; Ma = 10, Kn = 0.2.

r = 0.02 m

p ∞ = 2.3881 N/m2 T w = 208.4 K Pr = 0.67

Fig. 21. Density distribution along the stagnation streamline of a hemispherical nose; Ma = 5, Kn = 0.1.

Kn∞ = 0.1 T ∞ = 208.4 K

γ = 1.67 R = 208.13 m2





s2 K

(31)

Figs. 21 and 22 respectively show the density and temperature distribution along the stagnation streamline at Ma = 5. At Knudsen number Kn = 0.1, the flow is in continuum-transition regime and therefore there is significant difference in the results between the NS and SCB computations especially for the temperature distribution. As expected the density distribution results with NS and SCB equations are quite close. Figs. 23 and 24 respectively show the density and temperature distribution along the stagnation streamline at Ma = 15. At Knudsen number Kn = 0.1, the flow is in continuum-transition regime and therefore there is significant difference in the results between the NS and SCB computations especially for the temperature distribution. However, the difference between NS and SCB computations increases as the Mach number increases. These results clearly demonstrate the inability of NS equations to compute the hypersonic rarefied flow in continuumtransition regime; extended hydrodynamics equations (EHD) such as the SCB equations are needed when Kn ∼ O(1).

74

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

however are in closer agreement with the DSMC solutions. Based on the results presented in this paper, it is concluded that SCB equations are accurate and more efficient for computing hypersonic rarefied flows in continuum-transition regime compared to other set of Burnett equations namely the AB equations. Conflict of interest statement There is no conflict of interest with any person or entity concerning the results presented in this paper. Acknowledgements

Fig. 24. Temperature distribution along the stagnation streamline of a hemispherical nose; Ma = 15, Kn = 0.1.

The first author of this paper Wenwen Zhao was supported by the China Scholarship Council (No. 201306320107) and the research is also supported by the National Basic Research Program of China (2014CB340201). The authors would like to acknowledge the resources of Computational Fluid Dynamics laboratory at Washington University in St. Louis.

6. Computational efficiency of SCB equations References In the computations, a Dell computer with Intel Core i73720QM [email protected] with 16 GB memory was used. The computation time was 5.0313 × 10−6 , 6.6563 × 10−6 , and 8.8438 × 10−6 s per grid point per time step for NS, SCB and AB equations respectively. For the 2D cylinder case presented in this paper, on a 40 × 80 mesh, the computational time becomes 1.61, 2.13 and 2.83 s for 100 time steps for NS, SCB and AB equations respectively. Thus, SCB equations are computationally 33% more efficient than the AB equations and 25% less efficient than the NS equations. In addition to efficiency, it is important to note that the AB code is much more complex because of the addition of super-Burnett terms compared to the SCB code. Furthermore SCB equations are even simpler than the Conventional Burnett (CB) equations (which are unstable under small perturbations). The small differences among the results of NS, SCB, AB and DSMC are expected because difference in the shear stress and heat flux modeling in these approximations. These differences become smaller when Kn tends to zero in the continuum regime and becomes greater as the Kn increases and becomes of O(1) in the continuum-transition regime. 7. Conclusions A new set of Simplified Conventional Burnett (SCB) equations has been developed for calculating the hypersonic flow past two- and three-dimensional blunt bodies in continuum-transition regime with Kn ∼ O(1). The linearized stability analysis of 1D SCB equations shows that the equations are stable under small disturbances and do not violate the second law of thermodynamics. The solutions of the SCB equations are in close agreement with those obtained using the Augmented Burnett (AB) equations for a 2D cylindrical leading edge of a blunt body and for a 3D hemispherical nose. However, AB equations have many more terms for the stress tensor compared to SCB equations and therefore are relatively less efficient computationally. The SCB and AB solutions are also compared with NS and DSMC solutions. For small Knudsen number in the continuum regime, Kn ∼ O(0.01), the NS solutions and Burnett solutions are in close agreement for all flow variables as expected. However for Kn ∼ O(0.1) or greater, the NS solutions and Burnett solutions are in good agreement for density field but are significantly different for the temperature and velocity fields; the Burnett solutions

[1] R.K. Agarwal, K.Y. Yun, R. Balakrishnan, Beyond Navier–Stokes: Burnett equations for flows in the continuum-transition regime, Phys. Fluids 13 (2001) 3061–3085. [2] R. Balakrishnan, R.K. Agarwal, Numerical simulation of Bhatnagar–Gross– Krook–Burnett equations for hypersonic flows, J. Thermophys. Heat Transf. 11 (1997) 391–399. [3] R. Balakrishnan, R.K. Agarwal, K.Y. Yun, BGK–Burnett equations for flows in the continuum-transition regime, J. Thermophys. Heat Transf. 13 (1999) 397–410. [4] F.B. Bao, Z.H. Zhu, J.Z. Lin, Linearized stability analysis of two-dimension Burnett equations, Appl. Math. Model. 36 (2012) 1902–1909. [5] G.A. Bird, Molecular Gas Dynamics, Clarendon Press, Oxford, 1976. [6] G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press, Oxford, 1994. [7] A.V. Bobylev, The Chapman–Enskog and Grad methods for solving the Boltzmann equation, Sov. Phys. Dokl. 27 (1982) 29–31. [8] I.D. Boyd, G. Chen, G.V. Candler, Predicting failure of the continuum fluid equations in transitional hypersonic flows, Phys. Fluids 7 (1995) 210–219. [9] D. Burnett, The distribution of molecular velocities and the mean motion in a non uniform gas, Proc. Lond. Math. Soc. 40 (1936) 382–435. [10] K.A. Comeaux, D.R. Chapman, R.W. Maccormack, An analysis of the Burnett equations based on the second law of thermodynamics, in: 33rd Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, 1995, AIAA Paper 95-0415. [11] K.A. Fiscko, Study of continuum higher order closure models evaluated by a statistical theory of shock structure, PhD thesis, Stanford University, 1988, 124 pp. [12] K.A. Fiscko, D.R. Chapman, Comparison of Burnett, super-Burnett and Monte Carlo solutions for hypersonic shock structure, in: Intl. Symposium on Rarefied Gas Dynamics: Theoretical and Computational Techniques, Pasadena, CA, 1989, pp. 374–395. [13] L.S. García-Colín, R.M. Velasco, F.J. Uribe, Beyond the Navier–Stokes equations: Burnett hydrodynamics, Phys. Rep. 465 (2008) 149–189. [14] S. Jin, M. Slemrod, Regularization of the Burnett equations via relaxation, J. Stat. Phys. 103 (2001) 1009–1033. [15] J.C. Maxwell, On stresses in rarified gases arising from inequalities of temperature, Philos. Trans. R. Soc. Lond. 170 (1879) 231–256. [16] H.S. Tsien, Superaerodynamics, mechanics of rarefied gases, J. Aeronaut. Sci. 13 (1946) 653–664. [17] F.W. Vogenitz, G.Y. Takara, Monte Carlo study of blunt body hypersonic viscous shock layers, in: Rarefied Gas Dynamics, Pisa, Italy, 1971, pp. 911–918. [18] W.T. Welder, D.R. Chapman, R.W. Maccormack, Evaluation of various forms of the Burnett equations, in: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, American Institute of Aeronautics and Astronautics, 1993, AIAA Paper 93-3094. [19] S. Yoon, A. Jameson, Lower-upper symmetric-Gauss–Seidel method for the Euler and Navier–Stokes equations, AIAA J. 26 (1988) 1025–1026. [20] K.Y. Yun, R.K. Agarwal, Numerical simulation of 3-D augmented Burnett equations for hypersonic flow in continuum-transition regime, Technical report, Wichita State University, Dept. of Aerospace Engineering, 1999.

W. Zhao et al. / Aerospace Science and Technology 38 (2014) 64–75

[21] K.Y. Yun, R.K. Agarwal, R. Balakrishnan, Augmented Burnett and Bhatnagar– Gross–Krook–Burnett equations for hypersonic flow, J. Thermophys. Heat Transf. 12 (1998) 328–335. [22] X.L. Zhong, G.H. Furumoto, Solutions of the Burnett equations for axisymmetric hypersonic flow past spherical blunt bodies, in: 6th Joint Thermophysics and Heat Transfer Conference, American Institute of Aeronautics and Astronautics, 1994, AIAA Paper 94-1959. [23] X.L. Zhong, G.H. Furumoto, Augmented Burnett-equation solutions over ax-

75

isymmetrical blunt bodies in hypersonic flow, J. Spacecr. Rockets 32 (1995) 588–595. [24] X.L. Zhong, R.W. Maccormack, D.R. Chapman, Stabilization of the Burnett equations and application to high-altitude hypersonic flows, in: 29th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, 1991, AIAA Paper 91-0770. [25] X.L. Zhong, R.W. Maccormack, D.R. Chapman, Stabilization of the Burnett equations and application to hypersonic flows, AIAA J. 31 (1993) 1036–1043.