The unsteady shock wave diffraction by a square cylinder in gases of arbitrary particle statistics is simulated using an accurate and direct algorithm for solving the semiclassical Boltzmann equation with relaxation time approximation in phase space. The numerical method is based on the discrete ordinate method for discretizing the velocity space of the distribution function and high-resolution method is used for evolving the solution in physical space and time. The specular reflection surface boundary condition is employed. The complete diffraction patterns including regular reflection, triple Mach reflection, slip lines, vortices and their complex nonlinear manifestations are recorded using various flow property contours. Different ranges of relaxation times corresponding to different flow regimes are considered, and the equilibrium Euler limit solution is also computed for comparison. The effects of gas particles that obey the Maxwell–Boltzmann, Bose–Einstein and Fermi—Dirac statistics are examined and depicted.
The Boltzmann equation has been widely used to describe various transport phenomena in classical rarefied gas covering a wide range of flow parameters, such as Reynolds number, Mach number and Knudsen number. The Chapman–Enskog expansion method is usually applied to the Boltzmann equation to derive closed set hydrodynamic transport equations that apply to a broad range of flow regimes (Chapman & Cowling 1970). In analogy to the classical Boltzmann equation, a semiclassical Boltzmann equation, which generalizes the collision term to treat the collision of particles obeying quantum statistics, was developed by Uehling & Uhlenbeck (1933). The hydrodynamic behaviour of quantum gases has been the subject of some prominent research—see Nikuni & Griffin (1998), Wigner (1932), and the application of quantum Boltzmann hydrodynamic equations have been implemented in the analysis of electron flows in quantum semiconductor devices (Ancona & Iafrate 1989; Gardner 1994). In recent years, owing to the rapid advancements of microtechnology and nanotechnology, the device or structure characteristic length scales have become comparable to the mean free path and the wavelength of energy and information carriers (mainly electrons, photons, phonons and molecules), and some of the classical continuum transport laws are no longer applicable. It is generally believed that the microscopic description of the Boltzmann equation (classical and semiclassical) is adequate to treat transport phenomena in the mesoscale range. Besides, owing to the different types of carriers involved simultaneously in a single problem, it is desirable to have a method that can allow one to treat them in a unified and parallel manner (Chen 2005). With the semiclassical Boltzmann equation, it is possible to describe adequately the mesoscale transport of particles obeying arbitrary statistics.
The main difficulty encountered in solving the semiclassical Boltzmann equation as derived by Uehling & Uhlenbeck is the same as that encountered in the classical counterpart and is mainly due to the complicated integral nature of its collision term. The relaxation time approximation proposed by Bhatnagar et al. (1954) for the classical Boltzmann equation provides a much simpler form of collision term and retains the principal effects of particle collisions and enables more tractable solution methods. The Bhatnagar–Gross–Krook (BGK) relaxation time concept is quite general and is applicable to the semiclassical Boltzmann equation as well. The only change is that the equilibrium Maxwell–Boltzmann distribution in the classical case is replaced by the Bose–Einstein or Fermi–Dirac distribution depending on the types of carrier particles. In fact, the semiclassical Boltzmann–BGK equation has been widely used for electron carrier transport in semiconductors, see Lundstrom (2000), Fatemi & Odeh (1993), Carrillo et al. (2003a,b), Majorana & Pidatella (2002) and Pattamatta & Madnia (2009), and phonon energy transport in thermoelectric materials by Chen (2005). Recently, by following and extending the work of Yang & Huang (1995) developed for the classical Boltzmann–BGK equation, an accurate and direct numerical method for solving the semiclassical Boltzmann–BGK equation has been presented for initial value problems in one dimension (Muljadi & Yang 2010). Several aspects of the algorithm including the range of relaxation times and different range of Knudsen numbers have been tested. The novel algorithm enables one to treat the particles of three statistics on an equal footing and in a parallel manner.
In this work, we apply a newly developed accurate and direct solver for the semiclassical Boltzmann–BGK equation to simulate the complex unsteady shock wave diffraction by a two-dimensional cylindrical body placed in a rarefied gas of particles of arbitrary statistics. The unsteady shock diffraction phenomena generated by an incident shock wave on a body involve complex regular reflection, triple Mach reflection, slip lines, vortices and their complex interactions, and provide a rich manifestation of nonlinear waves in gas dynamics. In the past, many experimental and numerical studies for such shock wave diffraction problems have been reported for the classical gases (e.g. Bryson & Gross (1961) and Yang et al. (1987)) and they serve as good examples for testing the capability of a numerical method for solving gas dynamical problems. The total variation diminishing (TVD) scheme is adopted in this paper to describe the evolution of these equations. We consider the initial and boundary value flow problems in two space dimensions. In particular, the Maxwell-type boundary condition with the concept of an accommodation coefficient at the wall surface is adopted. The complete shock wave diffraction patterns are simulated and studied for the same flow conditions but with three different particle statistics.
Elements of the semiclassical Boltzmann–BGK equation are briefly described in §2. Its correlation with hydrodynamic equations is outlined. In §3, we describe the numerical method, including the use of the discrete ordinate method to discretize the particle distribution function into a set of hyperbolic conservation laws with source terms, the quadrature method and the second-order TVD implementation. In §4, the initial and boundary conditions to specify the physical flow problem are described. In §5, detailed numerical simulations of unsteady shock wave diffraction by a square cylinder are reported. Various flow contours are depicted to add better understanding of the phenomena. Finally, concluding remarks will be given in §6.
2. Semiclassical Boltzmann–Bhatnagar–Gross–Krook equation
In this section, the elements of the semiclassical Boltzman–BGK equation that are relevant to the present development are briefly described. We consider the extension of the Boltzmann equation to quantum systems due to Uehling & Uhlenbeck in which they took the Pauli exclusion principle into account. 2.1 where m is the particle mass, U is the mean field potential and f(p,x,t) is the distribution function that represents the average density of particles with momentum p at the space–time point x,t. The (δf/δt)coll. denotes the collision term and, according to Uehling & Uhlenbeck (1933), it takes the form 2.2 where K(p,q,Ω) denotes the collision kernel and Ω is the solid angle, and θ is a parameter that specifies the type of particle statistics. Here, for θ=+1, Bose–Einstein particles are considered; for θ=−1, Fermi–Dirac particles; and for θ=0, the Maxwell–Boltzmann classical particles are considered. It is well known that the collision integral of the semiclassical Boltzmann equation will automatically vanish when the gas flow is in equilibrium. The equilibrium distribution function for general statistics can be expressed as 2.3 where u(x,t) is the mean velocity, T(x,t) is temperature, kB is the Boltzmann constant and is the fugacity, where μ is the chemical potential. In (2.3), θ=−1, denotes the Fermi–Dirac statistics; θ=+1, the Bose–Einstein statistics; and θ=0, the Maxwell–Boltzmann statistics. We note that even with the case of θ=0, we still have chemical potential μ or fugacity z, which is rather different from the usual classical Maxwellian distribution. To avoid the mathematical difficulty caused by the nonlinear integral collision term, the relaxation time concept of Bhatnagar et al. is generally applied to replace the collision term of Uehling & Uhlenbeck, thus the semiclassical Boltzmann–BGK equation reads 2.4 Here, τ is the relaxation time and needs to be specified for each carrier scattering.
The macroscopic dynamic variables of interest such as number density, momentum density and energy density are the first few low-order moments of the distribution function and are defined by 2.5 2.6 and 2.7 Other higher-order moments such as the stress tensor Pij and the heat flux vector Φi(x,t) can also be defined accordingly as 2.8 and 2.9
The conservation laws of macroscopic properties can be obtained by multiplying equation (2.1), respectively, by 1, p and p2/2m 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 equation form for the conserved macroscopic quantities, i.e. number density n(x,t), momentum density ρu(x,t) and energy density ϵ(x,t) as follows: 2.10 2.11 and 2.12
(a) Fermi and Bose functions in relation to hydrodynamic properties
Before we proceed, without losing generality, we neglect the influence of an externally applied field U(x,t). To illustrate the present method, we formulate the two-dimensional equilibrium distribution function as 2.13 In a closed form in terms of quantum functions, we replace the distribution function f with f(0), which automatically reduces the source term in the BGK Boltzmann equation to zero. The macroscopic moments, i.e. number density n(x,y,t), momentum j(x,y,t) and energy density ϵ(x,y,t), are given by 2.14 2.15 2.16 and 2.17 where is the thermal wavelength and γ the ratio of specific heats, whereas β=1/kBT(x,t). The functions of order ν are, respectively, defined for Fermi–Dirac and Bose–Einstein statistics as 2.18 and 2.19 Here, applies for the Fermi–Dirac integral and for Bose–Einstein’s, whereas Γ(ν) is the gamma function. The definition of macroscopic quantities in terms of Fermi or Bose functions applies for both cases of quantum distributions.
Semiclassical Euler and Navier–Stokes equations can be obtained from the zeroth-order and first-order solutions of the Boltzmann–BGK equation, equation (2.4), via the Chapman and Enskog expansion (Chapman & Cowling 1970). Also, the transport coefficients such as viscosity η and thermal conductivity κ can be derived in terms of the relaxation time as 2.20 and 2.21 where P here is pressure.
In the following, we will also consider the case of the semiclassical Euler limit in which the particle distribution function is always in equilibrium, i.e. f=f(0), and the collision term of equation (2.1) vanishes automatically. In recent years, kinetic numerical methods for solving the ideal quantum gas flows based on f=f(0) have been developed (Yang & Shi 2006; Yang et al. 2007; Shi & Yang 2008).
(b) Classical limit
Note that for z≪1, both functions for all ν behave like z itself, i.e. . Physically, a large and negative value of chemical potential of a dilute system and high temperature environment do correspond to small z. The hydrodynamic properties of classical case can be acquired by replacing the function into z itself. For θ=0, the distribution function becomes 2.22 In this case, no approximation for is required and the macroscopic values can be obtained by 2.23 In other words, in the classical limit, 2.24 and 2.25
Before proceeding to discretize the equation, in this section, we introduce the characteristic properties of and for the purpose of normalization. The characteristic velocity and time can be defined as 2.26 with L defined as the characteristic length of the problem. Hence the definitions of non-dimensional variables are introduced as 2.27 Hence the normalized semiclassical Boltzmann–BGK equation, 2.28 with and as particle velocities. The normalized two-dimensional semiclassical equilibrium distribution function becomes 2.29
3. Numerical method
(a) Application of discrete ordinate method
By applying the discrete ordinate method to equation (2.28), the distribution function in phase space f(υx,υy,x,y,t) can be rendered into a set of hyperbolic conservation equations with source terms for fσ,δ(x,y,t) in the physical space, where σ=−N1,…,−1,1,…,N1 and δ=−N2,…,−1,1,…,N2. The resulting equations are a set of 3.1 with fσ,δ, υσ and υδ representing the values of f, υx and υy evaluating the discrete velocity points σ and δ, respectively.
(b) Quadrature methods
In the two-dimensional case, we may apply the Gauss–Hermite quadrature rule over the interval . The Gauss–Hermite quadrature rule reads, 3.2 or 3.3
The discrete points υα and weight Wα, with α=σ,δ, are tabulated in the table of the Gauss–Hermite quadrature, see Abramowitz & Stegun (1964) and can be found through 3.4 where l is number of quadrature points and υα are the roots of the Hermite polynomial Hl(υ).
Once the discretized functions fσ,δ(x,y,t) are solved for every time level, we can acquire and update the macroscopic moments in the physical space using the same quadrature method as described below 3.5 3.6 3.7 and 3.8
Here, we note that the basic criteria of choosing different quadratures is to guarantee that the macroscopic moments can be accurately computed; this requires an accurate representation of the distribution function with suitable coverage of velocity range. The equally spaced Newton–Cotes formula can also be used if one needs to cover the high energy tail of the distribution function.
(c) Application of total variation diminishing scheme
In this section, we describe the numerical algorithm to solve the set of equations (3.1). A class of TVD schemes with van Leer’s limiter (van Leer 1979; Harten 1983; Sweby 1984) is implemented to solve the set of conservation laws with source terms. We firstly discretize the space x,y and time t into a number of cells centred at i,j at time n; hence we approximate fσ,δ by . In terms of numerical fluxes, the evolution from nth time level to (n+1)th time level is expressed by 3.9 where the time step Δt is chosen to be less than the local relaxation time, τ. The numerical fluxes are defined by a flux-limiter method with 3.10 where are the low-order upwind flux given by 3.11 with the high-order Lax–Wendroff scheme given by 3.12
Here, and similarly, , whereas for the limiter function ϕ(θxi+1/2,j), van Leer’s limiter is chosen and is given by 3.13 The formulation of can be done in a similar fashion.
After acquiring the macroscopic moments, we need to update the values of z and T. To solve z(x,y,t), which is the root of equation 3.14 we can use numerical root finding methods. After z is acquired, we can immediately find T before advancing to the solution of the next time step.
4. Initial and boundary conditions
In this study, we consider an incident moving shock wave that is located in front and to the left of a square cylinder initially and starts to move to the right and impinge upon the square cylinder and encounter the complex shock diffraction processes. The initial condition at state 1 (ahead of the moving shock) and state 2 (behind the moving shock) for the incident shock wave shall obey the Rankine–Hugoniot condition, i.e. 4.1 4.2 and 4.3 where G=(γ+1)/(γ−1), Ms is the incident shock Mach number and c=(γP/n)1/2 is the local sound speed. For a two-dimensional case, we use γ=2.
The boundary condition needs to relate the distribution function of the incoming particle fI and that of the outgoing particle from the wall fO. When the specularly reflecting wall is considered, it is assumed that a particle does not exchange energy with the wall; the wall is adiabatic and no shear stress acts on the wall. Hence, the outgoing distribution function fO would follow analogously from the incoming distribution function with accordingly transformed normal velocities fI(−υn), i.e. 4.4 where υn=υ⋅n and n is the outward unit normal to the stationary solid surface. For the case of a perfectly thermalizing wall, the particle is subsequently emitted in a distribution fW that is determined by wall temperature TW and velocity uW, i.e. 4.5 Maxwell (1879) describes an intermediate between a specular type of reflection to a perfectly thermalizing/absorbing wall. From this model, the distribution function in the infinitesimal neighbourhood of the wall is given by 4.6 where χ is the accommodation coefficient ranging from 0 to 1. The boundary condition needs to satisfy a zero mass flux normal to the surface at the wall because the wall does not accumulate particles, i.e. 4.7 therefore, for a semiclassical gas, it follows that 4.8 This equation shall be solved to find zW. In the case of a classical gas, the value of nW needs to be solved instead (Yang & Huang 1995; Struchtrup 2005).
5. Computational results
In this section, simulation results of the complete shock wave diffraction processes are given and discussed. For all cases, a uniform Cartesian grid system is employed, and the simulations of shock wave diffraction by a square cylinder are given for Bose–Einstein, Fermi–Dirac and Maxwell–Boltzmann gases at various relaxation times. Although most results are in the near continuum regime (i.e. small relaxation time), results with finite relaxation times corresponding to the slip regime and transitional regime are also included. For the equilibrium Euler limit solution, we expect that the diffraction patterns are similar to existing experimental results and other computational results based on classical Euler equations. At t=0, z1,ux1,uy1 and T1 are assigned. Then, with the given γ value, we may calculate n1, P1 and ϵ1, respectively, through 5.1 5.2 and 5.3 By assigning the Ms number and using equations (4.1)–(4.3) we can acquire n2, u2, P2 and ϵ2. Hence, z2 can then be calculated by solving equation (3.14). In the examples discussed in this paper, we initially set a plane shock wave at x=0.4. The computational domain is constructed on a Cartesian coordinate at 0≤x≤10;0≤y≤10. A square cylinder is placed in the middle of the computational domain at 4.5≤x≤5.5;4.5≤y≤5.5. The shock wave is set to propagate at Mach number Ms=2.0 towards the square cylinder and experiences transient shock diffraction. The mesh is constructed with Δx=1/400 and Δy=1/400.
For the specular reflection wall boundary conditions, the numerical computation of the information at the wall is detailed here. At grid points around the boundary of the square cylinder facing the incident shock, the information required to analyse the total fluxes are incomplete. Suppose i=K+1 denote the grid point immediately inside the wall, in which the information of is not known a priori, then one would lack one point of information needed to calculate . For this, we insert the ghost cells relevant to the assumed boundary condition. For example, for a specular type of reflection, we insert the ghost cells . Hence, the flux can be calculated through 5.4 If the boundary conditions do not contain the temperature of the wall, i.e. χ=0, then the distribution function with any temperature satisfies equation (4.7) (Cercignani 1995). The same treatment can be applied at the right, bottom and top boundaries of the square cylinder. At the corners of the square cylinder, multi-valued ghost cells are assigned.
(a) Example 1: shock wave diffraction by a square cylinder: Fermi–Dirac gas
In this example, we set Ms=2.0 and the initial condition for the Fermi–Dirac gas at rest ahead of shock are set as (n1,ux1,uy1,P1)=(0.1497,0,0,0.038), which correspond to z1=0.1 and T1=0.5, whereas the conditions before the shock are set as (n2,ux2,uy2,P2)=(0.2994,0.7125,0,0.192), which correspond to z2=0.078. The relaxation time is set as τ=0.0005. In figure 1, the density contours of a Fermi–Dirac gas obtained using TVD scheme and van-Leer’s limiter with CFL=0.9 and 20×20 discrete velocity points are shown in a chronological order. The subsequent development of the diffraction process following the shock-wall collision can be seen and the transition from regular reflection to Mach reflection can be observed in detail. The primary incident shocks, reflected bow shocks, Mach shocks, contact discontinuities, triple points, and vortices can be clearly identified and the complicated flow interaction can be represented.
(b) Example 2: shock wave diffraction by a square cylinder: Maxwell–Boltzmann gas
For the case of a Maxwell–Boltzmann gas, the initial conditions are set as (n1,ux1,uy1,P1)=(0.157,0,0,0.0393), which correspond to z1=0.1 and T1=0.5. Assigning the same Mach number Ms=2.0 as in example 1, the conditions of state 2 are (n2,ux2,uy2,P2)=(0.314,0.7073,0,0.1964). The relaxation time is again set as τ=0.0005.
In figure 2, using 20×20 discrete velocity points, we can see the diffraction process in a chronological order. The complex flow field and interaction behind the square cylinder comprising primary incident shocks, reflected bow shocks, Mach shocks, contact discontinuities, triple points and vortices can be represented.
(c) Example 3: shock wave diffraction by a square cylinder: Bose–Einstein gas
For the case of a Bose–Einstein gas, the initial conditions are (n1,ux1,uy1,P1)=(0.1638,0,0,0.0395), which correspond to z1=0.1 and T1=0.5. Accordingly, behind the shock, the initial conditions are (n2,ux2,uy2,P2)=(0.3276,0.6944,0,0.1975). A constant relaxation time of τ=0.0005 is implemented for the moving shock wave of Ms=2.0. Like the other two statistics, the major flow structures of primary incident shocks, reflected bow shocks, Mach shocks, contact discontinuities, triple points and vortices can all be illustrated clearly in terms of density contours in figure 3.
To demonstrate further the present method and to assist a better analysis of the flow field, the comparison of pressure contours of the three statistics is given in figure 4. The legends (colour online) are given to illustrate the difference of pressure profiles between the three statistics depicting a Bose–Einstein gas with a qualitatively higher pressure region in front of the square cylinder followed by Maxwell–Boltzmann and Fermi–Dirac gases, respectively, at the same assigned values of z1 and Mach number. At the immediate cell in front of the square cylinder, the pressure values are P=0.315,0.312 and 0.298 for Bose–Einstein, Maxwell–Boltzmann and Fermi–Dirac. The enlarged views of fugacity contours around the square cylinder are given in figure 5. Note that along with Bose–Einstein and Fermi–Dirac statistics, with the current method, we also can generate information of fugacity for Maxwell–Boltzmann gas.
(d) Example 4: effects of different relaxation times
Lastly, we examine the effects of different values of relaxation time on the flow structures. Four cases were considered. We denote the results obtained by setting f=f(0) in the present method as the semiclassical Euler solution, which corresponds to τ=0. The case of τ=0.0005 should give results similar to those of the semiclassical Navier–Stokes solution; the case of τ=0.05 should give rarefied flow in the transitional regime. Finally, the case of τ=0.1 will exhibit even more rarefied flows. In figure 6, the density contours of Fermi–Dirac gases at different relaxation times are given for comparison and for demonstrating the robustness of the present method. The conditioning of the gas at the initial state is the same as given in example 1. Here, we compare the applications of τ=0.0005, τ=0.05, τ=0.1 and the Euler limit solution in which equilibrium is assumed at every time step. Physically, the Euler limit solution shall depict the condition where τ≈0. Hence, at the Euler limit solution, the sharp flow profiles of the incident shocks, reflected shocks, contact discontinuities and vortices are the most conspicuous. The flow structure at a relatively small relaxation time τ=0.0005 shows a structure closer to that of the Euler solution with slight appearance of diffusivity. This is due to the fact that, although not of a physical quantity, the application of relaxation time contributes to the addition of diffusivity to the flow field. Therefore, we can observe that the results of those with high relaxation times, τ=0.05 and τ=0.1, show very diffusive flow fields. At these τ(s), the flows can be associated to those in the transitional regimes.
6. Concluding remarks
In this work, the unsteady shock wave diffraction patterns by a square cylinder in gases of particles of arbitrary statistics are computed by using an accurate solver for the semiclassical Boltzmann–BGK equation in phase space. A series of diffraction patterns of different flow properties at various time instants are recorded for detailed flow fields analysis. The specular reflecting wall is assumed in all calculations. The numerical method is a combination of the discrete ordinate method and TVD scheme. Different ranges of relaxation times covering continuum and rarefied flow regimes are computed. For comparison purposes, the equilibrium limit solution, i.e. zero relaxation time and f=f(0) everywhere, is also computed. The effects of gas particles that obey the Maxwell–Boltzmann, Bose–Einstein and Fermi–Dirac statistics are examined and depicted. When the relaxation time values are small such that the flows are in the continuum regime, for a Maxwell–Boltzmann gas, all the expected flow features of gas dynamics comprising incident shocks, reflected shocks, triple Mach shocks, slip lines and their roll-up to form vortices can be seen with considerably good detail, and are in good agreement with available experimental observations. Moreover, even for gases with Maxwell–Boltzmann statistics, the present results are quite different from traditional classical results as the fugacity or chemical potential is present. The capability in describing the gas flows of arbitrary particle statistics in various flow regimes on equal footing have been illustrated without many major constraints. We believe the present work is the first result that gas dynamical characteristics of three particle statistics in an identical flow problem have been simultaneously exhibited. Owing to the intrinsic nature of the statistics, the flow characteristic of the Maxwell–Boltzmann particles always lies between those of Bose–Einstein and Fermi–Dirac particles. The present work emphasizes building the unified and parallel framework for treating semiclassical rarefied gas flows, and specific application to different individual carriers in a specific physical problem can be separately pursued, such as electrons transport (Carrillo et al. 2003b; Majorana & Pidatella 2002) and phonon energy transport (Hsieh & Yang 2010). More general wall boundary conditions such as the Maxwell-type boundary conditions, including partially specular and partially diffusive, can be directly employed.
This work is done under the auspices of National Science Council, TAIWAN through grants NSC-99-2922-I-606-002 and CQSE subproject no. 5 99R-80873.
- Received May 3, 2011.
- Accepted October 3, 2011.
- This journal is © 2011 The Royal Society