## Abstract

Computations of rarefied gas dynamical flows governed by the semiclassical Boltzmann ellipsoidal-statistical (ES) kinetic model equation using an accurate numerical method are presented. The semiclassical ES model was derived through the maximum entropy principle and conserves not only the mass, momentum and energy, but also contains additional higher order moments that differ from the standard quantum distributions. A different decoding procedure to obtain the necessary parameters for determining the ES distribution is also devised. The numerical method in phase space combines the discrete-ordinate method in momentum space and the high-resolution shock capturing method in physical space. Numerical solutions of two-dimensional Riemann problems for two configurations covering various degrees of rarefaction are presented and various contours of the quantities unique to this new model are illustrated. When the relaxation time becomes very small, the main flow features a display similar to that of ideal quantum gas dynamics, and the present solutions are found to be consistent with existing calculations for classical gas. The effect of a parameter that permits an adjustable Prandtl number in the flow is also studied.

## 1. Introduction

In recent years, owing to the rapid advance of nanotechnology, the device or structure characteristic length scale becomes comparable to the mean free path and the wavelength of energy and information carriers (mainly electrons, photons, phonons and molecules), and the use of classical continuum and kinetic descriptions to treat the transport phenomena may no longer be adequate. A semiclassical Boltzmann equation, which is analogous to the original Boltzmann equation and generalizes the collision term to treat collision of particles of quantum statistics, has been developed by Uehling & Uhlenbeck [1] (see [2]). It is commonly believed that this semiclassical Boltzmann equation is adequate to treat semiclassical transport phenomena in the mesoscale range. Hydrodynamic behaviour of quantum gases has been the subject of some prominent researches [3–5], and applications of quantum Boltzmann hydrodynamic equations have been implemented both in the analysis of electron flows in quantum semiconductor devices [6–8] and in the phonon energy transfer in nanocomposites [9,10]. Furthermore, as the different types of carriers may be involved simultaneously in a single problem, it is desirable to have a method that can allow them to be treated them in a unified and parallel manner. With the semiclassical Boltzmann equation, it is possible to tackle the non-equilibrium transport of particles of arbitrary statistics covering a wide range of flow parameters such as Reynolds number, Mach number and Knudsen number.

The major difficulty encountered in solving the semiclassical Boltzmann equation resembles that encountered in the classical counterpart and arises from the complicated nature of the scattering or collision term. The relaxation time approximation proposed by Bhatnagar *et al.* [11] for the classical Boltzmann equation provides a simplified and useful model of collision terms and retains certain essential features of particle collisions and enables more tractable solution methods. The concept of the Bhatnagar–Gross–Krook (BGK) relaxation time approximation, which conforms to the conservation laws of mass, momentum and energy of two elastically colliding particles and eventually reaches an equilibrium state after some time, is quite general and is applicable to the semiclassical Boltzmann equation. The only change needed is to replace the MB distribution in the classical case by the Bose–Einstein or Fermi–Dirac distribution depending on the types of carrier particles considered. Indeed, the semiclassical Boltzmann–BGK equation has been widely used in electron carrier transport in semiconductor [12–19] and phonon energy transport in thermoelectric materials [20]. The main shortcoming of the classical BGK model is that the Prandtl number is fixed and equal to one which is quite different from realistic problems. In order to improve the classical Boltzmann–BGK model equation which does not produce correct Prandtl numbers, an ellipsoidal-statistical (ES) kinetic model was proposed by Holway [21]. An H-theorem for the ES model has been proved [22]. Analogously, a semiclassical ES kinetic model for quantum gases in the normal phase has been recently developed by Wu *et al.* [23]. This new model mimics the classical ES model of Holway [21] and was derived by maximizing the corresponding entropy function under some constraints. It possesses all the desirable features of a kinetic model and yields adjustable Prandtl numbers through a parameter that relates the energy tensor and the pressure tensor. Physically, in the BGK model, the conservative macroscopic variables, density, mean velocity and temperature are involved in the post-collision state, whereas in the ES model, in addition to these conservative normal variables, the local stress tensor quantities are also involved in the post-collision state. In the semiclassical Boltzmann–BGK equation, if the Bose–Einstein (BE) or Fermi–Dirac (FD) distribution is replaced by the semiclassical ES distribution, then one has the semiclassical Boltzmann–ES–BGK equation [23]. It is emphasized here that in the FD or BE statistics, only the normal solution variables such as number density, mean velocity, temperature and chemical potential (or fugacity) are involved; however, in the semiclassical ES distribution, additional pressure tensor variables are also involved.

