2.1. Boltzmann–BGK equation

Since the scattering wavelength (frequency) could be comparable to the mean free path (collision frequency) of gas molecules, one should resort to the gas kinetic theory to describe the density fluctuation. The state of the gas is described by the one-particle distribution function, $f(\boldsymbol {r}, \boldsymbol {\xi }, t)$, in the phase space of $(\boldsymbol {r}, \boldsymbol {\xi })$, where $\boldsymbol {r}$ and $\boldsymbol {\xi }$ are, respectively, the location and molecular velocity, and $t$ is the time. In the present work, the following Boltzmann equation with the BGK model (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954) is used to describe the rarefied gas dynamics:

(2.5)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol{\xi}\boldsymbol{\cdot}\boldsymbol{\nabla} f = \varOmega \equiv -\frac 1{\tau}[\,f-f^{(eq)}], \end{equation}

where $\tau$ is the relaxation time, $f^{(eq)}$ is the local Maxwellian

(2.6)\begin{equation} f^{(eq)} = \frac{\rho}{(2{\rm \pi} R T)^{D/2}} \exp\left(-\frac{|\boldsymbol{\xi} - \boldsymbol{u}|^2}{2RT}\right), \end{equation}

and $D$, $\rho$, $\boldsymbol {u}$ and $T$ are the dimension of the space, macroscopic mass density, velocity and absolute temperature, respectively. According to the Chapman–Enskog expansion (Huang Reference Huang1987), the dynamic shear viscosity of the gas can be found as $\mu = \tau \rho RT = \tau p$, where $p=\rho RT$ is the hydrostatic pressure.

Non-dimensionalizing (Shan et al. Reference Shan, Yuan and Chen2006) using the characteristic velocity $u_0 \equiv \sqrt {RT_0}$ and the characteristic time $t_0 \equiv L_c/u_0$, the Maxwellian, after normalizing by $\rho _0(RT_0)^{-D/2}$, with $\rho _0$ an arbitrary reference number density, becomes

(2.7)\begin{equation} f^{(eq)} = \frac{\rho}{(2{\rm \pi}\theta)^{D/2}} \exp\left(-\frac{|\boldsymbol{\xi} - \boldsymbol{u}|^2}{2\theta}\right), \end{equation}

with $\theta = T/T_0$. Equation (2.5) remains unchanged except that all quantities are dimensionless with $x$ scaled by $L_c$, $\boldsymbol {\xi }$ and $\boldsymbol {u}$ scaled by $u_0$, and $t$ and $\tau$ scaled by $t_0$. In particular, the dimensionless relaxation time is exactly the Knudsen number,

(2.8)\begin{equation} \tau = \frac{\mu}{pt_0} = \frac{\mu k\sqrt{RT}}p = {\textit{Kn}}. \end{equation}

Hence, the dimensionless Boltzmann–BGK equation in one dimension can be written as

(2.9)\begin{equation} \frac{\partial f}{\partial t} + \xi_x\frac{\partial f}{\partial x} = -\frac 1 {\textit{Kn}} [\,f-f^{(eq)}], \end{equation}

where $f^{(eq)}$ is given by (2.7). To complete the equations, we note that the mass density $\rho$, fluid velocity $u_x$ and internal energy density $e$ are all velocity moments of the distribution function as follows:

(2.10a)\begin{gather} \rho = \int f \, \textrm{d} \boldsymbol{\xi}, \end{gather}

(2.10b)\begin{gather}\rho u_x = \int f\xi_x \, \textrm{d} \boldsymbol{\xi}, \end{gather}

(2.10c)\begin{gather}\rho e = \tfrac{1}{2}\int f(\xi-u)^2 \, \textrm{d} \boldsymbol{\xi}, \end{gather}

where, by the energy equipartition principle, $e$ is related to the temperature $\theta$ by $e = D\theta /2$.

2.2. The free-molecular limit

An analytical solution can be obtained in the free-molecular limit, where collisions are rare and hence $\varOmega =0$. Therefore, (2.9) is reduced to the following linear equation:

