Fast spectral solution of the generalized Enskog equation for dense gases

Article history: Received 10 July 2015 Received in revised form 13 September 2015 Accepted 18 September 2015 Available online 28 September 2015


Introduction
When the gas molecular mean free path (MFP) is comparable to the characteristic flow length, the Navier-Stokes equations based on the continuum-fluid assumption fail, and the Boltzmann equation is needed to study the rarefied gas dynamics [1]. Since in most situations the MFP is much larger than the molecular dimension, the binary collision modeled by the Boltzmann collision operator is localized in space. However, the collision is delocalized when the gas is compressed to such a density that the ratio of the molecular dimension to the MFP becomes appreciable; yet the continuum-fluid assumption is still inappropriate because of the small characteristic flow system size. This situation occurs, for instance, in gas flow in ultra-tight shale strata. For these dense gases, the Boltzmann equation is extended to the Enskog equation for hard-sphere molecules: where f (t, x, v) is the velocity distribution function (VDF) of molecular velocity v at spatial location x and time t, and a is the external acceleration. The right-hand side of Eq. (1) is the Enskog collision operator (ECO), where σ is the molecular diameter, d is the dimension of the velocity space, k is a unit vector associated with the relative position of two molecules at the time of their impact, and u = v − v * is the relative pre-collision velocity of two colliding molecules. The post-collision velocities v and v * are related to the pre-collision molecular velocities v and v * through v = v − (u · k)k and v * = v * + (u · k)k. Finally, the nonlocal collision is characterized by the product of the VDF and a pair correlation function at different locations: The velocity-independent pair correlation function χ is a given function of gas density. In the standard Enskog theory, χ coincides with the pair correction function in a state of uniform equilibrium evaluated at the contact point, while in the revised Enskog theory it is evaluated as the local equilibrium value in a nonuniform state [2]. In this paper, the former is adopted and the detailed expressions for the pair correlation function will be given below. The ECO has been extended to dense granular gases (i.e. materials composed of many small discrete grains), where the kinetic energy is not conserved during collisions [3][4][5]. A restitution coefficient 0 ≤ α ≤ 1, which is a function of the normal relative pre-collision velocity |u · k|, is introduced to characterize the loss of kinetic energy. Neglecting rotational degrees of freedom, the generalized Enskog collision operator (GECO) becomes: where ṽ and ṽ * are the pre-collision velocities producing post-collision ones v and v * after collision, and J is the Jacobian of the transformation from (ṽ, The generalized Enskog equation has applications ranging from astrophysics (e.g. stellar cloud formation) to industrial processes (e.g. fluidized beds and transport lines, handling of pharmaceuticals, and shale gas exploitation). The Enskog-Vlasov equation with the mean-field potential is even able to describe both liquid and vapor phases, and automatically capture liquid-vapor and fluid-solid interfaces [6,7], which has potential applications to multiphase flows in micro/nanoelectromechanical systems.
While various numerical methods have been proposed for the Boltzmann equation, there are only a few for the Enskog equation, let alone for the generalized equation. The most difficult part is the numerical solution of the collision operator. For elastic collisions, the ECO was first solved using a Monte Carlo quadrature method [8], and the study of shock propagation in a dense hard-sphere gas indicated that the Enskog equation can produce exact solutions comparable to molecular dynamics (MD) simulations [9]. Later, inspired by direct simulation Monte Carlo (DSMC) methods for the Boltzmann equation [10,11], several particle schemes were proposed for the simulation of dense gases [12][13][14][15], as well as the MD-DSMC hybrid method [16]. For inelastic collisions, the GECO has been solved by the DSMC [17,18] and spectral [19] methods for spatially-homogeneous relaxation problems.
As is well-known, the DSMC method has a computational cost proportional to the number of simulated particles (which may be far fewer than the actual number of molecules). However, particle-based methods are vulnerable to noise when macroscopic flow quantities are much smaller than the corresponding characteristic values (for example, when the flow velocity is much smaller than the sound speed). On the other hand, the spectral method is a deterministic numerical approach which solves the GECO with a computational cost O (N 2d ), where N is the number of discretized velocities in each direction [19]. It is good at resolving small flow signals (e.g. it can easily capture the high-energy tail of the VDF in a heated granular gas [19]) but has greater computational complexity than the DSMC method. A deterministic numerical method with high accuracy but small computational cost would be in strong demand.
Recently, a fast spectral method (FSM) with a computational cost of O (N d log N) has been proposed for solving the Boltzmann equation with spectral accuracy [20][21][22]. It has been successfully applied to many challenging cases for which the DSMC method is extremely time-consuming [23][24][25]. In this paper, we investigate the applicability of the FSM for solving the generalized Enskog equation.
This paper is organized as follows. The weak form of the GECO is expressed in the Carleman representation in Section 2, and is approximated by the FSM in Section 3. In Section 4, the accuracy of the FSM is assessed in several numerical tests. In Section 5, Fourier and Couette flows of a dense granular gas are investigated. Finally, concluding remarks are made in Section 6.

Weak form of the GECO in the Carleman representation
When the restitution coefficient α is a constant, we have ṽ = v − (1 + α)/(2α)(u · k)k, ṽ * = v * + (1 + α)/(2α)(u · k)k, and the Jacobian is J = 1/α. When α is a function of |u · k|, it is hard to express ṽ and ṽ * as functions of v and v * , and the Jacobian is very complicated. Therefore, it is convenient to write the GECO (3) in its weak form which does not involve the Jacobian and the pre-collision velocities ṽ and ṽ * .
For all smooth functions (v), the following equation holds for the GECO (we omit the location variable in what follows) To develop a FSM for the approximation of the GECO, it is useful to express its weak form in the Carleman representation. To achieve this, a unit vector satisfying the relation is introduced, and the post-collision velocity becomes By the relation 2 d−1 u · kdk = |u| sin 3−d (θ/2)d , where θ = arccos(u · /|u|) is the angle between u and , Eq. (4) is simplified to With the basic identity where δ is Dirac's delta function, the gain part of the weak form of GECO is written in the Carleman representation: Note that in the above manipulations we have used the transformations y = (|u| − u)/2 and z = v * − v − y = −u − y. Therefore, from Eq. (5) we obtain k = −y/|y|. Since the delta function requires that y be perpendicular to z, the restitution coefficient α is a function of |y| only. Also, cos θ = · u/|u| = −(y − z) · (y + z)/| y + z| 2 = (|z| 2 − |y| 2 )/(|y| 2 + |z| 2 ), which leads to |u| sin(θ/2) = |y|. Hence the weak form of the GECO for hard-disks (d = 2) and hard-spheres (d = 3) is written as

Fast spectral method for the GECO
The VDF is discretized on the truncated velocity d and N k being the number of velocity grid points in the k-th velocity direction.
The Fourier coefficients of the VDF are given bŷ where j = ( j 1 , j 2 , · · · , j d ), i is the imaginary unit, and ξ j = jπ /L v are the frequency components.
The weak form of the GECO given by Eq. (8) is also truncated [19,20], with the radius of y and z being R and is the kernel mode, with e, e the vectors in the unit sphere having the same direction as y and z, respectively.
The kernel mode β(l, m) should be decomposed into the form β(l, m) = p A p (l)B p (m), so that Eq. (10) can be calculated by FFT-based convolution, with a computational cost O (pN d log N) and p N d / log N. To achieve this, the integral with respect to ρ in Eq. (12) is approximated by Gauss-Legendre quadrature. Suppose ρ r and ω r (r = 1, 2, · · · , M r ) are the abscissas and weights of the Gauss-Legendre quadrature for ρ in the where The kernel mode can be further simplified. When d = 2, the two unit vectors e 1 and e 2 = −e 1 are perpendicular to e.
The integral with respect to e in Eq. (10) can be approximated by the trapezoidal rule. We divide the polar angle of the unit vector e into M sections, i.e. θ p = 2π /M with p = 1, 2, · · · , M. Then Eq. (10) is approximated by where e p = (cos θ p , sin θ p ). When d = 3, following the steps from Eq. (34) to Eq. (37) in Ref. [22] gives the final expression for the kernel mode as with J 0 and J 1 being the zeroth-and first-order Bessel functions of first kind, respectively. The integral with respect to e in Eq. (10) is again approximated by the trapezoidal rule. We divide the polar angle of the unit vector e into M 1 sections and the azimuthal angle of the unit vector e into M 2 sections, i.e. θ p = π/M 1 and ϕ q = 2π /M 2 with p = 1, 2, · · · , M 1 and q = 1, 2, · · · , M 2 . Then Eq. (10) becomes where e pq = (sin θ p cos ϕ q , sin θ p sin ϕ q , cos θ p ).
When the Fourier coefficients Q(x) have been obtained, the GECO Q(x) can then be calculated through the following FFT: According to the analysis above, the calculation of the GECO consists of three main steps. The first is to generate the spectrum of the VDF according to Eq. (9). The second step is the calculation of Q(x) according to Eqs. (17) or (19). The final step is to find the GECO in velocity space using Eq. (20). While the first and third steps have a computational cost O (N d log N), the calculation of Q(x) is time-consuming. If Eqs. (17) or (19) are solved by direct summation, the computational cost is O (N 2d ). However, the cost is reduced when Q(x) is calculated through FFT-based convolution, to a computational cost For elastic collisions, we have α = 1, and Eq. (15) does not depend on ρ r . Therefore, we introduce the function to simplify Eq. (13) to , which is N times faster than the inelastic collision case.
Finally, let us remark that although here we have only considered the hard-sphere potential between colliding molecules, the method applies also to general intermolecular potentials. In that case, the collision kernel k · u in Eq. (3) is replaced by general forms, say b(|u|, k · u), which, under the Carleman representation, is a function of |y| and |z| denoted by B(| y|, |z|).
Compared with the FSM for the Boltzmann equation for hard-sphere gases, the major difference is that the spectrum ĝ has to be collected from different locations when the unit vector e rotates, see Eqs. (11) and (19). Surprisingly, for elastic collisions, the computational cost for the Enskog equation is the same as for the Boltzmann equation when solved by the FSM. This is very interesting, considering the fact that the Boltzmann equation can be solved by other spectral methods with a cost of O (N 2d ), while the Enskog equation cannot be solved by conventional spectral methods.

Numerical accuracy tests
For spatially-homogeneous problems, the generalized Enskog equation reduces to the Boltzmann equation, and for elastic collisions it has been proven analytically and numerically that the FSM conserves mass, while momentum and energy are conserved with spectral accuracy. We now assess the accuracy of the FSM by considering three basic problems. First, the spatially-homogeneous relaxation of heated granular gases. Then, Poiseuille and Fourier flows (where we compare our results with MD and DSMC simulations).

A heated spatially-homogeneous granular gas
Consider a heated spatially-homogeneous granular gas, where the VDF is governed by the following Boltzmann equation: In this equation, the energy gain coefficient ≥ 0, the pair correlation function χ is a constant, and the term C 0 |u| γ represents the variable hard-sphere collision kernel: γ = 0 and 1 are for pseudo-Maxwellian and hard-disk (hard-sphere) gases, respectively.
For a pseudo-Maxwellian gas, the temperature relaxation can be obtained explicitly when the restitution coefficient α is a constant. Let us consider an initial datum such that the molecular number density n = f dv = 1 and the macroscopic velocity U = v f dv/n = 0. The time evolution of the temperature T = (2/d) |v − U | 2 f dv is governed by the following ordinary differential equation: where the first term on the right-hand side is the diffusion energy gain, while the second one is the collisional energy loss. Another exact expression for the kurtosis has also been found [26]; in a two-dimensional pseudo-Maxwellian gas it is . (25) Equations (24) and (25) can be used to assess the accuracy of the FSM. We first calculate the collisional loss of energy when the VDF takes the form f = v 2 exp(−v 2 )/π , for d = 2 and α = 0.5. In our numerical simulations, we take M r , M = 12.
When the number of velocity grid points increases from 16 2 to 24 2 , and finally to 32 2 , the relative error between the numerical and analytical results for the collisional energy loss decreases from 0.01 to 3.5 × 10 −4 , and finally to 2.8 × 10 −5 .
This demonstrates the spectral accuracy of the FSM. We also compare the time-relaxation of the temperature with the analytical solution, and we see in Fig. 1(a) that the FSM has a high accuracy. Fig. 1(b) shows that the kurtosis obtained from the FSM (solid line) compares well with the analytical expression (circles).
We also investigate the asymptotic behavior of the steady-state solution. Theoretically, the high-energy tail of the VDF in Eq. (23) is expected to behave as [18,[27][28][29][30][31][32][33] f (v, t = ∞) ∼ exp −A|v| 1+ γ 2 , |v| → ∞ (26) where A is a velocity-independent constant depending on α and γ . As shown in Fig. 2, these asymptotic high-energy tails are successfully captured in our FSM, even for the soft-potential gases for which γ is negative (i.e. molecules with smaller relative velocity have larger collision frequency). Since in DSMC the simulated particle population at high velocities is very small, spectral methods have the advantage of resolving these high-energy tails efficiently.

Force-driven Poiseuille flow of a hard-disk gas in a microchannel
Consider the flow of a large number of hard disks in a two-dimensional channel, subject to an external acceleration a in the x 1 direction. The two infinite plates of temperature T w are located at x 2 = −σ /2 and x 2 = L + σ /2, respectively. When a hard-disk hits a plate, it is diffusely reflected. For example, the boundary condition for the VDF at x 2 = 0 is where n w = −2 πm 0 /2k B T w v 2 <0 v 2 f (x 2 = 0, v)dv and m 0 is the mass of the hard-disk.
The acceleration is small so the flow is laminar, and we are interested in the velocity and temperature profiles when steady-state is reached. These profiles are determined by the reduced density calculated at the average number density n = n 0 , i.e. by the Knudsen number, and by the dimensionless acceleration, where k B is the Boltzmann constant, λ = [2 √ 2n 0 σ χ(n 0 )] −1 is the equilibrium MFP of a hard-disk, and is the pair correlation function for hard-disk systems [34]. Note that, although many other expressions for the paircorrelation function have been proposed [35,36], Eq. (31) is accurate up to the random packing limit η 2 0.74.
In our numerical simulations, the spatial region is discretized by 201 grid points, while the velocity space is discretized by 32 × 64 points; the polar angle of the unit vector e is divided into M = 12 sections. In the calculation of the ECO, as e rotates, the pair correlation function χ (x + σ e/2) and the spectrum f (x + σ e) in Eq. (11) may not fall exactly on the spatial grid points. In this case, a linear interpolation is employed, and when the location x + σ e is out of the computational domain, ĝ(x, e) is zero. The Enskog equation (1) with d = 2 is solved by an iterative numerical method, see Eq. (54) in Ref. [22]. The flow velocity and temperature are calculated as The iterations are terminated when the maximum relative difference in velocity and temperature between two consecutive steps is less than 10 −5 . Our numerical method is very efficient, and steady-state results are obtained in less than 5 minutes from our Matlab programme running on a PC with a single Intel Xeon 3.3 GHz CPU.
Comparable results of event-driven MD simulations are given in Ref. [37]. Since the hard-disk can only occupy regions half a diameter away from the walls, the reduced density η 2 and the dimensional acceleration Fr in Ref. [37] are multiplied by W /(W − 1) and (W − 1)/W , respectively, where W is the ratio between the wall distance and the molecular diameter.
Also, the Knudsen number is adjusted accordingly.
Comparisons of velocity and temperature profiles for various values of η 2 , Kn, and Fr between numerical results from the FSM and the published MD simulations are shown in Fig. 3. The agreement is satisfactory: the maximum relative difference in velocity is within 1.5%, while that in temperature is even smaller. Since the Enskog equation can be solved efficiently by the FSM, it can be used to investigate two-dimensional hard-disk flows; the generalized Enskog equation for dense granular gases will be investigated in Section 5 below.

Heat transfer in a hard-sphere gas
Consider the heat transfer through a dense hard-sphere gas between two parallel plates. The geometry is the same as in the force-driven Poiseuille flow case, but the plate at x 2 = L + σ /2 has a temperature of 2T w . In order to compare with the DSMC results in Ref. [38], we set L = 20λ, where the MFP is λ = [ √ 2πn 0 σ 2 χ (n 0 )] −1 , and the pair correlation function and reduced density η 3 are defined as [34,39] Note that Eq. (33) is accurate up to η 3 0.47, at which point the hard-sphere fluid undergoes a liquid to solid phase transition. Also, compared to the pair correlation functions in Refs. [35,36], χ is only 12.5% smaller at the random packing limit η 3 0.64.
In our numerical simulations, the three-dimensional velocity space is discretized into 32 × 48 × 32 grid points, while the spatial region is divided into 200 sections. In the discretization of the solid angle, we use M 1 = 8 and M 2 = 12. The Enskog equation (1) with d = 3 is solved by an iterative numerical method, see Eq. (55) in Ref. [22]. The iterations are terminated when the maximum relative difference in macroscopic flow quantities between two consecutive steps is less than 10 −5 . These macroscopic flow quantities include the molecular number density n, the temperature T , the kinetic pressure P kin 22 = m 0 v 2 2 f dv, and the kinetic heat flux Q kin Two values of the reduced density calculated at the average molecular number density, i.e. η 3 (n 0 ) = 0.1 (and 0.2) are considered, where the MFP is approximately equal to (and one-third of) the molecular diameter. The comparison of the results from the FSM and from DSMC is shown in Fig. 4, and good agreement can be seen. Note that when η 3 (n 0 ) = 0.2, the local value of η 3 (n) near the cooler plate reaches 0.6; when using the pair correlation function from Ref. [36], nothing changes in the macroscopic flow profiles.

Dense granular gases in spatially-inhomogeneous problems
Based on the excellent performance of the FSM for the Enskog equation, we now investigate dense granular gas flows that can be described by the generalized Enskog equation. We focus on the influence of the restitution coefficient on the macroscopic and microscopic properties in Fourier and Couette flows between two parallel plates. Note that in the dilute limit η → 0, the influence of inelasticity on the hydrodynamic profiles in the Couette and Fourier flows has already been extensively studied [40][41][42].

Fourier flow of a dense granular gas
The geometry is the same as the Fourier flow in Section 4.3, but the collisions between the hard-disks are inelastic. We set the channel width to be 20 times the disk diameter, and the reduced density to be η 2 (n 0 ) = 0.1. Therefore, the Knudsen number defined by Eq. (29) is 0.118, and the disk diameter is 0.422 times the MFP. The Maxwell diffuse boundary condition is used to describe interactions between the hard-disks and the two plates.  sions, the density in the bulk region (i.e. about one disk diameter away from the two plates) decreases monotonically from left to right, while the temperature varies oppositely. However, when α is below 0.9, the granular gas temperature becomes a convex function of x 2 , while the density in the bulk region is concave. As α decreases, the temperature becomes lower, while the hard disks increasingly concentrate around the central region of the channel. In the limiting case of α = 0, the  peak hard-disk density is about 5 times larger than the density near the warmer plate, while the minimum temperature (near the peak density) is about 17 times smaller than that of the cold plate. The concave density profile is probably because the granular gas has to be less dense near the plates in order to minimize the kinetic energy loss due to collisions, so that a heat flow path can be established across the channel. One of the most interesting questions is the asymptotic high-energy tail of the steady-state VDF. In Section 4.1 we showed that in the spatially-homogeneous relaxation of a heated granular gas the asymptotic behavior is determined by the collision kernel, irrespective of the restitution coefficient, see Eq. (26). For hard disks with γ = 1, we have log( f ) ∼ −|v| 3/2 as |v| → ∞. However, our numerical results in Fig. 6 show that the high-energy tail is strongly affected by the restitution coefficient. When the binary collision is elastic, we have log( f ) ∼ − v 2 throughout the channel (not shown). When the binary collision is inelastic, although the high-energy tail of the VDF in the middle of the channel can be fitted well by power-law functions, the exponent decreases from 2 as α decreases. In the limiting case α = 0, the high-energy behavior of the hard-disk flow is like that of the soft-potential gas in the heated spatially-homogeneous granular gas case (see Fig. 2).

Couette flow of a dense granular gas
We also consider Couette flow between two plates of temperature T w . The plate at x 2 = −σ /2 moves in the x 1 direction with speed 2k B T w /m 0 , while the plate at x 2 = L + σ /2 moves in the opposite direction with the same speed. We also assume L = 20σ and η 2 (n 0 ) = 0.1.
The density, temperature, and velocity profiles for various values of the restitution coefficient α are shown in Fig. 7. For elastic collisions, the temperature is higher than the wall temperature because of viscous heating. However, dissipation due to the kinetic energy loss during binary collisions quickly takes over when α ≤ 0.9, and the gas temperature is below the wall temperature: the smaller the α, the lower the temperature. The density variation is similar to the Fourier flow case: it is a concave function of x 2 in the bulk region when α ≤ 0.9. The flow velocity is a monotonically-decreasing function of x 2 , and the flow velocity decreases with α. The decrease of the flow velocity at the plate x 2 = 0 is because, as α decreases, the density becomes lower, so the local Knudsen number near the plates becomes larger, which results in higher velocity slip near the plates. The high-energy tail of the steady-state VDF is shown in Fig. 8. As in the Fourier flow case, the restitution coefficient has a remarkable influence on the asymptotic behavior.
Finally, we investigate the case when the average density of the hard-disk flow is high, i.e. η 2 (n 0 ) = 0.5, so that the hard-disk diameter is 5.40 times the MFP, and Kn = 0.0093 when L = 20σ . The behaviors of the macroscopic quantities shown in Fig. 9 are similar to those in the lower density case of η 2 (n 0 ) = 0.1, but their changes are more sensitive to the restitution coefficient because the dissipation is stronger at higher densities. Fig. 7. Density, temperature, and velocity profiles in the two-dimensional Couette flow of hard disks, when the reduced density is η 2 (n 0 ) = 0.1. Due to symmetry, the profiles for 10 ≤ x 2 /σ ≤ 20 are not shown.   9. Density, temperature, and velocity profiles in the two-dimensional Couette flow of hard disks, when the reduced density is η 2 (n 0 ) = 0.5.

Conclusions
Based on the Carleman representation of the collision operator in its weak form, and a Fourier-Galerkin discretization of velocity space, we have developed a fast spectral method to solve the generalized Enskog collision operator for dense granular gases. This method has a computational cost of O (M d−1 N d+1 log N), and works for any form of the intermolecular potential. In particular, the computational cost reduces to O (M d−1 N d log N) for elastic collisions with special collision kernels such as those of hard-disk and hard-sphere gases. The accuracy of the proposed method has been demonstrated through several numerical tests. To our knowledge, this is the first time that the generalized Enskog collision operator has been solved by a deterministic numerical method with both high accuracy and high efficiency.
The Fourier and Couette flows of dense granular gases have also been investigated. We found that the restitution coefficient has a strong influence on the molecular density and high-energy tail of the velocity distribution functions. Generally speaking, the smaller the restitution coefficient, the more molecules concentrate around the central region of the channel and the more molecules there are with high velocities.
The fast spectral method can be extended straightforwardly to gas mixtures, which could then help reveal the abundant exotic phenomena in dense granular gases. For granular flows with small values of the restitution coefficient, the temperature decreases drastically from the energy source due to the large dissipation, and thus a large number of discretized velocities is needed. In this case, velocity rescaling may be incorporated to improve the computational efficiency [43].