In theory, one applies the Chapman–Enskog expansion method to the full Boltzmann or model Boltzmann equation (classical or semiclassical) to yield sets of governing hydrodynamic transport equations, such as Euler, Navier–Stokes and Burnett equations which are, respectively, the zeroth-, first- and second-order expansions in the parameter Knudsen number [24]. The zeroth-order solution yields the classical (or quantum) equilibrium distribution in phase space, and its moment integration leads to the classical (or quantum) Euler equations which govern the ideal gas dynamics in physical space. The semiclassical ES equilibrium distribution is the zeroth-order solution of the semiclassical Boltzmann–ES–BGK equation. In recent years, several works on the ideal quantum gas dynamics as governed by the BE and FD equilibrium distribution have been presented [25–27]. Recently, ideal quantum gas dynamics governed by the semiclassical ES distribution was reported [28]. This work is an extension and sequel to Yang *et al*. [28] and the effects of finite relaxation time and non-equilibrium flows are studied. The objectives of this study are twofold. First, we present an accurate numerical solution strategy which includes a decoding method for solving the semiclassical Boltzmann–ES–BGK equation in phase space. Second, we investigate the effect of different ranges of values of relaxation time which correspond to different Knudsen numbers and thus different degrees of gas rarefaction. The accurate numerical method consists of three parts. First is the discrete-ordinate method to discretize the velocity (or momentum space), second is the high-resolution shock-capturing method in physical space and the third one is a decoding procedure to update the flow variables. In the classical rarefied gas flow computation, the implementation of a discrete-ordinate method and high-resolution total variation diminishing (TVD) [29] and essentially non-oscillatory (ENO) schemes [30] to nonlinear model Boltzmann equations has been developed [31] and a similar work using a high-order compact scheme has been proposed [32]. Extension to semiclassical Boltzmann–BGK equation using TVD and weighted ENO (WENO) schemes [33] has been reported [34]. Such a direct method will allow one to examine the same physical flow problems but with different gas or particles. In the past decade, there have been new developments in the numerical methods for solving the hyperbolic conservation law with and without a stiff source term such as the asymptotic preserving algorithm [35,36], discontinuous Galerkin Runge–Kutta methods [37] and recent high-order flux correction procedure via reconstruction (CPR) methods [38,39]. Such high-order methods can also be implemented for the present governing equation. An efficient asymptotic preserving scheme for the semiclassical Boltzmann–BGK model and the full Boltzmann equation has been presented [40]. Two other major methods for solving the Boltzmann–BGK-type equations are the gas-kinetic BGK method [41] and the lattice Boltzmann methods, particularly those that can be applied to high-speed compressible flows [42–44]. Our main purpose here is to explore the flow characteristics inherited in the semiclassical Boltzmann–ES–BGK equation especially in the transitional flow regime but not to exhaust all the possible numerical methods. So we will confine ourselves to our previously employed TVD and WENO methods. In contrast to the equilibrium ES flow presented in [28], the complex non-equilibrium flows will be emphasized especially those high-order moment quantities such as pressure tensors and heat flux vector. Also, if the Euler limit solution (when the relaxation becomes small and approaches zero) of the same flow problem is considered, then one expects to recover similar or identical flow structures as those obtained by solving the quantum Euler equations. Computations of several two-dimensional Riemann problems in gas flows of arbitrary statistics for different range of relaxation times, which correspond to different ranges of Knudsen numbers and thus different flow regimes, are shown to illustrate the complex rarefied gas dynamics as governed by the semiclassical Boltzmann–ES–BGK equation.

This paper is organized as follows. In §2, the semiclassical Boltzmann kinetic theory of Uehling and Uhlenbeck is first described and the Boltzmann–ES–BGK equation and its relation with hydrodynamic equations are outlined. In §3, the governing equations for non-equilibrium flows in two space dimensions and their non-dimensional forms are detailed. In §4, the present numerical solution method in phase space is described. A decoding procedure to determine the parameters for determining the equilibrium ES distribution is given. In §5, computations of two-dimensional semiclassical Riemann problems of rarefied gas dynamical flows are presented for several relaxation times to illustrate the present method as well as the rarefied gas flows. Finally, concluding remarks are given in §6.

## 2. Semiclassical Boltzmann ellipsoidal-statistical kinetic model equation