(2.11)\begin{equation} \frac{\partial f}{\partial t} + \xi_x\frac{\partial f}{\partial x} = 0, \end{equation}

an obvious solution of which is

(2.12)\begin{equation} f(x, \boldsymbol{\xi}, t) = f(x-\xi_x t, \boldsymbol{\xi}, 0). \end{equation}

Using the initial condition that the distribution at $t=0$ is a local Maxwellian with $\theta =1$ and $\rho =1+\varepsilon \delta (x)$, we have

(2.13)\begin{equation} f(x, \boldsymbol{\xi}, t) = \frac{1+\varepsilon\delta(x-\xi_x t)}{\sqrt{2{\rm \pi}}^3} \exp\left(-\frac{\xi^2}{2}\right), \end{equation}

where we let $D=3$. Then the density is

(2.14)\begin{equation} \rho(x,t) = 1 + \frac{\varepsilon}{|t|\sqrt{2{\rm \pi}}} \exp\left(-\frac{x^2}{2t^2}\right). \end{equation}

On substituting (2.14) into (2.1) and performing some standard manipulations, we have

(2.15)\begin{equation} G(x,t) = \frac{\varepsilon^2}{|t|\sqrt{2{\rm \pi}}}\exp{\left(-\frac{x^2}{2t^2}\right)}. \end{equation}

The spectral density in the free-molecular limit can be calculated as

(2.16)\begin{equation} S = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{\varepsilon^2}{|t|\sqrt{2{\rm \pi}}} \exp\left(-\frac{x^2}{2t^2}\right) \exp({\textrm{i}(x-\omega t)})\,\textrm{d} x\,\textrm{d} t = \varepsilon^2\sqrt{2{\rm \pi}}\exp\left(-\frac{\omega^2}{2}\right). \end{equation}

It should be emphasized that $k=1$ is only considered after normalization with itself in the present work. Therefore, the spectrum in the free-molecular limit is

(2.17)\begin{equation} S({\textit{Kn}},\omega) = \varepsilon^2\sqrt{2{\rm \pi}}\exp\left(-\frac{\omega^2}{2}\right). \end{equation}

2.3. The continuum limit

In the continuum flow regime where ${\textit {Kn}}\ll 1$, the classical Navier–Stokes–Fourier (NSF) equations are valid. To consider the small perturbation on a stationary base fluid, we write

(2.18a–c)\begin{equation} \rho = \rho_0 + \rho',\quad u_x = u_x',\quad \theta = \theta_0 + \theta', \end{equation}

where the subscript zero denotes the properties of the base flow and the prime denotes that of the perturbation. Assuming that the flow is 1-D and all spatial derivatives except $\partial /\partial x$ vanish, the linearized forms of the NSF equations are

(2.19a)\begin{gather} \frac{\partial \rho'}{\partial t} + \frac{\partial u_{x}'}{\partial x} = 0, \end{gather}

(2.19b)\begin{gather}\frac{\partial u_{x}'}{\partial t} + \frac{\partial \rho'}{\partial x} + \frac{\partial \theta'}{\partial x} - \frac{2(D-1){\textit{Kn}}}D\frac{\partial^2 u_{x}'}{\partial x^2} = 0, \end{gather}

