Numerical properties of high order discrete velocity solutions to the BGK kinetic equation

Numerical properties of high order discrete velocity solutions to the BGK kinetic equation

Applied Numerical Mathematics 61 (2011) 410–427 Contents lists available at ScienceDirect Applied Numerical Mathematics www.elsevier.com/locate/apnu...

2MB Sizes 2 Downloads 18 Views

Applied Numerical Mathematics 61 (2011) 410–427

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Numerical properties of high order discrete velocity solutions to the BGK kinetic equation A.M. Alekseenko Department of Mathematics, California State University, Northridge, CA 91330, USA

a r t i c l e

i n f o

Article history: Received 28 April 2009 Received in revised form 16 June 2010 Accepted 7 November 2010 Available online 16 November 2010 Keywords: Kinetic equations Discontinuous Galerkin methods Transient gas flows Normal shock wave Heat transfer

a b s t r a c t A high order numerical method for the solution of model kinetic equations is proposed. The new method employs discontinuous Galerkin (DG) discretizations in the spatial and velocity variables and Runge–Kutta discretizations in the temporal variable. The method is implemented for the one-dimensional Bhatnagar–Gross–Krook equation. Convergence of the numerical solution and accuracy of the evaluation of macroparameters are studied for different orders of velocity discretization. Synthetic model problems are proposed and implemented to test accuracy of discretizations in the free molecular regime. The method is applied to the solution of the normal shock wave problem and the one-dimensional heat transfer problem. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction High resolution simulations of gas flows in the non-continuum regime are emerging as an important area of practical interest. Applications of non-continuum flows include rocket propulsion, high-altitude flight, vacuum technology and gas flows in micro- and nano-devices. The most accurate description of non-continuum flows is given by the kinetic theory of gases (cf. [6,8]). Although there has been significant progress in understanding the properties of governing equations of non-continuum flow and different gas–gas and gas–surface interaction models (cf. [31,29,39,40,36]), the solvers are generally hampered by large computational costs and numerical uncertainty in the solution. The compelling need for increased efficiency and ability to capture high gradients and shocks in the numerical solution prompted the development of methods based on high order weighted essentially non-oscillatory (WENO) (cf. [7,21]) and discontinuous Galerkin (DG) (cf. [23]) spatial discretizations. However, in the design of such methods, the accuracy of the discretization in the velocity space, usually done by discrete ordinate methods, is often assumed. Because truncation errors introduced by discretizations in the velocity variable may limit the overall accuracy of the methods understanding the properties of the employed discretizations is important. Furthermore, the use of high order discretizations in both the physical and velocity variables can provide a more consistent approach to the design of accurate and efficient numerical solvers. In this work we propose high order discretizations of model Boltzmann equations in both spatial and velocity variables, based on DG methods. The DG discretizations in the velocity variable yield symmetric hyperbolic systems of equations that are solved numerically by Runge–Kutta discontinuous Galerkin (RKDG) techniques. We compare the performance of new high order discrete velocity (HODV) methods to some low order methods and study numerically the methods’ convergence with respect to the discretization in the velocity space. Whether the description of a flow is kinetic or continuum is determined by the flow’s Knudsen number, Kn = λ/Λ. Here λ is the molecular mean free path, i.e., the average distance that a molecule is expected to travel before colliding with

E-mail address: [email protected]. 0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.11.005

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

411