In this section, the elements of semiclassical Boltzmann kinetic model equations are briefly described. The full semiclassical Boltzmann equation and its related theory can be found in [1,2]. Here we follow the development in [28] and consider the semiclassical Boltzmann ES model equation by Wu *et al.* [23]
*f*(**p**,**x**,*t*) is the distribution function which represents the average number density of particles with momentum **p** at the space–time point (**x**,*t*), *m* is the particle mass, *U* is the mean field potential, *τ* is the relaxation time which needs to be specified for each carrier scattering and *f*^{ES} is the ES equilibrium distribution function given by
*z*(**x**,*t*)=*exp*(*μ*(**x**,*t*)/*k*_{B}*T*(**x**,*t*)) is the fugacity, where *μ*(**x**,*t*) is the chemical potential, *k*_{B} the Boltzmann constant, *T*(**x**,*t*) is the gas temperature, λ_{αβ} are elements of a matrix and its inverse matrix has elements *C*_{α}=(*p*_{α}/*m*)−*u*_{α} are the components of the thermal velocity and *u*_{α} are components of gas mean velocity. Here, for *θ*=+1, BE particles are considered, for *θ*=−1, FD particles, and for *θ*=0, the Maxwell–Boltzmann (MB) classical particles are considered.

The ES equilibrium distribution function *f*^{ES} is derived by the maximum entropy principle under some constraints that the mass, momentum and energy are conserved and the following relations are satisfied:
*D* is the dimension of momentum space and *W*_{αα}=*P*_{αα}=*Dp* is required for the conservation of energy, where *p* is the gas pressure. For further details, see [23]. Here we shall denote *W*_{αβ} as the energy tensor for convenience.

In principle, after solving the semiclassical Boltzmann–ES–BGK equation, equation (2.1), for the (non-equilibrium) distribution function *f*, all the macroscopic dynamic variables of interest such as number density *n*(**x**,*t*), momentum density *n*** u**(

**x**,

*t*), energy density

*ϵ*(

**x**,

*t*), pressure tensor

*P*

_{αβ}(

**x**,

*t*) and the heat flux vector

*Q*

_{α}(

**x**,

*t*) are moments of the distribution function and are defined by

*W*

_{αβ}to completely determine

*f*

^{ES}and require that

*W*

_{αα}=

*P*

_{αα}to insure conservation of energy, a natural choice leading to rotational invariance of the kinetic model equation is [21]

*p*is the hydrodynamic pressure defined as the average of the diagonal components of the pressure tensor

*P*

_{αβ}, i.e.

*p*=

*P*

_{αα}/

*D*, and

*b*is a parameter. Requiring

*W*

_{αβ}to be positive definite, the parameter

*b*satisfies

*D*=2 and

*D*=3. Finally, to determine the semiclassical ES distribution

*f*

^{ES}, one usually needs to calculate

*z*and the matrix λ

_{αβ}simultaneously through

*ν*.

It is worth noting that in the semiclassical ES model not only the conservative quantities number density *n* (or *ρ*=*mn*), mean velocity ** u** and temperature

*T*are involved in the post-collision equilibrium as in

*f*

^{(0)}(defined below), but also the other higher order quantities

*W*

_{αβ}(energy tensor) which in turn is related to pressure tensor

*P*

_{αβ}, involved in the post-collision equilibrium as in

*f*

^{ES}.

It is also noted that in equation (2.1), if *f*^{ES} is replaced by *f*^{(0)}, the standard BE or FD equilibrium distribution, then one has the usual semiclassical Boltzmann–BGK equation. Here, the standard BE or FD equilibrium distribution, *f*^{(0)}, is given by
*θ*=−1, denotes the FD statistics, *θ*=+1, the BE statistics and *θ*=0 denotes the MB statistics.

The conservation laws of macroscopic properties can be obtained by multiplying equation (2.1), respectively, by 1, **p** and **p**^{2}/2*m*, and integrating the resulting equations over all **p**. Consequently, the integrals of the collision terms in all three cases vanish automatically resulting in the conservation laws in differential equations form for the conserved macroscopic quantities, i.e. number density *n*(**x**,*t*), momentum density *mn***u**(**x**,*t*) and energy density *ϵ*(**x**,*t*).

We note that in the *semiclassical ES Euler* limit in which the particle distribution function is always in reference equilibrium, i.e. *f*=*f*^{ES} and the collision term of equation (2.1) vanishes automatically. In this case, the resulting set of governing macroscopic conservation equations are called the semiclassical ES Euler equations. We also point out that *f*^{ES} is not isotropic and has quantities like *P*_{αβ} (or *W*_{αβ}) and is quite different from *f*^{(0)}, the standard BE and FD distributions. In recent years, notable numerical methods to describe the ideal quantum gas flows have been developed [25–27]. Numerical experiments with *f*=*f*^{ES} in equation (2.1), corresponding to relaxation time *τ*=0, have been presented [28]. This work is a sequel to Yang *et al*. [28] and the effect of finite relaxation time which covers different non-equilibrium flow regimes is examined. It is worth pointing out that in the equilibrium flow (*τ*=0), one only has the unknown *f*^{ES} and the macroscopic quantities from *f*^{ES}, *n*(**x**,*t*),** u**(

**x**,

*t*) and

*W*

_{αβ}(

**x**,

*t*), are given by equation (2.3). On the other hand, for the non-equilibrium flows, one has to solve the unknown distribution function

*f*(

**p**,

**x**,

*t*) and the macroscopic quantities are given by equation (2.4) with

*W*

_{αβ}(

**x**,

*t*) defined by the parameter

*b*and the pressure tensor

*P*

_{αβ}(equation (2.5)).

## 3. Governing equations in two space dimensions

In this section, we describe in detail the theory and governing equation in two space dimensions. We also neglect the influence of the externally applied field *U*(**x**,*t*). The governing equation in two space dimensions is formulated directly without introducing the reduced distribution functions. It is noted that quantum statistics of ideal gases in two dimensions has striking features such as the specific heat for an ideal gas of Fermi particles is identical with that for an ideal Bose gas [45]. The semiclassical Boltzmann–ES–BGK equation in two space dimensions can be expressed as
*υ*_{x} and *υ*_{y} as particle velocities. The macroscopic properties such as number density, mean velocity, energy, pressure tensor and heat flux vector are the various moments of *f*(*υ*_{x},*υ*_{y},*x*,*y*,*t*)
*p*(*x*,*y*,*t*)=(*P*_{xx}+*P*_{yy})/2. For a specified value of parameter *b*, the tensor *W*_{αβ}(*x*,*y*,*t*) can be obtained through
*W*_{xx}+*W*_{yy}=*P*_{xx}+*P*_{yy} is required. To obtain the new *f*^{ES}, one needs *z* and λ_{xx},λ_{xy},λ_{yy} and these can be determined through solving the following equations simultaneously:
*b*=0, one recovers the semiclassical BGK model, i.e. *f*^{ES} becomes equal to *f*^{(0)}.

### (a) Governing equations in non-dimensional form

Before proceeding to discretize the equation, in this section we introduce the characteristic properties of *V* _{0}, *t*_{0} and *n*_{0} (or *z*_{0}) for the purpose of normalization. Choose some reference parameters as
*L* defined as the characteristic length of the problem. Hence the definitions of non-dimensional variables are introduced as

## 4. Solution algorithms in phase space

### (a) Implementation of discrete-ordinate method

In the two-dimensional case, we may apply the Gauss–Hermite quadrature rule over the interval *υ*_{α} and the corresponding weights *W*_{α}, with *α*=*σ*,*δ*, are tabulated in the table of the Gauss–Hermite quadrature [46]. Applying the discrete-ordinate method to equation (3.7) with (*υ*_{α}), then we have a set of partial differential equations in physical space (*x*,*y*,*t*) instead of a single equation in phase space. By casting the equations in conservation law form, we have
*f*_{σ,δ}(*x*,*y*,*t*)=*f*(*υ*_{σ},*υ*_{δ},*x*,*y*,*t*) and similarly for *f*^{ES}_{σ,δ}(*x*,*y*,*t*). Once the discrete distribution functions *f*_{σ,δ}(*x*,*y*,*t*) are solved, for every time level we can acquire and update the macroscopic moments in physical space using a selected quadrature method as described below
*α*=*x*, one sets *υ*_{α}=*υ*_{σ}, *u*_{α}=*u*_{x}, *C*_{σ}=*υ*_{σ}−*u*_{x} and when *α*=*y*, one has *υ*_{α}=*υ*_{δ}, *u*_{α}=*u*_{y}, *C*_{δ}=*υ*_{δ}−*u*_{y}. Once the values of *P*_{αβ} have been obtained, the gas pressure can be computed by *p*(*x*,*y*,*t*)=(*P*_{xx}+*P*_{yy})/2 and the tensor components *W*_{αβ} by equation (3.3).

Here, note that the selection criteria of a suitable quadrature rule is first to guarantee that the non-equilibrium distribution function *f* can be accurately represented with suitable velocity range being covered and second its macroscopic moments can be accurately computed. In this work, Gauss–Hermite quadrature with a suitable velocity range was employed, however, the equally spaced Newton–Cotes formula can also be used if one needs to cover the high-energy tail of the distribution function depending on the physical problem.

### (b) Implementation of high-resolution shock capturing schemes

In this section, we describe the numerical algorithm to solve the set of equation (4.3). A class of TVD schemes is implemented to solve the set of conservation laws. We firstly discretize the space *x*,*y* and time *t* into a number of uniform cells centred at *i*,*j* at time *n*, hence we approximate *f*_{σ,δ}(*x*,*y*,*t*) by *f*_{σ,δ} from *n*th time level to (*n*+1)th time level is expressed by
*x*- and *y*-directions, respectively. For the first-order upwind method, the numerical fluxes are defined by
*r*=3)) can be found in [34] and will not be repeated here. The class of asymptotic preserving schemes for the kinetic equations and related problems with stiff sources can be beneficially applied as well [36,40].

### (c) Determination of *z* and λ_{αβ}

After acquiring the macroscopic moments, we need to update the values of *z* and λ_{αβ} in order to obtain the *f*^{ES} distribution function for the next time step. Using equation (3.4) and for a specified *b*, we can solve *z*(*x*,*y*,*t*), which is the root of equation
*b*=0, one recovers the BGK model, i.e. *f*^{ES} becomes equal to *f*^{(0)}. After the new value of *z* is determined, we can immediately find λ_{αβ} before advancing the solution to the next time step.

After obtaining λ_{αβ}, we can calculate *f*^{ES} distribution function.

## 5. Numerical results

In this section, we present the results of selected computational tests. Here, we test our solution method by solving two-dimensional Riemann problems of semiclassical rarefied gas dynamics. In order to compare the present results of *τ*→0 with those of Yang *et al.* [28] which corresponding to *τ*=0, we have employed the same numerical method and the same mesh system as in [28]. Although the initial set-up of the present two-dimensional Riemann problem is the same as that in [28], however, due to finite relaxation time, the dissipative effects are not negligible and will cause the evolution and interaction of discontinuities to display more broad profiles in general depending on the value of relaxation time (or Knudsen number). The solution procedure for finite relaxation time non-equilibrium flow is different from that of *f*=*f*^{ES} (*τ*=0) equilibrium flows, namely, in addition to the usual conservative quantities *n*,** u**,

*T*and

*W*

_{αβ}, we also have pressure tensor

*P*

_{αβ}and heat flux vector

**and beyond. The degree of rarefaction of a gas can be characterized by the Knudsen number which is defined as**

*Q**L*is the characteristic length,

*a*is the

*s*-wave scattering length,

*n*

_{0}is the referenced number density. Applying the Chapman–Enskog expansion to the semiclassical Boltzmann equation, one can obtain the thermal conductivity and shear viscosity and the relaxation time can be shown to be related to the Knudsen number by

*z*=0, it reduces to

*Kn*number for a chosen value of relaxation time. If the temperature, pressure and fugacity are kept the same for a flow, then the relaxation time is proportional to the Knudsen number.

### (a) Two-dimensional Riemann problems

We numerically study the two-dimensional Riemann problems for rarefied quantum gas dynamics for several relaxation times using the present direct solver for the semiclassical Boltzmann–ES–BGK equation. In the limit *τ*→0, i.e. *Kn*→0 limit, then *f*→*f*^{ES}, we can recover the ideal gas dynamics governed by the quantum ES Euler equations. Following the works of earlier studies [47–49], we have selected several configurations to be tested among those 19 configurations classified. We also follow the notations used there, namely, *J*, denotes slip line, *R*, rarefaction wave and *S*, shocks and forward arrow indicates forward waves, and backward arrow indicates backward waves. In the Riemann problems of two-dimensional gas dynamics, initially the state in each of the four quadrants is in equilibrium and kept constant, i.e. *f*=*f*^{ES}=*f*^{(0)} in each quadrant, see the illustration in figure 1. The set-up is made so that only one elementary wave, a one-dimensional shock, ** S**, a one-dimensional rarefaction wave,

**or a two-dimensional slip line (or contact discontinuity),**

*R**J*, appears at each interface. The computational domain is (

*x*,

*y*)∈[0,1]×[0,1] and is divided into equal intervals with Δ

*x*=Δ

*y*. After some mesh refinement tests, we found that 200×200 grids with Δ

*x*=Δ

*y*=0.005 and the Gauss–Hermite quadrature rule with 20 abscissas both in

*υ*

_{x}and

*υ*

_{y}directions will yield good resolution of the major flow features and all the results below were obtained with these mesh and discrete velocity system. We examine the effects due to (i) the finite relaxation time which covers different flow regimes, (ii) the model parameter

*b*, which adjusts the Prandtl number in the flow and (iii) the three statistics. Specifically, we consider two configurations, namely the

*Configuration 5*and

*Configuration 17*among the 19 genuinely different configurations classified in the earlier studies [47–49]. To satisfy these constraints, the set of initial data must satisfy the Rankine–Hugoniot relations, i.e. for given a pair of left and right states (denoted by the subscripts

*l*and

*r*), we define the two configurations as follows.

### (b) Results of *Configuration 5*

*Configuration 5*. This is the four slip lines case with motion in the opposite directions. We have
*E*_{ij} denotes an elementary wave *E* between the *i*th and *j*th quadrant, and *E* can be either *J*, ** S** or

**. The initial constant states (**

*R**z*,

*u*

_{x},

*u*

_{y},

*T*) at the four quadrants are specified as

The general inviscid flow (*τ*=0 or very small) development is outlined below. After some time *t*>0, the slip lines *J*_{21} and *J*_{32} will meet the sonic circle of the constant state in the second quadrant and continue as almost straight lines so that a quarter of the sonic circle lies in between. Identical symmetric flow patterns for the slip lines *J*_{34} and *J*_{41} occur in the fourth quadrant. Inside the subsonic area the slip lines bend and end in spirals. Centred about the point *τ* is increasing, it is expected that the inviscid flow features will become smeared, even broaden and eventually hard to identify and delineate.

We first report the results for the *Configuration 5* for FD (*θ*=−1) statistics with three values of relaxation times. In figure 2, the number density, pressure, fugacity and pressure tensor contours are shown for the FD statistics for the relaxation time *τ*=1/10 000 with parameter *b*=0.5. The output time is at about the instant similar to that reported in [47,48]. Note that we are solving the non-equilibrium distribution function in contrast to those solve the classical Euler equations [47,48]. The results for the case of *τ*=1/10 000 approximate those of the quantum Euler solution as given in [28]. From the density and pressure contours, one can identify the interesting and complicated wave patterns. At the right top corner, the two slip lines *J*_{21} and *J*_{32} meet the sonic circle of the constant state in the second quadrant and continue as almost straight lines so that a quarter of the sonic circle lies in between. Our results of *τ*=1/10 000, which is close to the Euler solution ( *f*=*f*^{ES}), are consistent with calculations in [47,48] and with the theoretical conjecture in [49] for the same configuration and with the same initial data. It can be easily observed that the values of *P*_{xx} and *P*_{yy} are of the same order and are larger than *P*_{xy}. The pressure, *P*_{xx} and *P*_{yy} contours are of similar flow patterns. The tensor quantities *W*_{αβ} (not shown here) can also be calculated from the parameter *b* and pressure tensor *P*_{αβ}. The present results of the case of *τ*=1/10 000 are in good agreement with those given in [28,47,48] and serve as a validation example of this theory and solution method.

In figure 3, the number density, pressure and fugacity contours for the case of relaxation time *τ*=1/1000 with parameter *b*=0.5 are shown for BE (*a*–*c*), MB (*d*–*f*) and FD (*g*–*i*) statistics, respectively. The general flow patterns and structures are similar to that of *τ*=1/10 000 as shown in figure 2*a*–*c* for the same configuration except for slightly smeared profiles. There are detectable differences among the three statistics although the overall wave patterns are similar. Quantitatively, comparing the same contours for the three statistics, the numerical values of number density and pressure for the BE statistics are the largest among the three statistics, and the MB statistics always lie between the other two as dictated by the *θ* values. We further increase the relaxation time and examine the flow patterns. In figure 4, the number density, pressure and fugacity contours for the case of relaxation time *τ*=1/100 with parameter *b*=0.5 are shown for all three statistics, respectively. The flow patterns and structures are much broader when compared with those of figure 3. With increasing relaxation time, the corresponding Knudsen number is increasing and the effect of rarefaction becomes more evident and the discontinuities will exhibit progressively blurred profiles although the basic flow structures still can be identified. The effect of finite relaxation time can be partially examined from the results shown in figures 2–4. The two limiting extremes

### (c) Finer mesh solution and WENO3 solution

Here, we study the consistency and convergence of our solution method by increasing the number of grids in the computational domain and therefore reducing their size. Here we expect a decrease of the truncation error, and therefore observe a converging pattern to the solution of the chosen configuration under test. We report in figure 5 results for *Configuration 5* with relaxation time *τ*=1/1000 and parameter *b*=0.5 using a grid system of 400×400, i.e. four times that used in our initial test grid. We observe in figure 5*a*,*b* that slip line structures exhibit slightly sharper profiles when compared with figure 3*g*,*h*. This is also true in the vicinity of rarefaction fans and contact discontinuities. The higher order moment pressure tensor quantities also display the finer details of the flow patterns. With the same *τ*=1/1000 and same FD statistics, by comparing the number density and pressure contours of FD in figure 3 (200×200) with the present 400×400 results, a slightly finer resolution of flow profiles and values can be detected for the 400×400 results. Both results in figures 5 and 3 are computed with the *CFL* condition equal to 0.6 and corresponds to an output time of *t*=0.25. The trend of convergence of the solutions as mesh size is refined can be observed.

We also report some selective solutions using fifth-order WENO3 scheme and a grid system of 200×200 for both BE and FD statistics with *b*=0.5. The pressure contours together with the streamline patterns for the FD statistics for *Configuration 5* are shown in figure 6*a*, and the temperature contours together with the heat flux streamlines are shown in figure 6*b*. The corresponding results for the BE statistics for *Configuration 17* are shown in figure 6*c*,*d*. The description of *Configuration 17* is given below.

### (d) Results of *Configuration 17*

Next, we consider another different set-up of initial conditions, resulting in different flow structures, namely, *Configuration 17* of the two-dimensional Riemann problem.

*Configuration 17*. In this configuration, we have flow structures of two slips, 1 shock and 1 rarefaction initially

Here, the functions *Φ*_{41} and *Ψ*_{32} are the same as that defined in [48]. The initial conditions at the four quadrants are specified, respectively, by
*z*,** u**,

*T*) specified for each quadrant at time

*t*=0, one can then define the equilibrium distribution

*f*

^{ES}=

*f*

^{(0)}for each quadrant at

*t*=0 as the initial condition. Thus, the two-dimensional Riemann problem is completely specified on the distribution function and one starts marching in time the solution of the distribution function. This is rather different from the traditional two-dimensional Riemann problems, [47,48].

In this configuration, shock, slip line and rarefaction fan are conjectured to be propagating from interfaces of the four quadrants. A shock will be seen propagating through quadrants 3 and 2, rarefaction fan through quadrants 4 and 1 and slip discontinuities through quadrants 2 and 1 and through quadrants 4 and 1. Similar to the *Configuration 5*, as *τ* is increasing, it is expected that the inviscid flow features will become broader, even blurred and eventually hard to identify and delineate.

Now, we report the results for *Configuration 17*. The number density and pressure contours for the three statistics with relaxation time *τ*=1/1000 and parameter *b*=−0.5 are shown in figure 7. After some time (non-dimensional time *t*=0.2), as shown in the figure, the two slip lines *J*_{21} and *J*_{34} part the wave patterns into a left and a right part and connect with a vortex inside the subsonic region. For *Configuration 17*, again, the wave patterns shown in figure 7 are consistent with calculations in [47,48] and with the theoretical conjecture in [49] for the same configuration and with the same initial data. Notable differences among the three statistics can be observed although the overall wave patterns are similar. Quantitatively, the numerical values of number density and pressure contours for the BE statistics are the largest among the three statistics and the flow properties obtained by using MB statistics are always greater than those obtained with FD statistics and lower than those with BE statistics as it should be. The pressure tensor quantities *P*_{xx}, *P*_{xy} and *P*_{yy} for *Configuration 17* for the three statistics with *τ*=1/1000 and *b*=−0.5 are also studied and shown in figure 8. For each statistic, the *P*_{xx} and *P*_{yy} display similar flow patterns and values and their values are one order of magnitude larger than that of *P*_{xy}.

It is also interesting to investigate quantitatively and qualitatively the effects of parameter *b* in the ES distribution model, *f*^{ES}, on the wave flow patterns and properties. From qualitative observations on different configurations, it is shown that the parameter *b* can influence the results by increasing or decreasing the sharpness of the profiles of the two configurations. Quantitatively, we note that it is not certain how this parameter may affect the overall magnitude of flow properties. The results for *Configuration 17* with parameter *b*=0.5 for three statistics are shown in figures 9 and 10, which correspond to figures 7 and 8, respectively. One can first observe that the contours of all the flow variables in the case of *b*=0.5 are similar to their corresponding one in the above case of *b*=−0.5. In general, the upper (lower) bounds of number density and pressure are slightly lower (higher) for *b*=−0.5 than that of *b*=0.5 for three statistics and the upper bounds of *P*_{xx} and *P*_{yy} are lower and the lower bounds of *P*_{xy} are higher for *b*=−0.5 when compared with that for *b*=0.5.

Comparing the pressure contours with that of *P*_{xx} and *P*_{yy}, as before, we observe that the pressure and the main diagonal components, *P*_{xx} and *P*_{yy}, exhibit similar flow patterns but again slight numerical differences can be identified. Note that the main diagonal components will be equal to each other and equal to pressure, i.e. *P*_{xx}=*P*_{yy}=*p*, when the dimensionless parameter *b*=0 and *f*^{ES}=*f*^{(0)}. Based on the quantitative observations on several tested configurations, it is safe to conjecture that due to the slight differences on *P*_{xx} and *P*_{yy} and the contribution of *P*_{xy} components, as triggered by the dimensionless parameter *b* of the ES distribution function, the anisotropic ES model contributes to relax the fixed Prandtl number condition of the traditional equilibrium distribution model, *f*^{(0)}. It is also noted that when the flow is more towards the degenerate case then the difference between *P*_{xx} and *P*_{yy} will become larger (more anisotropic).

## 6. Concluding remarks

Numerical solutions of the semiclassical Boltzmann ES kinetic model equation proposed by Wu *et al.* [23] have been obtained for two-dimensional Riemann problem in rarefied gas of particles of arbitrary statistics. The semiclassical ES equilibrium distribution is anisotropic and was derived through the maximum entropy principle and involves not only the mass, momentum and energy, but also pressure tensor quantities and differs from the standard BE or FD distribution. Therefore, the gas dynamical behaviours away from ES equilibrium are not well known and are worth exploring. The equilibrium gas dynamical flows governed by the semiclassical ES equilibrium distribution have been explored in detail in [28]. Here, the unsteady non-equilibrium rarefied quantum gas dynamical flows are numerically studied which is a generalization and sequel to Yang *et al*. [28]. The computational method treats the governing equation in phase space and employs the discrete ordinate method and high-resolution shock capturing schemes. Both second-order TVD and fifth-order WENO schemes are employed. Specifically, we describe the solution method in details for the equation in two space dimensions. A different decoding procedure is devised for the semiclassical ES distribution. Computations of two-dimensional Riemann problems for rarefied gas flows of arbitrary particle statistics are presented for several relaxation times which correspond to various ranges of Knudsen numbers. A mesh refinement test has been conducted to show solution convergence. Our results for small Knudsen number which approximates the Euler limit are in good agreement with the calculations in [47,48] and consistent with the theoretical conjecture in [49]. The effect of the parameter *b* of the ES model on flow fields is examined. Results showed how this parameter interacts with the pressure tensor and influences the Prandtl number in the non-equilibrium flows and therefore the transport properties of classical and quantum gases. These computational examples serve the purpose of exploring the nonlinear manifestation of shock wave, contact line and rarefaction wave and testing the robustness of the present method. All the expected flow profiles comprising shock, rarefaction wave and contact discontinuities of semiclassical ideal gases and their nonlinear interactions can be observed with considerably good detail and are consistent with available results. The finite relaxation time introduces dissipative effects into the flows as well as non-equilibrium flows beyond the Navier–Stokes level. This work emphasizes the building of the unified and parallel framework for treating semiclassical rarefied gas dynamics of three statistics. The extension to boundary value problems with general partially specular reflection and partially diffusive reflection will be the next endeavour.

## Funding statement

This work was done under the auspices of National Science Council, Taiwan through grant nos. NSC-99-2221-E002-084MY3, NSC-100-2221-E002-106MY3 and CASTS 10R80909-4. The authors are grateful to Dr Bagus P. Muljadi for the early work on semiclassical Boltzmann–BGK equation. Z.L. is supported by National Nature Science Foundation, China under grant nos. 91016027 and 11325212.

## Acknowledgements

This work was initiated during the senior author's visit to No. 5 Institute of CARDC during the summer of 2012.

- Received January 23, 2014.
- Accepted May 22, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.