(2.19c)\begin{gather}\frac D2\frac{\partial \theta'}{\partial t} + \frac{\partial u_{x}'}{\partial x} - \frac{(D+2){\textit{Kn}}}{2{\textit{Pr}}}\frac{\partial^2 \theta'}{\partial x^2} = 0. \end{gather}

Hereinafter we set $D = 3$ to give the specific heat ratio of a monatomic gas.

On applying the following spatial and temporal Fourier transform to the quantities $M = [\rho ', u'_x, \theta ']^\textrm {T}$,

(2.20a,b)\begin{equation} M_k = \int_{-\infty}^{\infty} M\,\textrm{e}^{\textrm{i}x} \, \textrm{d} x, \quad \hat{M} = \int_{0}^{\infty} M_k \,\textrm{e}^{-\textrm{i}\omega t} \, \textrm{d} t, \end{equation}

where the integral limit of the temporal Fourier transform is from 0 to $+\infty$ because $G(x,t)$ is even in $t$, and noting that

(2.21)\begin{equation} \int_0^\infty \frac{\partial M_k}{\partial t} \, \textrm{e}^{-\textrm{i}\omega t} \, \textrm{d} t = -M_k(t=0) + \textrm{i}\omega\hat{M}, \end{equation}

where the initial condition has been considered in the above equation, equations (2.19) can be rewritten in the following matrix form (Wu & Gu Reference Wu and Gu2020):

(2.22)\begin{equation} \begin{bmatrix} \textrm{i}\omega & -\textrm{i} & 0\\ -\textrm{i} & \textrm{i}\omega + \dfrac{4}{3}{\textit{Kn}} & -\textrm{i}\\ 0 & -\textrm{i} & \dfrac{3}{2}\textrm{i}\omega + \dfrac{5}{2{\textit{Pr}}}{\textit{Kn}}\end{bmatrix} \begin{bmatrix} \hat{\rho}' \\ \hat{u}_x' \\ \hat{\theta}' \end{bmatrix} =\begin{bmatrix} \varepsilon \\ 0 \\ 0 \end{bmatrix}. \end{equation}

From this, the SRBS spectrum predicted by the NSF equations can be obtained as

(2.23) \begin{equation} S(\textit{Kn},\omega) = \varepsilon^2\mathrm{Re}\left\{\frac{6\textit{Kn}(5 + 4\textit{Pr})\omega - 2{\rm i}(20\textit{Kn}^2 + 6\textit{Pr} - 9\textit{Pr}\,\omega^2)}{5(4\textit{Kn}^2 + 3\textit{Pr})\omega - 9\textit{Pr}\,\omega^3 -3{\rm i}\textit{Kn} [5 - (5 + 4\textit{Pr})\omega^2]}\right\}. \end{equation}

It should be noted that the expression above is only valid when ${\textit {Kn}}\ll 1$.

On the other hand, by multiplying the numerator and denominator in (2.23) by the complex conjugate of the denominator, neglecting the lower- and higher-order terms with respect to ${\textit {Kn}}$ and $\omega$, respectively, and taking the real part, the asymptotic solution of (2.23) at ${\textit {Kn}}\gg 1$ and $\omega \leq 1$ can be obtained:

(2.24)\begin{equation} S({\textit{Kn}},\omega) = \frac{24\varepsilon^2{\textit{Kn}}}{9+16{\textit{Kn}}^2\omega^2}. \end{equation}

Obviously, the dependence of $S$ on $\omega$ does not tend to be Gaussian function dictated by the SRBS spectrum (2.17) in the free-molecular limit. This difference is due to the failure of the NSF equations to describe flows with large Knudsen numbers.

Below, in § 4.2, the spatial Fourier transform of the perturbation density $\rho '_{k}(t)$ will be used to extract the kinematic viscosity $\nu$ and thermal diffusivity $\kappa$. Note that $\kappa$ is related to the thermal conductivity $\varGamma$ by $\kappa = \varGamma /(c_p\rho)$, where $c_p$ is the heat capacity at constant pressure. It is obvious that $\rho '_{k}(t)$ can be derived by use of the inverse Fourier transform:

(2.25)\begin{equation} \rho'_{k}(t) = \frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty} \hat{\rho}'(\omega)\,\textrm{e}^{\textrm{i}\omega t}\,\textrm{d} \omega. \end{equation}

It should be noted that the denominator of the right-hand side of (2.23) is the characteristic polynomial of the linear hydrodynamic equations (Shan & Chen Reference Shan and Chen2007). The corresponding characteristic equation has three roots: one real root, $\omega _t$, that gives the thermal diffusion mode; and a pair of complex ones, $\omega _{\pm }$, giving the two acoustic modes (Shan Reference Shan2019). At the small-${\textit {Kn}}$ limit, one can write the dispersion relations of the three modes as

(2.26a)\begin{gather} \omega_t = \textrm{i} \frac{Kn}{Pr}, \end{gather}

(2.26b)\begin{gather}\omega_{\pm} = \pm\sqrt{\frac{5}{3}} + \textrm{i}\frac{(1 + 2{\textit{Pr}}){\textit{Kn}}}{2{\textit{Pr}}}. \end{gather}

Therefore, the denominator of the right-hand side of the (2.23) can be expressed as

(2.27)\begin{equation} -9{\textit{Pr}}(\omega - \omega_t)(\omega - \omega_+)(\omega - \omega_-). \end{equation}

By means of the residue theorem and dropping some high-order terms in ${\textit {Kn}}$, $\rho '_{k}(t)$ can be obtained in an analytical form when the Knudsen number is small:

(2.28)\begin{equation} \rho'_{k}(t) = \frac{2\varepsilon}{5}\exp\left(-\frac{{\textit{Kn}}}{{\textit{Pr}}}t\right) + \frac{3\varepsilon}{5}\exp\left[-\frac{(1+2{\textit{Pr}}){\textit{Kn}}}{3{\textit{Pr}}}t\right] \cos\left(\sqrt{\frac{5}{3}}t\right). \end{equation}

It can be seen that the initial disturbance in the density can propagate but will be eventually damped out by the viscous and thermal dissipative processes. In fact, the non-dimensional kinematic viscosity $\nu$ and thermal diffusivity $\kappa$ are, respectively,

(2.29a,b)\begin{equation} \nu = {\textit{Kn}}\quad\textrm{and} \quad\kappa = {\textit{Kn}}/{\textit{Pr}}. \end{equation}

Figure 1 shows the typical SRBS spectra (2.23) as predicted by the NSF equations. When ${\textit {Kn}}\ll 1$, the spectrum shows three peaks (the third one is not shown due to the symmetry about $\omega =0$) corresponding to different mechanisms (Reichl Reference Reichl2016). The Rayleigh peak centred at $\omega =0$ is due to the thermal mode of gas molecules. The other two Brillouin peaks centred at $\omega = \pm c_s = \pm \sqrt {5/3}$ are related to the sound mode. As ${\textit {Kn}}$ increases, the increased dissipation leads to larger width of all peaks. For comparison, the spectrum (2.17) in the free-molecular limit is also shown. As ${\textit {Kn}}\rightarrow \infty$, the NSF spectrum does not converge to it due to the failure of the continuum hypothesis.

Figure 1.

The normalized SRBS spectra calculated by the NSF equations for ${\textit {Kn}}=0.01$, $0.1$, $0.3$, $0.5$, $1.0$ and $5.0$, and ${\textit {Pr}}=1$. The horizontal and vertical axes are the angular frequency and $S(\omega )/S(0)$, respectively. The diamonds represent the spectrum (2.17) in the free-molecular flow regime. Only half of the spectrum is shown due to the symmetry with respect to $\omega = 0$.

3.1. Reduction to one dimension

We examine the capability of the LBM to capture the rarefaction effects by solving (2.9) numerically, where the molecular velocity space is three-dimensional. The reduced velocity distribution functions are introduced to save computational cost (Chu Reference Chu1965; Nie, Shan & Chen Reference Nie, Shan and Chen2008b; Sharipov & Kalempa Reference Sharipov and Kalempa2008). Specifically, we solve the following two equations:

(3.1)\begin{equation} \frac{\partial \psi_{i}}{\partial t} + \xi_x\frac{\partial \psi_i}{\partial x} = -\frac 1 {\textit{Kn}} [\psi_i - h_i\psi^{eq}], \quad i = 1, 2, \end{equation}

where $h_1 = 1$, $h_2 = 2\theta$, $\psi _i$ are the reduced distributions

(3.2a)\begin{gather} \psi_{1}(x, \xi_x, t) = \int f \, \textrm{d} \xi_{y}\, \textrm{d} \xi_{z}, \end{gather}

(3.2b)\begin{gather}\psi_{2}(x, \xi_x, t) = \int f(\xi_{y}^{2} + \xi_{z}^{2})\, \textrm{d} \xi_y\, \textrm{d} \xi_z, \end{gather}

and $\psi ^{eq}$ is the Hermite expansion of the Maxwellian. According to Meng et al. (Reference Meng, Zhang and Shan2011), the following second-order expansion is sufficient for linear problems:

(3.3)\begin{equation} \psi^{eq} = \rho\omega(\xi_x) [1 + u_x\xi_x + \tfrac 12(u^2_x+\theta-1)(\xi^2_x-1)]. \end{equation}

In the equations above,

(3.4)\begin{equation} \left.\begin{gathered} \displaystyle\rho = \int\psi_1\, \textrm{d} \xi_x,\\ \displaystyle\rho u_x = \int\xi_x\psi_1\, \textrm{d} \xi_x,\\ \displaystyle\rho e = \dfrac{3\rho\theta}2 = \dfrac 12\left[\int(\xi_x-u_x)^2\psi_1\, \textrm{d} \xi_x + \int\psi_2\, \textrm{d} \xi_x\right].\end{gathered}\right\} \end{equation}

Hereinafter the subscript $x$ is dropped for simplicity.

Following Shan et al. (Reference Shan, Yuan and Chen2006), the above equations are discretized in velocity space. Let $w_j$ and $\xi _j$, $j = 1, \ldots , d$, be the weights and abscissas of a degree-$q$ quadrature rule such that the identity

(3.5)\begin{equation} \int \omega(\xi)p(\xi)\, \textrm{d} \xi = \sum_{j=1}^dw_jp(\xi_j) \end{equation}

holds for any polynomial $p(\xi )$ of an order not exceeding $q$. Defining

(3.6)\begin{equation} \psi_{ij}(x,t) \equiv \frac{w_j\psi_i(\xi_j)}{\omega(\xi_j)}, \quad i = 1, 2, \quad j = 1,\ldots, d, \end{equation}

we have the following 1-D lattice BGK equations:

(3.7)\begin{equation} \frac{\partial \psi_{ij}}{\partial t} + \xi_j\frac{\partial \psi_{ij}}{\partial x} = -\frac 1{\textit{Kn}}[\psi_{ij} - h_i\psi^{eq}_j], \end{equation}

where

(3.8)\begin{equation} \psi^{eq}_j = \rho w_j [1 + u\xi_j + \tfrac 12(u^2+\theta-1)(\xi^2_j-1)]. \end{equation}

The macroscopic quantities are defined by

(3.9a–c)\begin{equation} \rho = \sum_{j=1}^d\psi_{1j},\quad \rho u = \sum_{j=1}^d \xi_j\psi_{1j},\quad 3\rho\theta+u^2 = \sum_{j=1}^d(\xi_j^2\psi_{1j} + \psi_{2j}). \end{equation}

As shown previously (Shan et al. Reference Shan, Yuan and Chen2006; Nie et al. Reference Nie, Shan and Chen2008a), the equations above preserve hydrodynamic moments up to the order of $q-N$ for a distribution truncated to the $N$th order. It is therefore desirable to use a quadrature rule that is as accurate as possible.

3.2. The spatial and temporal discretization

Equation (3.7) is a set of $2\times d$ partial differential equations in the 1-D space–time $(x, t)$, representing a significant simplification of (3.1) in two dimensions. The continuous space–time has to be discretized for the equations to be solved numerically. When the discrete velocities form a Bravais lattice in the velocity space (also known as ‘on-lattice’), the space–time discretization renders the Boltzmann equation a simple streaming-collision process similar to the d'Alembert solution to the wave equation. The fact that the advection of this process is exact greatly reduces numerical dissipation. However, as the abscissas of the extremely efficient Gauss quadrature hardly form a Bravais lattice except for the lowest orders, it is sometimes desirable to use them for the higher algebraic orders at the same number of collocation points. The latter approach (often referred to as ‘off-lattice’) involves some form of interpolation, which introduces additional numerical dissipation. In the present work, we use the finite-volume (FVM) lattice Boltzmann formulation (Kim et al. Reference Kim, Pitsch and Boyd2008; Meng & Zhang Reference Meng and Zhang2011a,Reference Meng and Zhangb; Ambruş & Sofonea Reference Ambruş and Sofonea2016a,Reference Ambruş and Sofoneab) with Gauss quadrature to benchmark its numerical performance against that of the on-lattice streaming-collision scheme. Below, we give the space–time discretization of both schemes.

To use an off-lattice set of discrete velocities, let the 1-D space be covered by cells of uniform size $\delta x$ and centred at $x_k$. Integrating equation (3.7) over the $k$th cell and denoting the cell averaging by

(3.10)\begin{equation} \bar{\psi}_{ij,k} = \frac{1}{\delta x}\int_{x_k-\delta x/2}^{x_k+\delta x/2}\psi_{ij}\, \textrm{d} x, \end{equation}

equation (3.7) can be cast into the FVM form

(3.11)\begin{equation} \frac{\partial \bar{\psi}_{ij,k}}{\partial t} + \frac{1}{\delta x}[F_{k+1/2} - F_{k-1/2}] = -\frac 1{\textit{Kn}}[\bar{\psi}_{ij,k} - h_i\bar{\psi}^{eq}_{j,k}], \end{equation}

where the fluxes of $\psi _{ij}$ through the left and right boundaries of the $k$th cell are

(3.12)\begin{equation} F_{k\pm 1/2} \equiv \xi_j\psi_{ij}(x_k\pm\delta x/2). \end{equation}

Using the forward Euler time stepping, we have

(3.13)\begin{equation} \bar{\psi}_{ij,k}^{t+\delta t} = \bar{\psi}_{ij,k}^{t} - \frac{\delta t}{\delta x} [F_{k+1/2}^t - F_{k-1/2}^t] - \frac{\delta t}{{\textit{Kn}}}[\bar{\psi}_{ij,k}^{t} - h_i^t\bar{\psi}_{j,k}^{eq,t}]. \end{equation}

Following Leveque (Reference Leveque1996), the fluxes are

(3.14)\begin{equation} \left. \begin{aligned} F_{k-1/2}^{+} &= \xi_j\bar{\psi}_{ij,k-1} + C(\xi_j)(\bar{\psi}_{ij,k} - \bar{\psi}_{ij,k-1}) \chi(\vartheta^{+}),\\ F_{k-1/2}^{-} &= \xi_j\bar{\psi}_{ij,k} + C(\xi_j)(\bar{\psi}_{ij,k} - \bar{\psi}_{ij,k-1}) \chi(\vartheta^{-}), \end{aligned}\right\} \end{equation}