another molecule, and Λ is the characteristic length scale of the physical system. Flows with Kn  1 are continuum flows. For such flows, small volumes of gas contain sufficiently many molecules so that the local velocity distribution is well approximated by the Maxwellian equilibrium distribution (the so-called continuum assumption). Navier–Stokes equations are traditionally used to model flows in the continuum and near-continuum regimes with slip boundary conditions applied to flows with 10−3 < Kn < 10−1 (cf. [31,34]). The continuum assumption, however, does not hold for flows with Kn > 10−1 . For, in this case the frequency of molecular collisions is not sufficiently large to maintain thermodynamic equilibrium through energy exchange. It is now common to use the so-called hybrid approach in which the equations of hydrodynamics are solved in the continuum regions and kinetic equations are solved in the non-continuum regions, cf. [9,20]. Non-continuum  ), which is defined by the gas can be modeled on the microscopic level using the molecular distribution function f (t , x, u  )x  following property: f (t , x, u u gives the number of molecules found at time t in the volume x = x  y  z contain = (u , v , w ). The ing the point x = (x, y , z) whose velocity belongs to the volume d v = u  v  w containing the point u molecular distribution function is a solution to the Boltzmann kinetic equation

∂  x f (t , x, u ) = Q ( f ),  ) + u · ∇ f (t , x, u ∂t

(1)

where the contribution to the change of the distribution function due to molecular collisions is modeled by the non-linear operator Q ( f ). The complete form of the Boltzmann collision operator involves integration over five-dimensional space at each point of six-dimensional phase space, thus making the solution of the Boltzmann equation computationally very challenging. Statistical methods, such as the direct simulation Monte Carlo (DSMC) methods, have emerged as powerful techniques for the simulation of flows with Kn > 101 ([19,27], see also modifications of the DSMC method for transitional flows in [24,3]). However, for flows with 10−1 < Kn < 101 , such as those in microdevices, statistical methods become prohibitively time-consuming due to the high signal-to-noise ratio and to the length of time required to reach a steady state. A deterministic solution to the Boltzmann equation can be obtained using some approximate models, the most widely known of which are the Bhatnagar–Gross–Krook (BGK) [5] and the ellipsoidal-statistical (ES-BGK) [17] model equations. In the BGK model, the collision operator is replaced with a simpler operator

Q ( f ) = ν ( f 0 − f ), where

ν (t , x) is the collision frequency, and f 0 is the Maxwellian equilibrium distribution function  ¯ (t , x))2   −3/2 ( u−u  ) = n(t , x) 2π RT (t , x) . f 0 (t , x, u exp − 2RT (t , x)

(2)

¯ and temperature, T , are defined as follows The gas density, n, bulk velocity, u,



 ) du , f (t , x, u

n(t , x) = R3



¯ (t , x) = n(t , x)u

 f (t , x, u ) du , u R3

n(t , x) T (t , x) =

1



3R

¯ )2 f (t , x, u ) du , ( u−u

(3)

R3

where R is the normal gas constant. The BGK model satisfies conservation of moments and has positive entropy generation rate [28]. It has been shown [2,35] to produce qualitatively correct results for flows for which the distribution functions are close to the Maxwellian equilibrium distribution. The ES-BGK model has much of the simplicity of the BGK model and has the advantage of reproducing the correct Prandtl number. Similar to the BGK model, the ES-BGK satisfies the conservation properties and entropy inequality [4]. Thus, it is known to produce solutions that are more physically accurate than the BGK model (see, however, [26]). Experimenting with the ES-BGK model will be the subject of a following paper. However, we expect that the stability properties of both ES-BGK and BGK models are the same. Thus most of the results obtained in this work for the BGK model are expected to apply to the ES-BGK model as well. One approach to the solution of model kinetic equations is to use the discrete ordinate method, which consists of replacing the integrals in (3) with quadrature formulas and enforcing (1) only for a discrete set of velocities used as abscissas for the quadrature formulas. The method’s relative simplicity and the guaranteed symmetric hyperbolic form of the resulting equations makes it very attractive for numerical implementation. Discrete ordinate methods based on Gauss-type quadrature formulas have been introduced for the BGK model (cf. [18,32]). In [32], a priori error estimates for the discrete ordinate solution were obtained for the linearized plane Couette flow problem. Discrete ordinate methods based on quadrature formulas on uniform velocity grids have been used as well (cf. [23,1]). One obvious disadvantage of the discrete ordinate method is that the conservative relationships that hold for the continuous solution are satisfied by the numerical solution

412

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

only up to truncation errors. Moreover, accurate prediction of such truncation errors can be complicated. A number of schemes have been proposed for the solution of model Boltzmann equations (cf. [25]) that employ small corrections of the discrete equilibrium distribution function so as to satisfy the conservation of mass, momentum and energy in the discrete solution. While the use of conservative techniques usually results in a significant improvement in the simulations, these techniques, however, do not control truncation errors in the numerical solution, but rather project the solution to some suitable physical space. Additionally, the conservative techniques tend to become more time-expensive as resolution increases in the velocity variable. Instead of employing conservative techniques, we plan to achieve a reasonable approximation to the solution of (1) by taking the advantage of high order DG approximations in the velocity space. We study the convergence of the discrete solution with respect to the velocity variable and attempt to quantify the amount of error involved in different discretizations of the BGK equations. We are particularly focused on accuracy estimates for evaluating macroparameters of the numerical solution. For this purpose, we develop a series of experiments in the free molecular regime, i.e., in the case of vanishing collision frequency, using exact solutions to the transport equation. We apply HODV discretizations of the BGK equations to the solution of the normal shock wave problem. Our results have shown that HODV methods are equally or more efficient than discrete ordinate methods. Unless the ordinates are chosen in a such a way as to enable high order quadratures in the velocity space, we expect that the HODV methods are more efficient. In [15], RKDG discretizations in the temporal and spatial variables and spectral Galerkin discretizations in the velocity variable are introduced for the Boltzmann equation. The particular form of basis functions and the quadrature formulas used in [15] result in a high order numerical scheme that is virtually equivalent to the discrete ordinate method for which the ordinates are chosen at the Gauss–Hermite nodes (up to a necessary scaling). The resulting method is fast and simple in implementation. The difference between the HODV methods and the spectral Galerkin method, however, is that DG discretizations have considerably more flexibility in the design of adaptive schemes. The generalization of the method of [15] to DG techniques will be the author’s future work. DG methods emerged as a powerful approach to the solution of the Boltzmann–Poisson equation describing the electron transport, cf. [7,14,10]. In fact, the theoretical study of the DG methods presented in [14] applies to the methods of this paper after some straightforward modifications are made. However, the equations studied in this paper do not contain derivatives with respect to the velocity variable and have a different collision operator. This makes practical implementation of the methods essentially different. The sections below are organized as follows. In Section 2 we introduce the HODV methods. Sections 3 and 4 deal with the discontinuous Galerkin formulations of the HODV models. In Section 5 we summarize the Runge–Kutta integration used in our approach. Section 6 describes the one-dimensional reduction of the BGK equations. Lastly, Sections 7, 8, 9 and 10 are dedicated to numerical results. 2. Discontinuous Galerkin high order velocity approximation Here we develop the discretization of (1) in the velocity variable. For simplicity, we will only be concerned with approximations with finite support. Thus we consider a bounded region Υ ⊂ R3 in the velocity space with the property that the  ) and its integrals are small in R3 \ Υ . In a typical application, f (t , x, u ) decreases exponentially distribution function f (t , x, u if the velocity variable approaches infinity. Therefore such Υ can be found. We note, however, that the method extends to infinite domains by using a mixture of finite and infinite elements with appropriate basis functions. We introduce a partition of Υ into a finite number of polyhedral regions V i , i = 1, . . . , M, and consider basis functions λl,i ( u ), l = 1, . . . , s, defined on V i . On each V i we seek a solution in the form

 )| V i = f (t , x, u

s 

f l,i (t , x)λl,i ( u ).

(4)

l =1

Substituting (4) into (1), multiplying by the basis function λm,i ( u ) and integrating over V i we obtain s  l =1

D ml,i ∂t f l,i (t , x) + ∂x

s 

u T ml x) + ∂ y ,i f l,i (t , 

l =1





f 0 λm,i −

s 



s 

v T ml x) + ∂ z ,i f l,i (t , 

l =1

s 

w T ml x) ,i f l,i (t , 

l =1

D ml,i f l,i (t , x) ,

(5)

l =1

Vi

u v w where the components of the matrices T ml ,i , T ml,i , T ml,i and D ml,i are given by



u T ml ,i =



u λl,i λm,i , Vi

v T ml ,i =



v λl,i λm,i , Vi

w T ml ,i =

w λl,i λm,i , Vi

 D ml,i =

λl,i λm,i .

(6)

Vi

u u v v w w Introducing the vector notation fi = f l,i and the matrix notations T ml ,i = Ti , T ml,i = Ti , T ml,i = Ti and D ml,i = Di we rewrite (5) in compact form using matrix multiplication:

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

413

 f 0 λm,i − Di fi (t , x) .

 Di ∂t fi (t , x) + ∂x Tui fi (t , x) + ∂ y Tiv fi (t , x) + ∂z Tiw fi (t , x) = ν

(7)

Vi

We notice that the matrices Tui , Tiv , Tiw and Di are symmetric. In addition, Di is positive definite. Thus (7) has the form of a symmetric hyperbolic system. Well-posedness and different discretization techniques have been extensively studied for this type of systems. The class of admissible basis functions λl,i ( u ) that can be used in DG velocity discretizations is quite large. Some suitable choices for both structured and non-structured grids are described in [11,12]. In Section 6, we use basis functions that are products of Legendre polynomials on rectangular partitions. The resulting Di are diagonal. A spectral Galerkin method using products of Lagrange polynomials on Gauss–Hermite nodes is presented in [15]. This discretization has a remarkable feature in that all of the matrices Tui , Tiv , Tiw are diagonal. However, this discretization uses the entire velocity space as a single velocity cell. It is anticipated that multi-cell DG approximations will be more efficient in capturing sharp gradients in the solution. (Therefore it will be interesting to see if the method of [15] can be extended to DG approximations.) By multiplying (7) by (Di )−1 we obtain the discrete velocity model

   ∂t fi (t , x) + ∂x Tˆ ui fi (t , x) + ∂ y Tˆ iv fi (t , x) + ∂z Tˆ iw fi (t , x) = ν (Di )−1 f 0 λm,i − fi (t , x) ,

(8)

Vi

where Tˆ iv = (Di )−1 Tiv , Tˆ ui = (Di )−1 Tiv and Tˆ iw = (Di )−1 Tiv . 3. Discontinuous Galerkin weak formulation—the upwind flux We consider a partition of the spatial domain Ω into polyhedral elements K j . As is customary in DG approximation, we assume that the numerical solution fi (t , x) to (8) belongs to a space H of smooth functions within each K j . However, the solution may be discontinuous at the faces of K j . The weak formulation follows by the multiplication of (8) by a test function ϕ (x) ∈ H and integration over K j . After integrating by parts in the second, third and fourth terms, we obtain



 ∂t fi (t , x)ϕ (x) − Kj

 =

Kj





ν (Di )−1

Kj

Tˆ ui fi (t , x)∂x ϕ (x) −



Tˆ iv fi (t , x)∂ y ϕ (x) −

Kj

 f 0 λm,i − fi (t , x) ϕ (x),

 Kj

Tˆ iw fi (t , x)∂z ϕ (x) +



Tˆ nij fi (t , x)ϕ (x)

∂K j

(9)

Vi

where the matrix Tˆ nij is given by

Tˆ nij (x) = nx, j (x)Tˆ ui + n y , j (x)Tˆ iv + n z, j (x)Tˆ iw .

(10)

Here nx, j , n y , j and n z, j are the components of the unit normal vector to ∂ K j pointing outward. Consider the boundary integral in (9). Since the numerical solution can be discontinuous on ∂ K j , the value fi |∂ K j is defined in a non-local sense only. Different approximations for the boundary integral, called the numerical fluxes, are considered cf. [11,12]. This current paper is concerned with the upwind flux which amounts to the substitution of the incoming components of fi with their limits on ∂ K j from the outside and the outgoing components with their limits from the inside.

First, we notice that Tˆ nij is diagonalizable since matrix Di is positive definite. Thus we consider the eigenvalues of Tˆ nij ,

βl , and the corresponding unit eigenvectors ξl , l = 1, . . . , s and introduce the matrix Ai j whose l-th column is ξl . Next we define the diagonal matrices

  Λ− = diag min(0, β1 ), . . . , min(0, βs ) ,   Λ+ = diag max(0, β1 ), . . . , max(0, βs )

(11)

and the projection operators 1 L− = Ai j Λ− A− , ij ij

1 L+ = Ai j Λ+ A− . ij ij

Operators L− and L+ discern modes of the solution fi propagating inside and outside K j respectively. Notice that ij ij

Tˆ nij = L− + L+ . ij ij

(12)

414

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

The upwind flux approximation of the boundary term takes the form





 L− fi (t , x)OUT + L+ fi (t , x)IN ϕ (x),

(13)

∂K j

where fi (t , x)OUT and fi (t , x)IN are the limits of fi (t , x) at the boundary ∂ K j taken from the outside and inside, correspondingly. 4. Discontinuous Galerkin ansatz in the spatial variable On each elements K j , we introduce basis functions

fi (t , x) = f l,i (t , x) =

μ(l) 

ϕ p, j (x), p = 1, . . . , k. We seek coefficients fi (t , x) in the form

f l,i ; p , j (t )ϕ p , j (x).

(14)

p =1

Notice that μ(l)  k, the total number of basis functions ϕ p , j (x) in (14), can be different for different components f l,i (t , x).  ) obtained The reason for introducing μ(l) becomes clear when we consider the combined spectral decomposition of f (t , x, u by a substitution of (14) into (4):

 )| K j × V i = f (t , x, u

μ(l) s  

f l,i ; p , j (t )ϕ p , j (x)λl,i ( u ).

(15)

l =1 p =1

If λl,i ( u ) are polynomials of degree at most s and ϕ p , j (x) are polynomials of degree at most k, then an efficient discretization is obtained by setting μ(l) = min(max_degree − l, k), where max_degree is the maximum degree of polynomials in (15), e.g., max_degree = max(k, s). We substitute representation (14) into (9) and use the discrete flux (13). Using the test function ϕq, j (x), q = 1, . . . , μ(l), we obtain μ(l) 

∂t f l,i ; p , j (t )C pq, j −

p =1

s 

u Tˆ lm ,i

m =1



s 

w Tˆ lm ,i

μ (m) 

m =1



 =

ν

f m,i ; p , j (t )C xpq, j −

p =1 s 

z f m,i ; p , j (t )C pq ,j +

( D lm,i )−1

s 

v Tˆ lm ,i

m =1

p =1

s 

m =1

Kj

μ (m) 

m =1



− Llm ,i j

y

f m,i ; p , j (t )C pq, j

p =1

f m,i ; p , j ∗ (t )C ∂pqK, j ∗ +

p =1



f 0 λm,i − f l,i (t , x)

μ (m) 

μ (m) 

s 

+ Llm ,i j

m =1

μ (m) 

f m,i ; p , j (t )C ∂pqK, j

p =1

ϕq, j ,

(16)

Vi

where j ∗ is the index corresponding to the exterior element K j ∗ that shares the particular part of the boundary with K j , and

 C j = C pq, j =

 Cxj = C xpq, j =

ϕ p, j (x)ϕq, j (x), Kj

y

y

Kj



C j = C pq, j =

ϕ p, j (x)∂z ϕq, j (x), Kj



C∂j K = C ∂pqK, j =



z Czj = C pq ,j =

ϕ p, j (x)∂ y ϕq, j (x), Kj

ϕ p, j (x)∂x ϕq, j (x),

ϕ p, j (x)ϕq, j (x) dσ ,



C∂j ∗K = C ∂pqK, j ∗ =

∂K j

ϕ p, j∗ (x)ϕq, j (x) dσ .

∂K j

A simpler but perhaps less efficient form of (16) is obtained when the index p have the same upper limit,

(17)

μ(l)

p =1



k

p =1 .

μ(l) = k. In this case, all summations with respect to

Introducing the matrix of coefficients

fi j (t ) = f l,i ; p , j (t ), with a slight abuse of notation we rewrite (16) using matrix multiplication

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

y ∂t fi j (t )C j − Tˆ ui fi j (t )Cxj − Tˆ iv fi j (t )C j − Tˆ iw fi j (t )Czj + L− f ∗ (t )C∂j ∗K + L+ f (t )C∂j K ij ij ij ij    s  = ν ( D lm,i )−1 f 0 λm,i − fi (t , x) ϕq, j . Kj

m =1

415

(18)

Vi

Finally, we multiply (18) with the matrix inverse of C j from the left to obtain: y ∂t fi j (t ) − Tˆ ui fi j (t )Cˆ xj − Tˆ iv fi j (t )Cˆ j − Tˆ iw fi j (t )Cˆ zj + L− f ∗ (t )Cˆ ∂j ∗K + L+ f (t )Cˆ ∂j K ij ij ij ij    s  = ν ( D lm,i )−1 f 0 λm,i − fi (t , x) ϕq, j (C j )−1 , Kj

m =1

(19)

Vi

where

Cˆ xj = Cxj (C j )−1 ,

Cˆ j = C j (C j )−1 , y

y

Cˆ zj = Czj (C j )−1 ,

Cˆ ∂j K = Cxj (C j )−1 .

(20)

In the general case of μ(l), formulas (16) need to be used. The equations, however, are somewhat simpler if the basis functions in spatial variables are orthogonal. In this case C j is diagonal and the necessary inverses can be computed immediately for any choice of μ(l). An example of this approach is presented in Section 6 in one spatial dimension using Legendre’s polynomials. 5. Runge–Kutta integration in time The discretization in time is done by Runge–Kutta methods. In the experiments below we used optimal Runge–Kutta methods of second and third orders [13], a fourth order method and a consistent fifth order method [30]. We notice that the results shown in Sections 7 and 8 do not use adaptive time stepping techniques, for which consistent Runge–Kutta methods are very well suited. Adaptive techniques will be implemented in future. Total Variation Diminishing (TVD) methods [37,38] and the Total Variation Bounded (TVB) methods [11] are known to limit growth of spurious oscillatory modes. However, we did not observe a significant difference between solutions obtained by optimal Runge–Kutta methods of orders two and three and the TVD methods. The fact that solutions to the BGK equation do not exhibit strong spurious oscillations may be due to the diffusive nature of the non-linear collision term. 6. One-dimensional BGK equations We apply HODV methods to the solution of gas flow between two parallel infinite plates. The plates are parallel to the Y Z plane. We assume that density, bulk velocity, temperature etc. are constant on each of the plates. As a result, the  ) is constant in y and z variables. Moreover, we assume that the properties of gas molecular distribution function f (t , x, u are isotropic. Since conditions are not changing on the plates, no information is lost if the distribution function is averaged with respect to the velocity components v and w. This suggests to introduce reduced distribution functions

∞ ∞ f 1 (t , x, u ) =

 ) dv dw , f (t , x, u

(21)

 ) dv dw . v 2 f (t , x, u

(22)

−∞ −∞ ∞ ∞

f 2 (t , x, u ) = −∞ −∞

The governing equations for f 1 (t , x, u ) and f 2 (t , x, u ) can be deduced from (1) as:

     −1/2 (u − u¯ (t , x))2 − f 1 (t , x, u ) , ∂t f 1 (t , x, u ) + u ∂x f 1 (t , x, u ) = ν n(t , x) 2π RT (t , x) exp − 2RT (t , x)     1/2  RT (t , x) (u − u¯ (t , x))2 − f 2 (t , x, u ) . ∂t f 2 (t , x, u ) + u ∂x f 2 (t , x, u ) = ν n(t , x) exp − 2π 2RT (t , x)

(23) (24)

where n, u¯ and T are defined by (3) and depend only on t and x. HODV methods are used to discretize (23) and (24). We select a rectangular domain in x and u variables, e.g., [− X , X ] × [−U , U ] and partition it into rectangles X j × U i := [x j −1/2 , x j +1/2 ] × [u i −1/2 , u i +1/2 ], j = 1, . . . , N, i = 1, . . . , M. We introduce one-dimensional basis functions ϕ p , j (x), p = 0, . . . , k, and λl,i (u ), l = 0, . . . , s, on intervals X j and U i , respectively:



ϕ p, j (x) = P p

2( x − x j )

x j

 ,



λl,i (u ) = P l

2( u − u i )

u i

 ,

416

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

where x j = (x j +1/2 + x j −1/2 )/2, x j = x j +1/2 − x j −1/2 , u i = (u i +1/2 + u i −1/2 )/2, u i = u i +1/2 − xi −1/2 and P l (t ) is the l-th Legendre’s polynomial. On each rectangle X j × U i we seek for a solution in the form

f 1,2 (t , x, u ) =

μ(l) s  

1,2

f l,i ; p , j (t )ϕ p , j (x)λl,i (u ).

l =0 p =0

(25)

Here the number of terms in the sum is determined by s, k and the function μ(l) = min(max_deg − l, k). Our results suggest using s  k, s = max_deg. This perhaps is surprising given the fact that the transport part of the Boltzmann equation does not contain derivatives in the velocity variable. One might, indeed, expect using a higher order of interpolation in x rather than in u. However, numerical experiments presented in the next section, indicate that, at least in the free molecular regime, the solution quickly looses smoothness in the velocity variable due to the dependence of the transport operator on the velocity. Unless high order approximation is used in u, the numerical solution rapidly looses accuracy. We substitute representation (25) into (23) and (24) and follow the steps described in Sections 2, 3 and 4 to obtain the HODV discontinuous Galerkin discretization:

∂t ˜f m1 ,i ;q, j (t ) − =

x j

s 

S qp ˜f m1 ,i ; p , j (t ) +



2q + 1

x j









αm,i Φ ˜f m1 ,i; p, j , ˜f m1 ,i; p, j∗ , ϕq, j



−1 1 A ml ,i G ( f ) l ,i ,

l =0

Xj

2q + 1

x j

αm,i

μ (m)  p =0



2q + 1

μ (m)  p =0

ϕq, j (x)

∂t ˜f m2 ,i ;q, j (t ) − =

x j

αm,i



2q + 1

x j

2q + 1

ϕq, j (x)

s  l =0

Xj

S qp ˜f m2 ,i ; p , j (t ) +



2q + 1

x j

αm,i Φ ˜f m2 ,i; p, j , ˜f m2 ,i; p, j∗ , ϕq, j



−1 2 A ml ,i G ( f ) l ,i ,

(26)

for all m = 0, . . . , s, q = 0, . . . , μ(m), i = 1, . . . , M, j = 1, . . . , N. Here j ∗ = j − 1 is numerical flux Φ is defined as

Φ



˜f 1,2 , ˜f 1,2 ∗ , m ,i ; p , j m ,i ; p , j

ϕq, j =

p =1

2p + 1

u i

m,i ; p , j



⎩ ˜f 1,2

(t , x j +1/2 )ϕq, j (x j +1/2 ) − ˜f m1,,2i ; j (t , x j −1/2 )ϕq, j (x j −1/2 ), αm,i < 0.

(27)

(t )ϕ p , j (x). The matrix Aml,i consists of unit eigenvectors ξm,i of

u λl,i (u )λm,i (u ),

(28)

Ui

that correspond to eigenvalues

˜f m,i ;q, j (t ) =

⎧ ⎨ ˜f 1,2 (t , x j +1/2 )ϕq, j (x j +1/2 ) − ˜f 1,2 (t , x j −1/2 )ϕq, j (x j −1/2 ), αm,i > 0, m ,i ; j m , i ; j −1 m , i ; j +1

μ(m) ˜ 1,2 f

1, 2 where ˜f m,i ; j (t , x) =

T ml,i =



αm,i > 0 and j ∗ = j + 1 if αm,i < 0. The

αm,i . Finally ˜f m1,,2i; p, j (t ) are related to f m1,,2i; p, j (t ) through the identities

min(max_deg−q,s)



−1 A ml f l,i ;q, j (t ),

q = 0, . . . , k ,

l =0

f m,i ;q, j (t ) =

s 

A ml ˜f l,i ;q, j (t ),

q = 0, . . . , min(max_deg − m, k).

l =0

7. Numerical experiments Our first group of experiments is concerned with the estimates of approximation errors in the developed HODV method. We begin from verification of the order of convergence in the regime free of molecular collisions, that is, when the right sides of (23) and (24) are replaced with zeros. We notice that high order DG velocity schemes are not very efficient for simulating collisionless flow because their accurate approximation in the velocity space is not used. Instead, various particle methods (cf. [6,16,27,24,33,22] and references therein) may be employed for such simulations. We however will verify the accuracy of HODV approximations of the transport part of kinetic equations and of some boundary conditions. This will be the first step to ensure the validity of HODV simulations of the BGK equation.

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

417

Fig. 1. The exact solution (29) corresponding to the profile (30) with a = 1024, b = 4, r1 = −0.1 and r2 = 1, shown at different times: (a) t = 0; (b) t = 0.4.

We use (26) to reconstruct the exact solution to the transport operator,





f (t , x, u ) = ˜f −ut ± (x + cu ) ,

(29)

where ˜f (ξ ) is a differentiable function of single variable and c is a constant parameter. Since we are aimed at eventually using (26) for the solution of the normal shock wave problem, we test accuracy of (26) on a function which has gradients comparable to that of the shock wave solution. A profile with exponential gradients is given, e.g., by

˜f (ξ ) = a exp



 −b2 , (ξ − r1 )(ξ − r2 )

r1  x  r2 .

(30)

The function is extended by zero outside the interval [r1 , r2 ]. Fig. 1 shows the exact solution (29), (30) at different values of time. Fig. 2 shows the reconstructed solution for different orders of approximating polynomials and Runge–Kutta methods. Fig. 3 shows the relative error for methods of orders one, three and five: the solution converges in agreement with the selected order of the spatial and temporal discretizations (see also Table 1 for the detailed analysis). Tables 1 and 2 summarize the convergence results and the CPU time requirements, respectively, for methods of different order. The time in seconds is measured on an AMD 32-bit 2.3 GHz processor, Windows Vista operational system. Notice that times for shorter calculations are less accurate since they are contaminated by IO disk writing operations. However, it is the longer times that are most important. Tables 1 and 2 demonstrate the superiority of the high order methods versus low order methods. According to Table 1, it will require at least 128 × 128 uniform cells to achieve the 5% accuracy and about 640 × 640 uniform cells to achieve the 1% accuracy by using the first order method (k = s = 0). Table 2 suggests that these calculations will require 5.6 minutes and 10.5 hours, respectively. Also, according to Table 1 the solution obtained by the fifth order method (k = s = 4) is within the 1% accuracy bound for only 4 × 4 uniform cells. The comparison of the CPU time in Table 2 suggests that by replacing the first order method with a fifth order method we achieve a thousands times improvement in speed. 8. Approximation of the initial data and errors in macroparameters When a direct discretization is applied to (23) and (24) (and (1)), the numerical solution satisfies the continuous conservation laws only up to the truncation errors. Problems of gas dynamics, including the normal shock wave problem considered below, are known to be sensitive to violations of the conservative relationships. By enforcing the discrete conservation laws (cf. [25,26]) it is possible to eliminate the non-physical features of the solution. Furthermore, following a delicate argument, one can show that the discrete equilibrium function introduced in [25] gives the minimum to the discrete entropy [28,4,2].

418

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

Fig. 2. The reconstructed solution at t = 0.4. The degrees k and s of interpolating polynomials in variables x and u, and the numbers N and M of the uniform mesh cells for variables x and u, respectively, are as follows: (a) k = s = 0, N = M = 64; (b) k = s = 2, N = M = 16; (c) k = s = 4, N = M = 4. The order of the Runge–Kutta method is k + 1 in each case.

Fig. 3. Relative L 2 -error, f − f d L 2 / f L 2 , for methods of order one, three and five. The numbers of the mesh cells N and M (N = M) for variables x and u, respectively, are 4, 8, 16, 32 and 64. The degrees k and s of interpolating polynomials in variables x and u, respectively are as follows: (a) k = s = 0; (b) k = s = 2; (c) k = s = 4. Table 1 Comparison of the relative L 2 -error. The asterisk (∗) stands for the improvement ratio. N, M

s, k = 0



s, k = 1



s, k = 2



s, k = 3



s, k = 4



4 8 16 32 64

6.4E–1 4.4E–1 3.0E–1 1.8E–1 1.0E–1

1.4 1.5 1.6 1.8

2.6E–1 1.0E–1 2.5E–2 5.7E–3 1.4E–3

2.6 4.1 4.3 4.1

1.1E–1 1.5E–2 2.0E–3 2.5E–4 3.1E–5

7.3 7.8 8.0 7.9

3.0E–2 2.4E–3 1.6E–4 1.0E–5 6.5E–7

12.5 15.2 15.8 15.6

9.4E–3 3.6E–4 1.2E–5 3.9E–7 1.2E–8

26.0 29.8 31.3 31.9

In this work, however, we avoid using conservative techniques and notice that non-linear equations for correction coefficients are becoming in a sense degenerate as accuracy of the approximation increases. We expect that high order velocity discretizations will give us sufficient accuracy inexpensively. Our interest is, therefore, to obtain estimates on the errors in the evaluation of density, bulk velocity and temperature (some of which are sometimes the conserved quantities). We begin from verifying accuracy of approximation of a linear superposition of two Maxwellian equilibrium distributions used as the initial data for the normal shock wave problem. As is common in gas dynamics we will not be comparing the

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

419

Table 2 Comparison of the CPU time in seconds for AMD 32-bit 2.3 GHz processor. The asterisk (∗) stands for the lapsed time ratio for consecutive resolutions. N, M

s, k = 0



s, k = 1



s, k = 2



s, k = 3



s, k = 4



4 8 16 32 64

0.03 0.13 0.78 5.91 44.80

3.9 6.2 7.6 7.6

0.06 0.38 2.73 20.58 91.74

6.1 7.3 7.5 4.5

0.20 1.11 8.35 41.18 233.35

5.5 7.5 4.9 5.7

0.44 3.23 17.32 108.40 732.39

7.4 5.4 6.3 6.8

1.34 7.80 52.99 396.41 3025.64

5.8 6.8 7.5 7.6

Fig. 4. Initial data for a Mach 10 normal shock wave problem: (a) and (b) are plots of the reduced distribution function f 1 (0, x, u ); (c) is the normalized initial density, (n(0, x) − n1 )/(n2 − n1 ); (d) is the bulk velocity, (u¯ (0, x) − u¯ 2 )/(u¯ 2 − u¯ 1 ); and (e) is the temperature, ( T (0, x) − T 1 )/( T 2 − T 1 ).

physical data to its DG approximation. Instead, we will compare macroparameters of density, bulk velocity and temperature evaluated from the data to those evaluated from the DG approximation. Gauss numerical quadratures are used for the evaluation of macroparameters that are exact on the DG solutions up to roundoff errors. Therefore, errors in macroparameters are dominated by errors in the DG approximations. It is important to understand, however, that the most natural norm to use for evaluating errors of DG approximations, the L 2 -norm, can be large while the errors in the macroparameters are small. We consider simulations of a Mach 10 shock wave in Ar gas (the normal gas constant is R = 208.243, the atomic mass is m = 6.6337E–26 kg). Down the stream, the gas has density d1 = 1.0668E–4 kg/m3 , temperature T 1 = 300 K and bulk velocity u¯ 1 = 3226.61 m/s. Up the stream, the gas has density d2 = 4.1428E–4 kg/m3 , temperature T 2 = 9636.938 K and bulk velocity u¯ 2 = 830.852 m/s. Fig. 4 shows the initial data for the reduced distribution function f 1 (t , x, u ). On interval x ∈ [−0.3, 0.03], f 1 (0, x, u ) coincides with the reduced Maxwellian equilibrium distribution corresponding to n1 = d1 /m, u¯ 1 and T 1 . Similarly, on x ∈ [0.07, 0.13], f 1 (0, x, u ) is determined by n2 = d2 /m, u¯ 2 and T 2 . In the transition region, x ∈ (0.03, 0.07), f 1 (0, x, u ) is a linear superposition of the two equilibrium distributions. Figs. 4(c), (d) and (e), show the normalized number density, |n(0, x) − n1 |/|n2 − n1 |, bulk velocity, |u¯ (0, x) − u¯ 2 |/|u¯ 1 − u¯ 2 |, and temperature, | T (0, x) − T 1 |/| T 2 − T 1 |, respectively. The convergence of the relative error in macroparameters is shown in Fig. 5. Here plots (a), (b) and (c) correspond to DG approximations in the velocity variable by fourth degree polynomials, s = 4, and plots (d), (e) and (f) correspond to approximations by piecewise constant functions, s = 0. The latter coincides with a discrete ordinate method in which the ordinates are taken at the center points of the velocity cells and integration with respect to the velocity variable is replaced by a Riemann sum. The spatial approximation is the same in each case and is by forth degree polynomials, k = 4, on N = 32 uniform grid cells. In each case, the error appears to decrease faster than the order of the DG approximation used. A remarkable feature of the Maxwellian equilibrium distribution is that the density, bulk velocity and temperature obtained

420

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

Fig. 5. Relative errors in the evaluation of macroparameters in the initial data for Mach 10 shock wave. Figures (a) and (d) show log10 (|n(x) − nexact (x)|/|nexact (x)|), (b) and (e) log10 (|u¯ (x) − u¯ exact (x)|/|u¯ exact (x)|) and (c) and (f) log10 (| T (x) − T exact (x)|/| T exact (x)|), evaluated for different numbers of uniform velocity cells, M. Plots (a)–(c) are obtained by a fifth and (d)–(f) by a first order HODV methods. The spatial approximation is by a fifth order DG method on 32 uniform cells in each case.

using low order quadrature formulas, e.g., the trapezoidal rule and the uniform Riemann sum, converge faster than any power of the mesh size. This can be observed in Figs. 5(d), (e) and (f), respectively. A study of macroparameters after a hundred of time steps suggests that the fast convergence of the low order velocity schemes does not hold during the evolution. This may be due to the loss of smoothness in the numerical solution due to truncation errors in the spatial and temporal discretizations. In Fig. 6, solutions obtained by a fifth order HODV method (plots (d)–(f)) and a first order method (plots (g)–(i)) are compared to the solution obtained by a ninth order method, s = 8, on M = 64 velocity cells at time t = 10−6 seconds. The time is about a half of a percent of the time necessary for the shock wave to reach a steady state. It can be seen on plots (a)–(c) that macroparameters do not differ much from their initial values in Figs. 4(c)–(e). However, the fast convergence of the low order velocity schemes can no longer be observed in Figs. 6(g) and (h). Instead, the fifth order method plotted in Figs. 6(d)–(f) performs better. Accidental fast convergence of low order quadrature formulas on the initial data coupled to its spatial non-smoothness made the comparison of accuracy difficult in previous experiments. We therefore consider a simpler problem that consists of evaluating the discrete moments of the exact solution to (23) in the free molecular regime,

f (t , x, u ) = 1 − cos





π (−ut + x) .

(31)

Notice that although this solution remains smooth at all times, as time increases, the solution gains steeper gradients in the direction of u (see Figs. 7(a)–(c)). It is anticipated that when the velocity gradients in the solution reach certain value, the discrete integration techniques will fail. Eq. (23) is solved in the domain (x, u ) ∈ [0, 1] × [0, 10] in the free molecular regime using the exact initial and boundary data (31). The same numerical integration method as in the shock wave problem is used to compute the total mass of the solution,

1 10 m(t ) =

f (t , x, u ) du dx, 0

0

at each moment of time and compare it to the known values. The results are presented in Fig. 7. The graphs of relative errors suggest that the fifth order HODV method (Fig. 7(d)) is significantly more accurate than the first order method (Fig. 7(e)) in computing the total mass.

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

421

Fig. 6. Convergence of macroparameters at an early time in a Mach 10 shock wave problem. Plots (a), (b) and (c) show the normalized number density, bulk velocity and temperature, respectively, at time t = 10−6 . Figures (d) and (g) give log10 (|n(x) − nref (x)|/|nref (x)|), (e) and (h) log10 (|u¯ (x) − u¯ ref (x)|/|u¯ ref (x)|) and (f) and (i) log10 (| T (x) − T ref (x)|/| T ref (x)|), evaluated for different numbers of uniform velocity cells, M. Figures (d)–(f) are obtained by the fifth and (g)–(i) by the first order HODV methods. The reference solution is computed by a ninth order HODV method on 64 velocity cells. The spatial approximation is by a fifth order DG method on 32 uniform cells in each case.

9. The normal shock wave problem We apply scheme (26) to the solution of the one-dimensional normal shock wave problem. We consider the cases of Mach 1.2 and Mach 10 shock waves. The parameters for the Mach 1.2 wave are as follows: the gas is nitrogen; down the stream, the gas has density d1 = 1 kg/m3 , temperature T 1 = 300 K and bulk velocity u¯ 1 = 769.668768 m/s; up the stream, the gas has density d2 = 2.28571424 kg/m3 , temperature T 2 = 623.437512 K and bulk velocity u¯ 2 = 336.73009 m/s. The solutions are presented in Figs. 8 and 9. It is known that the location of the shock wave is sensitive to the initial data and therefore to the truncation errors in the initial data. As a result, locations of shock waves in numerical solutions were different for different spatial and velocity resolutions. To compare the solutions we shift them in Fig. 9 horizontally some small distances so as to make the density match up to twelfth digit at the middle of the wave. This explains a sudden drop in error in Figs. 9(a), (d). The convergence is studied in comparison to the reference solution which is computed by the fifth order HODV method on M = 128 uniform velocity cells. A loss of convergence is observed for M = 32 and 64 velocity cells in Figs. 9(a)–(c). However, since errors propagate from the boundary and grow exponentially toward the interior of the domain, we can conclude that the errors are due to the insufficient accuracy of the spatial discretization, which is fifth order on N = 16 uniform cells. (The temporal discretization is fifth order RK scheme with the time step proportional to the step of the spatial discretization.) Indeed, simulations on N = 32 spatial cells shown in Figs. 9(d)–(f) demonstrate an improvement in convergence. The solution shown in Fig. 8 coincides with the solution obtained in [40] by a second order conservative discretization of

422

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

Fig. 7. Relative error in the total mass of f (t , x, u ) = 1 − cos(π (−ut + x)) on (x, u ) ∈ [0, 1] × [0, 10]. Figures (a), (b) and (c) show the solution at times t = 0, t = 0.25 and t = 1, respectively. Figure (d) shows the relative error in evaluating the total mass using fifth order discrete velocity approximation, s = 4, on uniform grid. Similarly, figure (e) shows the error in the first order method, s = 0. The spatial discretization is by a fifth order DG method on 16 cells in each case.

the BGK equation. In contrast, our HODV methods are currently not conservative. However, a better control over the accuracy in the evaluation of macroparameters allowed the HODV methods to produce a physically meaningful solution. Our next example is the simulation of Mach 10 normal shock wave in argon gas. The parameters of the shock wave are as follows: down the stream, the gas has density d1 = 1.0668E–4 kg/m3 , temperature T 1 = 300 K and bulk velocity u¯ 1 = 3226.61 m/s; up the stream, the gas has density d2 = 4.1428E–4 kg/m3 , temperature T 2 = 9636.938 K and bulk velocity u¯ 2 = 830.852 m/s. The reduced distribution function f 1 (t , x, u ) and the macroparameters are shown in Fig. 10. A comparison of the fifth order HODV, s = 4, and the first order, s = 0, methods is presented in Fig. 11. The spatial discretization is by a fifth order DG method on 32 cells in each case. The reference solution is obtained by a ninth order HODV, s = 8, method on 64 velocity cells. The spatial discretization of the reference solution is by a fifth order method on 64 cells. In this figure the numerical solutions are again shifted horizontally so that values of the density coincide up to twelfth digit at the center of the wave. We recall that the first order HODV method is equivalent to a discrete ordinate method in which ordinates are taken at the center points of the cells. It can be seen that the convergence of the fifth order method is better than that of the first order method. In particular, errors on Figs. 11(a)–(c) are dominated by the errors of spatial approximation in the cases of M = 64 and M = 128. This definitely cannot be achieved in the case of s = 0. Notice that in each simulation the temperature converged to roughly one digit less than the density and average velocity. Consider the solutions whose density and velocity converged within 0.5% (and whose temperature converged within 5%). The first such solution for s = 4 is M = 32 and for s = 0 is M = 512. In Table 3, the computational time in seconds are summarized. Here M is the number of the uniform grid cells in the velocity variable. The final time is t = 2.5 × 10−4 s. According to Table 3 high order methods give at least five times improvement in speed. We only show simulations that converged to a physical solution. In other words, s = 4, M = 16 and s = 0, M = 128 are the minimum resolution requirements for simulations to reach a physical steady state. We also conjecture that the use of conservative techniques in either s = 4 or s = 0 case would allow to re-distribute the error that is currently consolidated at the center of the wave thus providing a reasonable approximation to solutions with coarser resolutions.

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

423

Fig. 8. Mach 1.2 shock wave in N2 . Plots (a)–(c) show the reduced distribution function f 1 (t , x, u ). Plots (d), (e) and (f) give the normalized initial density, bulk velocity and temperature, respectively.

Fig. 9. Convergence of macroparameters in a Mach 1.2 shock wave problem. Figures (a) and (d) give log10 (|n(x) − nref (x)|/|nref (x)|), (b) and (e) log10 (|u¯ (x) − u¯ ref (x)|/|u¯ ref (x)|) and (c) and (f) log10 (| T (x) − T ref (x)|/| T ref (x)|), evaluated for different numbers of uniform velocity cells, M. All figures are obtained by a fifth order HODV method in both space and velocity. The reference solution is computed by a fifth order HODV on 128 velocity cells. Figures (a)–(c) correspond to N = 16 uniform spatial cells and (d)–(f) to N = 32 cells, respectively.

424

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

Fig. 10. Mach 10 shock wave in Ar. Plots (a)–(c) show the reduced distribution function f 1 (t , x, u ). Plots (d), (e) and (f) give the normalized initial density, bulk velocity and temperature, respectively.

Fig. 11. Convergence of macroparameters in a Mach 10 shock wave problem. Figures (a) and (d) give log10 (|n(x) − nref (x)|/|nref (x)|), (b) and (e) log10 (|u¯ (x) − u¯ ref (x)|/|u¯ ref (x)|) and (c) and (f) log10 (| T (x) − T ref (x)|/| T ref (x)|), evaluated for different numbers of uniform velocity cells, M. Figures (a)–(c) are obtained by fifth order HODV approximation in the velocity variable, s = 4, and (d)–(f) by a first order method, s = 0. The spatial discretization is by a fifth order HODV method on 32 cells in each case. The spatial and velocity discretizations of the reference solution are by a fifth order DG method on 64 cells and by a ninth order HODV method on 64 cells, respectively.

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

425

Table 3 Comparison of the CPU time for Mach 10 shock wave. s=4

s=0

M

CPU time, sec

M

CPU time, sec

16 32 64 128

2127 4230 8469 17 384

128 256 512 1024

6182 11 954 23 725 46 738

Fig. 12. The heat shock propagation problem. Figures (a), (b) and (c) show the development of density, bulk velocity and temperature, respectively. The graphs are given for time increments of approximately 3.0E–5 s. Figures (d)–(f) give log10 (|m(t ) − mexact |)/mexact where m(t ) is the total mass of gas in the numerical solution. Figure (d) corresponds to simulations for s = 0, (e) to s = 4 and (f) to s = 8.

10. Heat shock between two parallel plates Our last experiment is concerned with the problem of one-dimensional heat transfer between parallel infinite plates. This problem is sensitive to violations of conservation laws and therefore it is a good test for HODV methods. Let argon gas be located between two stationary plates. The distance between the plates is 0.1 m. The velocity space in the numerical method is spanned from −4000 m/s to 4000 m/s. The gas and the plates are initially at 300 K and the initial gas density is d = 1.0E–4 kg/m3 . This corresponds to Knudsen number Kn = 0.01. Fully diffuse reflection is assumed at the surfaces of the plates. At time t = 0.0 s the plate on the right is instantly heated to 1000 K. We consider heat propagation from the moment of instant heating until the moment when the heat perturbation reaches the other plate. On each time step, we evaluate the total mass of the gas between the plates. The results of simulations are presented in Fig. 12. The development of density, velocity and temperature profiles are shown in Figs. 12(a), (b) and (c), respectively. Each line corresponds to an increment in time of approximately 3.0E–5 s. The final time is 3.0E–4 s which is approximately 1/30 of the time necessary to reach the steady state. The lines are chronologically developing from right to left. In Figs. 12(d)–(f) the graphs of relative error in the total mass of the gas are shown. Fig. 12(d) corresponds to simulations by first order HODV method, s = 0, (e) to fifth order method, s = 4, and (f) to ninth order method, s = 8. The spatial discretization is the same in each case and is by a fifth order DG method on 16 cells. It can be seen that very high order discrete velocity methods, s = 8, which are shown on Fig. 12(f) conserve the total mass the best. For only 8 velocity cells the total mass is accurate within the sixth digit. A good conservation is expected in the case of s = 4. At the same time, in the case of s = 0 the mass growth is the worst. Although the growth rate is only slightly higher in this case than in s = 4 and s = 8 cases, the low order convergence of the total mass with resolution in the velocity variable makes s = 0 methods inferior to the other two.

426

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

Given the fact that simulations for (s = 8, M = 8), (s = 4, M = 16) and (s = 0, M = 64) take approximately the same time, the use of very high order methods is strongly preferable in this problem. 11. Conclusions We develop and implement HODV DG methods for the solution of the BGK equation. Simulations of synthetic problems, the normal shock wave problem and the heat transfer problem have shown that HODV methods are more accurate in the evaluation of macroparameters of the numerical solution than low order methods. The ability of HODV methods to achieve high accuracy inexpensively makes the use of conservative techniques unnecessary in one-dimensional problem. However, conservative HODV techniques may be required to provide meaningful solutions in multidimensional cases when mesh resolution is limited by the computer memory and the computational time considerations. It is interesting to further study the fast convergence of low order quadrature formulas for solutions to the BGK equation. A better understanding of fast convergence can potentially lead to the development of improved techniques for integration in velocity space. It will be interesting to perform a detailed comparison of HODV and conservative techniques in the situations when conservation laws are enforced on a coarse velocity mesh. Also, it is interesting to apply HODV methods to the solution of other test problems that can help evaluate different aspects of the discretizations. Examples of such problems include the spatially homogeneous case, heat transfer problem and plane Couette flow problem. The author plans to implement multidimensional HODV methods for BGK and ES-BGK models and use the models for simulations of slow transitional flows. Acknowledgements I want to thank the U.S. Air Force SFFP program for supporting this work and my advisor, Dr. Ingrid Wysong, for her interest in this work and useful discussions. I thank Dr. Sergey Gimelshein for his involvement and consultation in the duration of this work. I thank Natasha Gimelshein for the critical help in the development of the computer code. I most cordially thank Dr. Alina Alexeenko for her help in the many engineering aspects of the problem. Finally, I want to thank Dr. L.L. Foster for her help in the preparation of the final version of this paper. References [1] A. Alexeenko, C. Galitzine, A.M. Alekseenko, High-order discontinuous Galerkin method for Boltzmann model equations, in: 40th AIAA Thermophysics Conference, Seattle, WA, http://www.aiaa.org/agenda.cfm?lumeetingid=1820P&dateget=25-Jun-08#session9731, 2008. [2] P. Andries, K. Aoki, B. Perthame, A consistent BGK-type model for gas mixtures, http://hal.inria.fr/inria-00072389/en/. [3] P. Andries, J.-F. Bourgat, P. le Tallec, B. Perthame, Numerical comparison between the Boltzmann and ES-BGK models for rarefied gases, Comput. Methods Appl. Mech. Engrg. 191 (31) (2002) 3369–3390. [4] P. Andries, P. le Tallec, J.-P. Perlat, B. Perthame, The Gaussian-BGK model of Boltzmann equation with small Prandtl number, Eur. J. Mech. B Fluids 19 (6) (2000) 813–830. [5] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. i. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (3) (1954) 511–525. [6] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flow, Oxford Engineering Science Series, Clarendon, 1994. [7] J.A. Carrillo, I.M. Gamba, A. Majorana, C.-W. Shu, A WENO-solver for the transients of Boltzmann–Poisson system for semiconductor devices: performance and comparisons with Monte Carlo methods, J. Comput. Phys. 184 (2) (2003) 498–525. [8] C. Cercignani, The Boltzmann Equation and its Applications, Applied Mathematical Sciences, vol. 67, Springer-Verlag, New York, 1988. [9] S. Chen, W. E, Y. Liu, C.-W. Shu, A discontinuous Galerkin implementation of a domain decomposition method for kinetic-hydrodynamic coupling multiscale problems in gas dynamics and device simulations, J. Comput. Phys. 225 (2) (2007) 1314–1330, http://www.sciencedirect.com/science/article/ B6WHY-4N0HJKH-3/2/1a286ba61fbc1e882e65ec7a5400dc22. [10] Y. Cheng, I.M. Gamba, A. Majorana, C.-W. Shu, A discontinuous Galerkin solver for Boltzmann–Poisson systems in nano devices, Comput. Methods Appl. Mech. Engrg. 198 (37–40) (2009) 3130–3150, http://www.sciencedirect.com/science/article/B6V29-4WFG5PG-4/2/3b6e2295a1d5ff9c8d8f7881db6221b6. [11] B. Cockburn, Discontinuous Galerkin methods for convection-dominated problems, in: T. Barth, H. Deconink (Eds.), High-order Methods for Computational Physics, in: Lect. Notes Comput. Sci. Eng., vol. 9, Springer-Verlag, Berlin, 1999, pp. 69–224. [12] B. Cockburn, C.-W. Shu, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (3) (2001) 173–261. [13] L.V. Fausett, Applied Numerical Analysis Using MATLAB, Prentice Hall, 1999. [14] Y. Cheng, I.M. Gamba, J. Proft, Positivity-preserving discontinuous Galerkin schemes for linear Vlasov–Boltzmann transport equations, Mathematics of Computation (2010), in press. [15] M.K. Gobbert, T.S. Cale, A Galerkin method for the simulation of the transient 2-D/2-D and 3-D/3-D linear Boltzmann equation, J. Sci. Comput. 30 (2) (2007) 237–273. [16] J.M. Haile, Molecular Dynamics Simulation: Elementary Methods, Wiley, 1997. [17] L.H. Holway, New statistical models for kinetic theory: methods of construction, Phys. Fluids 9 (9) (1966) 1658–1673. [18] B. Huang, D.P. Giddens, The discrete ordinate method for the linearized boundary value problems in kinetic theory of gases, in: C. Brundin (Ed.), Rarefied Gas Dynamics, vol. 1, Academic Press, New York, 1967, pp. 481–504. [19] M.S. Ivanov, S.F. Gimelshein, Computational hypersonic rarefied flows, Annual Rev. Fluid Mech. 30 (1) (1998) 469–505, http://arjournals.annualreviews. org/doi/abs/10.1146/annurev.fluid.30.1.469. [20] E. Josyula, R. Arslanbekov, V. Kolobov, S. Gimelshein, Testing of the unified flow solver (UFS) for nozzle and plume flows, AIAA paper 2007-209. [21] C. Kim, K. Xu, L. Martinelli, A. Jameson, Analysis and implementation of the gas-kinetic BGK scheme for computational gas dynamics, Internat. J. Numer. Methods Fluids 25 (1) (1997) 21–49. [22] S. Li, W.K. Liu, Meshfree and particle methods and their applications, Appl. Mech. Rev. 55 (2002) 1–34. [23] H. Liu, K. Xu, A Runge–Kutta discontinuous Galerkin method for viscous flow equations, J. Comput. Phys. 224 (2) (2007) 1223–1242.

A.M. Alekseenko / Applied Numerical Mathematics 61 (2011) 410–427

427

[24] M.N. Macrossan, ν -DSMC: A fast simulation method for rarefied flow, J. Comput. Phys. 173 (2) (2001) 600–619, http://www.sciencedirect.com/science/ article/B6WHY-45BC24W-1V/2/272ec05258e1fedde19f6ce80a4b3224. [25] L. Mieussens, Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries, J. Comput. Phys. 162 (2000) 429–466. [26] L. Mieussens, H. Struchtrup, Numerical comparison of Bhatnagar–Gross–Krook models with proper Prandtl number, Phys. Fluids 16 (8) (2004) 2797– 2813. [27] E.S. Oran, C.K. Oh, B.Z. Cybyk, Direct simulation Monte Carlo: Recent advances and applications, Annual Rev. Fluid Mech. 30 (1) (1998) 403–441, http://arjournals.annualreviews.org/doi/abs/10.1146/annurev.fluid.30.1.403. [28] B. Perthame, Global existence to the BGK model of Boltzmann equation, J. Differential Equations 82 (1) (1989) 191–205. [29] B. Perthame, Mathematical tools for kinetic equations, Bull. Amer. Math. Soc. (NS) 41 (2) (2004) 205–244 (electronic). [30] A. Ralston, P. Rabinovich, A First Course in Numerical Analysis, McGraw-Hill, New York, 1978. [31] J.M. Reese, M.A. Gallis, D.A. Lockerby, New directions in fluid dynamics: non-equilibrium aerodynamic and microsystem flows, R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci. 361 (1813) (2003) 2967–2988 (mathematics, physics and engineering). [32] G. Richter, On the convergence of the discrete ordinate method in kinetic theory, SIAM J. Appl. Math. 25 (2) (1973) 149–157. [33] D.H. Rothman, S. Zaleski, Lattice-Gas Cellular Automata: Simple Models of Complex Hydrodynamics, Cambridge University Press, Cambridge, UK, 1997. [34] S. Roy, R. Raju, H.F. Chuang, B. Cruden, M. Meyyappan, Modeling gas flow through microchannels and nanopores, J. Appl. Phys. 93 (8) (2003) 4870– 4879. [35] L. Saint-Raymond, From the BGK model to the Navier–Stokes equations, Ann. Sci. École Norm. Sup. 36 (2) (2003) 271–317, http://www.sciencedirect. com/science/article/B6VKH-48HS9DK-5/2/4b7102c9ed9f501112dc9b08b7c9ae3d. [36] N. Selden, C. Ngalande, N. Gimelshein, S. Gimelshein, A. Ketsdever, Origins of radiometric forces on a circular vane with a temperature gradient, J. Fluid Mech. 634 (2009) 419–431, doi:10.1017/S0022112009007976. [37] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [38] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, II, J. Comput. Phys. 83 (1989) 32–78. [39] Y. Sone, Molecular Gas Dynamics, Modeling and Simulation in Science, Engineering and Technology, Birkhäuser Boston Inc., Boston, MA, 2007 (theory, techniques, and applications). [40] D.C. Wadsworth, N.E. Gimelshein, S.F. Gimelshein, I.J. Wysong, Assessment of translational anisotropy in rarefied flows using kinetic approaches, in: American Institute of Physics Conference Series, vol. 1084, 2008, pp. 206–211.