where the superscripts $+$ and $-$ are the sign of $\xi _j$,

(3.15)\begin{equation} C(\xi) = \frac{|\xi|}2\left(1-\frac{|\xi|\delta t}{\delta x}\right) \end{equation}

and

(3.16a,b)\begin{equation} \vartheta^{+} = \frac{\bar{\psi}_{ij,k-1} - \bar{\psi}_{ij,k-2}}{\bar{\psi}_{ij,k} - \bar{\psi}_{ij,k-1}} \quad\textrm{and}\quad \vartheta^{-} = \frac{\bar{\psi}_{ij,k+1} - \bar{\psi}_{ij,k}}{\bar{\psi}_{ij,k} - \bar{\psi}_{ij,k-1}}, \end{equation}

with $\chi (\vartheta )$ being the flux limiter of the form

(3.17)\begin{equation} \chi(\vartheta) = \begin{cases} 0, & \vartheta<0, \\ 2\vartheta, & 0\leq\vartheta< \tfrac{1}{3}, \\ \tfrac12 (1+\vartheta), & \tfrac{1}{3}\leq\vartheta<3, \\ 2, & \vartheta\geq 3.\end{cases} \end{equation}

Since the convective and collision terms in (3.13) are discretized into explicit form, the time step in the FVM is restricted by both the Courant–Friedrichs–Lewy (CFL) condition and relaxation time:

(3.18)\begin{equation} \delta t = \min\left({\textit{Kn}},\frac{\alpha\delta x}{\xi_m}\right),\end{equation}

where $\alpha$ is the CFL number and $\xi _m$ is the maximum discrete velocity.

For an on-lattice set of discrete velocities, the streaming-collision algorithm can be implemented and the collision term is discretized in a semi-implicit form. Specifically, on integrating equation (3.7) from $t$ to $t+1$ and applying the trapezoidal rule for the collision term, we have

(3.19)\begin{equation} \psi_{ij}(x+\xi_j,t+1) - \psi_{ij}(x,t) = \tfrac{1}{2}[\varOmega_{ij}(x+\xi_j,t+1) + \varOmega_{ij}(x,t)]. \end{equation}

Clearly, the right-hand side is implicit. The implicity can be removed by the change of variable (He, Shan & Doolen Reference He, Shan and Doolen1998; Dellar Reference Dellar2001)

(3.20)\begin{equation} \varphi_{ij} = \psi_{ij} - \frac{\varOmega_{ij}}{2},\end{equation}

and (3.19) can then be rewritten as

(3.21)\begin{equation} \varphi_{ij}(x+\xi_j, t+1) - \varphi_{ij}(x,t) = \frac{1}{{\textit{Kn}} + \tfrac12}[h_i\varphi_{j}^{eq}(x,t) - \varphi_{ij}(x,t)], \end{equation}

where

(3.22)\begin{equation} \varphi^{eq}_j = \rho w_j [1 + u\xi_j + \tfrac 12(u^2+\theta-1)(\xi^2_j-1)]. \end{equation}

In fact, we directly track the distribution functions $\varphi _{1}$ and $\varphi _{2}$ instead of $\psi _{1}$ and $\psi _{2}$ in the on-lattice LBM.

3.3. The 1-D on-lattice quadrature

There have been several approaches to recover the well-known lattices and develop new lattices for the LBM (Philippi et al. Reference Philippi, Hegele, Dos Santos and Surmas2006; Chen & Shan Reference Chen and Shan2008; Chikatamarla & Karlin Reference Chikatamarla and Karlin2009; Surmas et al. Reference Surmas, Ortiz and Philippi2009; Karlin & Asinari Reference Karlin and Asinari2010; Yudistiawan et al. Reference Yudistiawan, Kwak, Patil and Ansumali2010). In this study, as discussed before, quadrature in the form of (3.5) with abscissas forming Bravais lattices in velocity space is often desirable. Here we apply the same procedure (Shan et al. Reference Shan, Yuan and Chen2006; Shan Reference Shan2010, Reference Shan2016) to a very high order in one dimension. In fact, for any given abscissas, one can always find a Hermite quadrature rule by solving the weights directly from (3.5) for a set of base functions of the functional space of polynomials. For instance, we can obtain the weights of a degree-$q$ rule by solving (3.5) for $p(\xi ) = 1, \xi , \xi ^2, \ldots , \xi ^q$, which yields ${q+1}$ linear equations and requires $q+1$ points to have a solution. An equivalent but much simpler way (Shan et al. Reference Shan, Yuan and Chen2006) is to use orthogonal polynomials, in this case the Hermite polynomials $H^{(n)}$, as the base functions. The sufficient and necessary condition for $w_i$ and $\xi _i$ to be the weights and abscissas of a degree-${q}$ quadrature rule is

(3.23)\begin{equation} \sum_{i=1}^{d}w_{i}H^{(n)}(\xi_i) = \int\omega(\xi)H^{(n)}(\xi)\, \textrm{d} \xi = \begin{cases} 1, & n=0, \\ 0, & n=1, 2, \ldots, {q}, \end{cases}\end{equation}

where the second equality is due to the orthogonality of $H^{(n)}$.

We note that, if $\xi _i$ can be freely chosen, by the celebrated Gauss integration theorem, the optimal choice is the set of $N$ roots of $H^{(N)}(\xi )$, which yields the GH rule of degree $q=2N-1$. The number of points required for degree $q$ is thus $(q+1)/2$, half of that for fixed abscissas. If the abscissas are required to be symmetric and form a Bravais lattice, the set of abscissas and weights can be written as $\{(\pm ck, w_k): k=0, 1, \ldots , m\}$. Equation (3.23) is automatically satisfied for all odd $n$, and we only need to solve $w_i$ and $c$ for even $n$ with the additional condition that all $w_i$ must be positive to ensure numerical stability.

The mathematical structure of the solution has been extensively discussed previously up to the ninth degree (Shan Reference Shan2016). Here we use the same procedure in one dimension up to $m=19$, corresponding to $q=2m+1=39$. Exploiting the degree of freedom associated with the lattice constant, $c$, one of the weights can be made zero to eliminate two of the $2m+1$ points, yielding degree-$q$ rules with $q-2$ points, three fewer than the $q+1$ points for fixed abscissas. It was also observed without proof that, for each given number of abscissas, $d=2m-1$, two solutions of $c$ exist. Shown in figure 2 are the lattice constants versus $d$. They are clearly segregated into two groups, of which the one with smaller $c$ will be referred to as P1 and the other as P2. Again, it should be emphasized that the P1 and P2 rules are constructed based on the belief of obtaining the highest-order on-lattice quadratures with the fewest abscissas. Actually, there are infinite sets of prescribed-degree on-lattice quadratures where all weights are greater than zero and the lattice constant is between the minimum and maximum lattice constants, indicated by P1 and P2 rules, respectively. This will be discussed in the future.

Figure 2.

Lattice constants of the 1-D on-lattice quadrature rules versus number of quadrature points, $d$. For each $d$, two rules with clearly separated lattice constants exist. Within each of the two groups, denoted as P1 and P2, the lattice constant decreases smoothly with $d$.

The lattice constant stretches the velocity space. To show the effective distribution of the abscissas in the velocity space and the magnitudes of their weights, we plot in figure 3 the weights of the two $d$-point rules, D1Q$d$P1 and D1Q$d$P2, against $\xi =ck$, the discretized microscopic velocity. Here $d = 17$, $27$ and $37$ for simplicity. Meanwhile, the weight of the $d$-point GH rule, GH$d$, is also plotted for comparison. Note that GH quadrature is of a much higher degree than the other two with a same number of abscissas. Rather surprisingly, the weights coincide with each other quite well. However, in the case of the same number of abscissas, the abscissas of P1 are more compactly distributed around zero velocity that those of the corresponding GH, which in turn is more compact than those of P2. We suspect that the relative compactness might explain the difference in numerical performance reported below in § 4.

Figure 3.

Abscissa distribution in the velocity space and the magnitudes of the associated weights for the two $d$-point on-lattice rules, D1Q$d$P1 and D1Q$d$P2, together with the $d$-point GH rule, GH$d$. Here $d=17$, $27$ and $37$. Only half of the discrete velocities are shown due to the symmetry about $\xi =0$. It is seen that these sets of weights almost collapse together, and the abscissas of P1 are more compactly distributed than those of the corresponding GH around zero velocity, which is in turn more compact than those of P2.

The full details of D1Q5 and D1Q37 are listed in tables 1 and 2, respectively.

Table 1.

Lattice constant and weights of the five-point on-lattice quadrature rules D1Q5P1 and D1Q5P2. The abscissas corresponding to the weight $w_k$ are $\pm ck$ for $k=0,1,2,3$.

Table 2.

Lattice constant and weights of the 37-point on-lattice quadrature rules D1Q37P1 and D1Q37P2. The abscissas corresponding to the weight $w_i$ are $\pm ck$ for $k=0,1, 2,\ldots ,19$.