Half-range lattice Boltzmann models for the simulation of Couette flow using the Shakhov collision term

TL;DR

This paper develops and validates half-range lattice Boltzmann models with the Shakhov collision term for simulating rarefied Couette flow, demonstrating accurate results for Knudsen numbers up to 0.1 and analyzing deviations at higher Kn.

Contribution

The paper introduces a novel mixed lattice Boltzmann approach using half-range and full-range quadratures with the Shakhov collision model for simulating rarefied gas flows.

Findings

Accurate simulation of Couette flow for Kn up to 0.1.

High quadrature orders needed for non-linear velocity profiles at Kn > 1.

Good agreement with DSMC results for Maxwell molecules, Helium, and Argon.

Abstract

The three-dimensional Couette flow between parallel plates is addressed using mixed lattice Boltzmann models which implement the half-range and the full-range Gauss-Hermite quadratures on the Cartesian axes perpendicular and parallel to the walls, respectively. The ability of our models to simulate rarefied flows are validated through comparison against previously reported results obtained using the linearized Boltzmann-BGK equation for values of the Knudsen number (Kn) up to $100$. We find that recovering the non-linear part of the velocity profile (i.e., its deviation from a linear function) at ${\rm Kn} \gtrsim 1$ requires high quadrature orders. We then employ the Shakhov model for the collision term to obtain macroscopic profiles for Maxwell molecules using the standard $\mu \sim T^\omega$ law, as well as for monatomic Helium and Argon gases, modeled through ab-initio potentials, where the viscosity is recovered using the Sutherland model. We validate our implementation by comparison with DSMC results and find excellent match for all macroscopic quantities for ${\rm Kn} \lesssim 0.1$. At ${\rm Kn} \gtrsim 0.1$, small deviations can be seen in the profiles of the diagonal components of the pressure tensor, the heat flux parallel to the plates, and the velocity profile, as well as in the values of the velocity gradient at the channel center. We attribute these deviations to the limited applicability of the Shakhov collision model for highly out of equilibrium flows.

Tables

Table 1. Table 1: Simulation parameters for the results shown in Fig.  3 (a).

$k$	$N_x$	$Q_x$	$S$
$0.03$	$3$	$4$	$64$
$0.1$	$6$	$7$	$32$
$1$	$6$	$21$	$32$
$10$	$6$	$80$	$32$

Table 2. Table 2: Reference values for the length and time for various values of Kn Kn {\rm Kn} , in the context of the simulations discussed in Sec.  V .

$\mathrm{Kn}$	${\tilde{l}}_{\mathrm{ref}}$	${\tilde{t}}_{\mathrm{ref}}$
$0.01$	$0.8833m$	$3.71\mathrm{ms}$
$0.05$	$0.17666m$	$0.741\mathrm{ms}$
$0.1$	$0.08833m$	$0.371\mathrm{ms}$
$0.25$	$0.035332m$	$0.148\mathrm{ms}$
$0.5$	$0.017666m$	$74.1\mu s$
$1.0$	$0.008833m$	$37.1\mu s$

Table 3. Table 3: Simulation parameters for the results shown in Fig.  8 and in Tables  4 , 5 and 6 . The simulations were conducted on a grid with S = 16 𝑆 16 S=16 nodes, stretched according to Eq. ( 46 ) with A = 0.98 𝐴 0.98 A=0.98 . The expansion order of the equilibrium distribution with respect to the half-range Hermite polynomials was always kept at N x = 6 subscript 𝑁 𝑥 6 N_{x}=6 .

$\delta$	$Q_x$	$\delta t$
$\ge 1$	$7$	$5\times {10}^{-4}$
$0.1$	$30$	$2.5\times {10}^{-4}$
$0.01$	$100$	$2.5\times {10}^{-4}$

Table 4. Table 4: Comparison of the results obtained for the velocity gradient ν 𝜈 \nu ( 41 ) using the LB model employed in this paper and the DSMC data reported in Ref.  Sharipov and Strapasson ( 2013 ) , for u w = 2 subscript 𝑢 𝑤 2 u_{w}=\sqrt{2} and various values of δ 𝛿 \delta .

	$\nu$ (Argon)			$\nu$ (Helium)
$\delta$	LB	Ref. Sharipov and Strapasson (2013)	$\epsilon \left({\nu }_{\mathrm{Ar}}\right)$	LB	Ref. Sharipov and Strapasson (2013)	$\epsilon \left({\nu }_{\mathrm{He}}\right)$
$0.01$	$0.027$	$0.048$	$43.8\%$	$0.028$	$0.041$	$31.7\%$
$0.1$	$0.145$	$0.173$	$16.2\%$	$0.149$	$0.172$	$13.40\%$
$1$	$0.484$	$0.486$	$0.42\%$	$0.493$	$0.494$	$0.20\%$
$10$	$0.824$	$0.819$	$0.62\%$	$0.831$	$0.826$	$0.61\%$
$20$	$0.874$	$0.873$	$0.12\%$	$0.880$	$0.880$	$0.0\%$
$40$	$0.905$	$0.904$	$0.12\%$	$0.911$	$0.914$	$0.33\%$

Table 5. Table 5: Comparison of the results obtained for the shear stress Π Π \Pi ( 42 ) using the LB model employed in this paper and the DSMC data reported in Ref.  Sharipov and Strapasson ( 2013 ) , for u w ∈ { 2 , 0.1 ​ 2 } subscript 𝑢 𝑤 2 0.1 2 u_{w}\in\{\sqrt{2},0.1\sqrt{2}\} and various values of δ 𝛿 \delta .

	$\Pi$ (Argon)			$\Pi$ (Helium)
$\delta$	LB	Ref. Sharipov and Strapasson (2013)	$\epsilon \left({\Pi }_{\mathrm{Ar}}\right)$	LB	Ref. Sharipov and Strapasson (2013)	$\epsilon \left({\Pi }_{\mathrm{He}}\right)$
$u_w=\sqrt{2}$
$0.01$	$0.5619$	$0.5612$	$0.12\%$	$0.5618$	$0.5615$	$0.05\%$
$0.1$	$0.5358$	$0.5319$	$0.73\%$	$0.5344$	$0.5315$	$0.55\%$
$1$	$0.3720$	$0.3663$	$1.56\%$	$0.3671$	$0.3616$	$1.52\%$
$10$	$0.09787$	$0.09777$	$0.10\%$	$0.09585$	$0.09551$	$0.36\%$
$20$	$0.05328$	$0.05316$	$0.23\%$	$0.05223$	$0.05191$	$0.62\%$
$40$	$0.02769$	$0.02766$	$0.11\%$	$0.02718$	$0.02704$	$0.52\%$
$u_w=0.1\sqrt{2}$
$0.01$	$0.5594$	$0.5575$	$0.34\%$	$0.5594$	$0.5585$	$0.16\%$
$0.1$	$0.5225$	$0.5167$	$1.12\%$	$0.5225$	$0.5191$	$0.65\%$
$1$	$0.3392$	$0.3365$	$0.81\%$	$0.3392$	$0.3382$	$0.30\%$
$10$	$0.08324$	$0.08320$	$0.05\%$	$0.08322$	$0.08324$	$0.03\%$
$20$	$0.04546$	$0.04531$	$0.34\%$	$0.04545$	$0.04540$	$0.12\%$
$40$	$0.02383$	$0.02381$	$0.09\%$	$0.02382$	$0.02383$	$0.04\%$

Table 6. Table 6: Comparison of the results obtained for the temperature T 0 subscript 𝑇 0 T_{0} at the center of the channel using the LB model employed in this paper and the DSMC data reported in Ref.  Sharipov and Strapasson ( 2013 ) for u w = 2 subscript 𝑢 𝑤 2 u_{w}=\sqrt{2} and various values of δ 𝛿 \delta .

	$T_0$ (Argon)			$T_0$ (Helium)
$\delta$	LB	Ref. Sharipov and Strapasson (2013)	$\epsilon \left(T_0^{\mathrm{Ar}}\right)$	LB	Ref. Sharipov and Strapasson (2013)	$\epsilon \left(T_0^{\mathrm{He}}\right)$
$0.01$	$1.663$	$1.667$	$0.24\%$	$1.663$	$1.665$	$0.12\%$
$0.1$	$1.646$	$1.661$	$0.91\%$	$1.645$	$1.660$	$0.91\%$
$1$	$1.561$	$1.587$	$1.64\%$	$1.558$	$1.583$	$1.58\%$
$10$	$1.357$	$1.360$	$0.22\%$	$1.355$	$1.356$	$0.08\%$
$20$	$1.315$	$1.316$	$0.08\%$	$1.314$	$1.313$	$0.08\%$
$40$	$1.291$	$1.291$	$0\%$	$1.291$	$1.289$	$0.16\%$

Table 7. Table 7: Butcher tableau for the third-order Runge-Kutta integration summarized in Eq. ( 45 ).

0
1	1
1/2	1/4	1/4
	1/6	1/6	2/3

Table 8. Table 8: The limiting values of ω ¯ q subscript ¯ 𝜔 𝑞 \overline{\omega}_{q} ( 51 ) when any combination of indicator of smoothness functions σ i subscript 𝜎 𝑖 \sigma_{i} have vanishing values.

	${\overline{\omega }}_1$	${\overline{\omega }}_2$	${\overline{\omega }}_3$
${\sigma }_1={\sigma }_2={\sigma }_3=0$	$0.1$	$0.6$	$0.3$
${\sigma }_2={\sigma }_3=0$	$0$	$2/3$	$1/3$
${\sigma }_3={\sigma }_1=0$	$1/4$	$0$	$3/4$
${\sigma }_1={\sigma }_2=0$	$1/7$	$6/7$	$0$
${\sigma }_1=0$	$1$	$0$	$0$
${\sigma }_2=0$	$0$	$1$	$0$
${\sigma }_3=0$	$0$	$0$	$1$

Equations

$${\partial }_{\tilde{t}}\tilde{f}+\frac{\tilde{\mathbf{p}}}{\tilde{m}}\cdot \tilde{\nabla }\tilde{f}=-\frac{1}{\tilde{\tau }}\left[\tilde{f}-{\tilde{f}}^{(\mathrm{eq})}\left(1+\mathbb{S}\right)\right],$$

$${\tilde{f}}^{(\mathrm{eq})}=\frac{\tilde{n}}{(2\pi \tilde{m}{\tilde{k}}_B\tilde{T}{)}^{3/2}}\exp \left[-\frac{(\tilde{\mathbf{p}}-\tilde{m}\tilde{\mathbf{u}}{)}^2}{2\tilde{m}{\tilde{k}}_B\tilde{T}}\right],$$

$$\mathbb{S}=\frac{1-\mathrm{Pr}}{\tilde{n}{\tilde{K}}_B^2{\tilde{T}}^2}\left(\frac{{\tilde{\mathbf{\xi }}}^2}{5\tilde{m}{\tilde{K}}_B\tilde{T}}-1\right)\tilde{\mathbf{q}}\cdot \tilde{\mathbf{\xi }},$$

$$\tilde{n}=$$

$$\tilde{\mathbf{u}}=$$

$${\tilde{T}}_{ij}=$$

$$\tilde{\mathbf{q}}=$$

$$\tilde{\mu }=\tilde{\tau }\tilde{n}{\tilde{K}}_B\tilde{T},{\tilde{\kappa }}_T=\frac{1}{\mathrm{Pr}}\frac{5{\tilde{K}}_B}{2\tilde{m}}\tilde{\tau }\tilde{n}{\tilde{K}}_B\tilde{T},$$

$$\mathrm{Pr}=\frac{{\tilde{c}}_p\tilde{\mu }}{{\tilde{\kappa }}_T},{\tilde{c}}_p=\frac{5{\tilde{K}}_B}{2\tilde{m}},$$

$${\tilde{c}}_{\mathrm{ref}}=\sqrt{\frac{{\tilde{K}}_B{\tilde{T}}_{\mathrm{ref}}}{{\tilde{m}}_{\mathrm{ref}}}}=\frac{1}{\sqrt{\gamma }}{\tilde{c}}_{s;\mathrm{ref}},$$

$$\mathrm{Ma}=u\sqrt{\frac{3}{5}}\simeq 0.775u.$$

$$f=\frac{\tilde{f}}{{\tilde{n}}_{\mathrm{ref}}}({\tilde{m}}_{\mathrm{ref}}{\tilde{k}}_B{\tilde{T}}_{\mathrm{ref}}{)}^{3/2},$$

$$f^{(\mathrm{eq})}=\frac{n}{(2\pi mT{)}^{3/2}}\exp \left[-\frac{(\mathbf{p}-m\mathbf{u}{)}^2}{2mT}\right],$$

$${\partial }_{t}f+\frac{\mathbf{p}}{m}\cdot \nabla f=-\frac{1}{\tau }\left[f-f^{(\mathrm{eq})}\left(1+\mathbb{S}\right)\right],$$

$${\partial }_{t}f+\frac{p_x}{m}{\partial }_{x}f=-\frac{1}{\tau }\left[f-f^{(\mathrm{eq})}(1+\mathbb{S})\right].$$

$$\varphi ={\int }_{-\infty }^{\infty }d{p}_{z}f,\chi ={\int }_{-\infty }^{\infty }d{p}_{z}f\frac{p_z^2}{m},$$

$${\partial }_t\left(\begin{array}{c}\varphi \\ \chi \end{array}\right)+\frac{p_x}{m}{\partial }_x\left(\begin{array}{c}\varphi \\ \chi \end{array}\right)=-\frac{1}{\tau }\left(\begin{array}{c}\varphi -{\varphi }^{(\mathrm{eq})}(1+{\mathbb{S}}_{\varphi })\\ \chi -{\chi }^{(\mathrm{eq})}(1+{\mathbb{S}}_{\chi })\end{array}\right),$$

$${\varphi }^{(\mathrm{eq})}=n{g}_x{g}_y,g_{\alpha }=\frac{1}{\sqrt{2\pi mT}}\exp \left[-\frac{(p_{\alpha }-m{u}_{\alpha }{)}^2}{2mT}\right].$$

$${\mathbb{S}}_{\varphi }=$$

$$=$$

$${\mathbb{S}}_{\chi }=$$

$$=$$

$$\varphi =$$

$${\varphi }^{(\mathrm{eq})}(1+{\mathbb{S}}_{\varphi })=$$

$$\omega ({\overline{p}}_y)=\frac{1}{\sqrt{2\pi }}{e}^{-{\overline{p}}_y^2/2}.$$

$${\mathcal{F}}_{\ell }=$$

$${\mathcal{F}}_{\ell }^{{\mathbb{S}}_{\varphi }}=$$

$${\int }_{-\infty }^{\infty }d{\overline{p}}_y\omega ({\overline{p}}_y)H_{\ell }({\overline{p}}_y)H_{{\ell }^{\prime }}({\overline{p}}_y)=\ell !{\delta }_{\ell {\ell }^{\prime }}.$$

$${\partial }_t{\mathcal{F}}_{\ell }+\frac{p_x}{m}{\partial }_x{\mathcal{F}}_{\ell }=-\frac{1}{\tau }\left({\mathcal{F}}_{\ell }-{\mathcal{F}}_{\ell }^{{\mathbb{S}}_{\varphi }}\right),$$

$${\partial }_t\varphi +\frac{p_x}{m}{\partial }_x\varphi =-\frac{1}{\tau }[\varphi -{\varphi }^{(\mathrm{eq})}],$$

$$\tau =\frac{k}{\sqrt{2}},$$

$$u_{\mathrm{slip}}=u_w-u_y(L/2),$$

$$\begin{gathered} u_y(L/2)=\frac{(L/2-x_{S-1})(L/2-x_{S-2})u_{y;S}}{(x_S-x_{S-1})(x_S-x_{S-2})} \\ +\frac{(L/2-x_S)(L/2-x_{S-2})u_{y;S-1}}{(x_{S-1}-x_S)(x_{S-1}-x_{S-2})} \\ +\frac{(L/2-x_S)(L/2-x_{S-1})u_{y;S-2}}{(x_{S-2}-x_S)(x_{S-2}-x_{S-1})}. \end{gathered}$$

$$\begin{gathered} u_y^{\prime }(0)=\frac{x_2^2{x}_3^2{u}_{y;1}}{x_1(x_1^2-x_2^2)(x_1^2-x_3^2)} \\ +\frac{x_1^2{x}_3^2{u}_{y;2}}{x_2(x_2^2-x_1^2)(x_2^2-x_3^2)}+\frac{x_1^2{x}_2^2{u}_{y;3}}{x_3(x_3^2-x_1^2)(x_3^2-x_2^2)}, \end{gathered}$$

Full text

I Introduction

It is generally recognized that the Navier-Stokes-Fourier equations are not appropriate to investigate the flow phenomena in highly rarefied gases where the continuum hypothesis is no longer valid. To investigate such far from equilibrium fluids, the Boltzman equation, which governs the evolution of the one-particle distribution function $f\equiv f(\bm{x},\bm{p},t)$ in a seven-dimensional space, can be employed instead Grad (1958); Kogan (1969); Harris (1971); Gatignol (1975); Cercignani (1988); Bird (1994); Sone (2002); Liboff (2003); Karniadakis et al. (2005); Shen (2005); Sone (2007); Tabeling (2011); Sharipov (2016). Finding solutions of the Boltzmann equation is a challenging task due to the complexity of the collision term, which requires the evaluation of $5$-dimensional integrals over the momentum space. Effective approaches to solve the Boltzmann equation include the celebrated direct simulation Monte Carlo (DSMC) method Bird (1994); Sharipov (2016), where the collision integral is sampled by considering a sufficiently large ensemble of representative particles which are evolved individually; the discrete velocity method (DVM) Broadwell (1964); Sone (2002); Aoki et al. (2003); Sone (2007); Sharipov (2016); and more recently, the fast spectral method, which relies on the projection of the collision term on orthogonal functions Mouhot and Pareschi (2006); Filbet (2012); Wu et al. (2013); Gamba et al. (2017). In all approaches mentioned above, the evaluation of the collision term still remains the most time-consuming part of the numerical algorithm, placing severe constraints on the size and complexity of the systems which can be analyzed numerically.

In the early ’50s, Bhatnaghar, Gross and Krook introduced their single relaxation time approximation of the collision term of the Boltzmann equation describing the ideal gas. This approximation, known as the BGK model, was derived under the assumption that the deviation of the gas from the local (Maxwellian) equilibrium is small Bhatnagar et al. (1954). The severe drawback of the BGK model is that the transport coefficients are governed by the single relaxation time $\tau$, and in particular, the Prandtl number ${\rm Pr}$ is fixed at $1$. This limitation was overcome through the collision term model proposed by Shakhov Shakhov (1968a, b); Yang et al. (1995); Li et al. (2005), who extended the single relaxation time (BGK) model to allow ${\rm Pr}$ to be controlled independently from the relaxation time $\tau$. Other extensions of the BGK model, which allow ${\rm Pr}$ to be controlled, include the ellipsoidal BGK (E-BGK) model Holway, Jr. (1966); Meng et al. (2013a); Zhang et al. (2018) and the multi-relaxation time (MRT) models widely employed in lattice Boltzmann simulations d’Humières (1992); d’Humières et al. (2002). The Shakhov model (also known as the S-model) was subsequently extended by Rykov et al. to account for rotational degrees of freedom Rykov (1975). In the linearized regime, the Gross-Jackson Gross and Jackson (1959) and the McCormack McCormack (1973) models extended the relaxation time paradigm to account for realistic interaction cross-sections for single gases and for gaseous mixtures, respectively.

The simplicity of the relaxation time formulation motivated researchers to implement and develop such models and a surprising range of effects turned out to be correctly recovered. It is now generally accepted that the relaxation time approach can be used to simulate flows which are not far from equilibrium, provided the transport coefficients are correctly recovered Ambru\cbs and Sofonea (2012); Sharipov (2016); Titarev (2018). Solutions of these so-called Boltzmann model equations can be obtained using a variety of numerical methods, amongst which we recall the DVM Ho and Graur (2015); Sharipov (2016), the discrete unified gas-kinetic scheme (DUGKS) Guo et al. (2013, 2015); Peng et al. (2018); Zhu et al. (2017, 2016), the lattice Boltzmann (LB) models Aidun and Clausen (2010); Chen and Doolen (1998); Sukop and Thorne (2006); Shan et al. (2006); Succi (2001); Wolf-Gladrow (2000); Krüger et al. (2017); Succi (2018), the off-lattice Boltzmann models (OLBM) Nannelli and Succi (1992); Lee and Lin (2003); Cheng and Hung (2004); Min and Lee (2011); Fakhari and Lee (2015); Krämer et al. (2017), and the discrete Boltzmann models (DBM) He et al. (1998); Li and Kwok (2003); Xu et al. (2012, 2016); Lin et al. (2017); Xu et al. (2018).

In this paper, we reasses the capabilities of the Shakhov collision model in the context of the three-dimensional ($3D$) Couette flow between parallel plates by employing lattice Boltzmann (LB) models based on the Gauss quadrature method. From a historical perspective, the LB models emerged as successors of the lattice gas automata Succi (2001); Wolf-Gladrow (2000) introduced more than three decades ago and were originally designed to recover the Navier-Stokes equations using a fast and simple algorithm. The efficiency of this algorithm relies on the implementation of the advection and time stepping according to the collide-and-stream paradigm involving the discretization of the momentum space using a relatively small number of momentum vectors which exactly connect neighboring lattice sites Aidun and Clausen (2010); Chen and Doolen (1998); He and Luo (1997a, b); Succi (2001). Besides their numerical efficiency, the key to the success of these early LB models is the implementation of the collision term using a polynomial series which allows its moments to be exactly recovered up to a certain order $N$ Shan et al. (2006). While astonishingly successful at the Navier-Stokes level Aidun and Clausen (2010); Chen and Doolen (1998); Sukop and Thorne (2006); Shan et al. (2006); Succi (2001); Wolf-Gladrow (2000); Krüger et al. (2017); Succi (2018), the LB models based on the collide-and-stream paradigm obtained only limited success when applied to flows of rarefied gases Shi et al. (2011); Meng and Zhang (2011a); Shi et al. (2015). In particular, it was noted that in order to achieve accurate results at non-negligible values of the Knudsen number ${\rm Kn}$, the velocity set must be enriched to account for higher order moments of the distribution function. This brings about a series of complications, since during a single time step, the particles must hop over an increasing number of lattice sites Chikatamarla and Karlin (2009); Namburi et al. (2016); Philippi et al. (2006); Sbragaglia et al. (2007). This also renders the implementation of boundary conditions cumbersome Meng and Zhang (2011a, 2014).

A straightforward alternative to the collide-and-stream paradigm comes from the DVM implementations of the Boltzmann equation, where the advection and time stepping are implemented using finite differences. Retaining the simplified polynomial truncation of the collision term, these implementations can be referred to as finite difference lattice Boltzmann (FDLB) models Cao et al. (1997); Chen (1998); He (1997); Mei and Shyy (1998); Ubertini et al. (2003); Watari and Tsutahara (2004); Bardow et al. (2008); Gan et al. (2011); Meng and Zhang (2011b); Meng et al. (2011), or discrete Boltzmann models He et al. (1998); Li and Kwok (2003); Xu et al. (2012, 2016); Lin et al. (2017); Xu et al. (2018). Amongst the first high-order FDLB models are the 2D shell-based models introduced by Watari and Tsutahara Watari and Tsutahara (2003); Watari (2007, 2010, 2016) and their 3D generalizations Watari and Tsutahara (2006); Watari (2009); Romatschke et al. (2011); Ambru\cbs and Sofonea (2012), the models based on the tensor Hermite polynomials Shan et al. (2006); Meng and Zhang (2011b); Meng et al. (2011), and, more recently, the models based on the Cartesian split of the velocity space discussed in Refs. Ambru\cbs and Sofonea (2014a, b, c, 2016a, 2016b).

The Couette flow between parallel plates has become a benchmark problem for fluid dynamics simulations. In the context of rarefied gases, the linearized Boltzmann-BGK equation has been solved in a semi-analytic manner to high numerical precision in Ref. Jiang and Luo (2016) (see also Refs. Sone (2007); Yap and Sader (2012); Li et al. (2015) for previous results). The LB results obtained within the collide-and-stream paradigm reported excellent agreement with DVM or DSMC results for small values of ${\rm Kn}$, but their performance quickly deteriorated as ${\rm Kn}$ was increased towards the transition regime Ansumali et al. (2007); Yudistiawan et al. (2008, 2010); Shi et al. (2011); Meng and Zhang (2011a); Shi et al. (2015). There has been significant effort devoted to deriving macroscopic equations which account for the non-equilibrium features appearing at non-negligible ${\rm Kn}$, from which we only mention the regularized 13 moments system of equations (R13) Struchtrup and Torrilhon (2007); Torrilhon and Struchtrup (2008); Taheri et al. (2009); Struchtrup and Taheri (2011); Rana and Struchtrup (2016). The limited success of the R13 system in the transition regime suggests that even more moments should be taken into account. From a lattice Boltzmann perspective, this is equivalent to extending the quadrature order of the model. Indeed, this was confirmed in Ref. Ambru\cbs and Sofonea (2012), when it was shown that FDLB models based on the spherical factorization of the momentum space exhibit a slow but steady convergence with respect to the quadrature order at ${\rm Kn}=0.5$.

Recent FDLB studies of the Couette flow showed that accurate solutions of the S-model equation can be obtained by enriching the velocity set Ambru\cbs and Sofonea (2012, 2014a, 2014c, 2016a). The crucial piece which makes the simulations much more efficient is to take into account the discontinuity in the distribution function, which is induced by the boundary conditions prescribed at the solid walls Gross et al. (1957); Takata and Funagane (2013). For this reason, appropriate half-range quadratures should be employed on the direction perpendicular to the walls Yang et al. (1995); Li and Zhang (2003, 2004, 2009); Lorenzani et al. (2007); Frezzotti et al. (2009, 2011); Gibelli (2012); Guo et al. (2013); Ghiroldi and Gibelli (2014); Ambru\cbs and Sofonea (2014c, b, a); Guo et al. (2015); Ghiroldi and Gibelli (2015); Shi et al. (2015); Ambru\cbs and Sofonea (2016a, b). In this paper, we will employ the mixed quadrature LB models introduced in Ref. Ambru\cbs and Sofonea (2016a) for the study of the 2D Couette flow.

Since we are interested in obtaining accurate results for the temperature field, as well as for the heat flux, the number of degrees of freedom in the momentum space must be equal to three. In the context of flows which are effectively two-dimensional (e.g., in the $xy$ plane), this can be achieved by working with reduced distributions Li and Zhang (2004); Meng et al. (2013b). These distributions are obtained by analytically integrating the momentum space degree of freedom along the direction which is perpendicular to the normal to the walls and to the direction of the flow. The ensuing LB models employ the half-range Gauss-Hermite quadrature on the axis perpendicular to the walls (the $x$ axis) and the full-range Gauss-Hermite quadrature for the axis parallel to the wall (the $y$ axis). The advantage of these LB models with mixed quadratures is described as follows. The discontinuity induced by the diffuse reflection boundary conditions imposed on the walls perpendicular to the $x$ axis warrants the use of the half-range Gauss-Hermite quadrature, which requires $Q=N+1$ points on each Cartesian semiaxis to ensure the accurate recovery of the moments of the distribution function $f$ up to $N$th order. The $2Q$ quadrature points on the whole axis are twice the number of points required by the full-range Gauss-Hermite quadrature to achieve the same degree of accuracy when considering the full-space moments of $f$ Ambru\cbs and Sofonea (2016a). Thus, it is more convenient to employ the full-range Gauss-Hermite quadratures on the $y$ direction (i.e., the direction parallel to the walls), where no discontinuities in the distribution function arise, resulting in an overall smaller velocity set.

The aim of this paper is to demonstrate that the Shakhov collision term model can be successfully employed in LB simulations to match Direct Simulation Monte Carlo (DSMC) results. In particular, we consider comparisons with the data obtained for the Couette flow of Maxwell molecules Struchtrup and Torrilhon (2007, 2008); Torrilhon and Struchtrup (2008); Taheri et al. (2009); Schuetze , as well as for the flow of Helium and Argon modeled using ab initio potentials Sharipov and Strapasson (2013).

For the spatial advection, we employed the fifth order weighted essentially non-oscillatory (WENO-5) scheme described in Refs. Jiang and Shu (1996); Gan et al. (2011); Ambru\cbs and Blaga (2018); Hejranfar et al. (2017); Busuioc and Ambru\cbs , while the time stepping was performed using a third-order Runge-Kutta algorithm Shu and Osher (1988); Gottlieb and Shu (1998); Henrick et al. (2005); Trangenstein (2007); Ambru\cbs and Blaga (2018); Busuioc and Ambru\cbs . The implementation of the diffuse reflection boundary conditions on the channel walls is identical to that presented in Ref. Busuioc and Ambru\cbs .

The outline of this paper is as follows. In Sec. II, the Shakhov collision term is briefly described. Section III is dedicated to discussing the Boltzmann equation for the $3D$ Couette flow between parallel plates, in the context of reduced distributions. Section IV validates our models in the incompressible (low Mach) regime, by comparison with the semi-analytic benchmark data reported in Ref. Jiang and Luo (2016). In Sec. V, our models are validated against the DSMC results for Maxwell molecules reported in Refs. Struchtrup and Torrilhon (2007, 2008); Torrilhon and Struchtrup (2008); Taheri et al. (2009); Schuetze , at various values of ${\rm Kn}$ and of the wall velocity $u_{w}$. In Sec. VI, we consider a comparison with the DSMC simulation results based on ab-initio potentials for Helium and Argon Sharipov and Strapasson (2013). Section VII concludes this paper. The details regarding our numerical scheme are presented in Appendix A. We warn our readers that, in order to facilitate the comparison of our LB results with various data in the literature, the notation used to refer to the degree of rarefaction varies between the sections where these results are reported, as follows: in Sec. IV, the notation $k$ is employed to refer to the Knudsen number; in Sec. V, the Knudsen number is denoted using the familiar notation ${\rm Kn}$ and is linked to $k$ through ${\rm Kn}=k/\sqrt{2}$; in Sec. VI, the degree of rarefaction is characterized using the rarefaction parameter $\delta$, which is linked to $k$ (from Sec. IV) and ${\rm Kn}$ (from Sec. V) through $\delta=1/k$ and $\delta=1/{\rm Kn}\sqrt{2}$, respectively.

II Boltzmann equation with the Shakhov collision term

In this section, we establish our notation by introducing the Shakhov collision term model, as well asour non-dimensionalization convention.

II.1 Shakhov model

The Boltzmann equation for a force-free flow with the Shakhov collision term is given by Shakhov (1968a, b); Guo et al. (2015); Titarev (2007); Graur and Polikarpov (2009); Ambru\cbs and Sofonea (2012):

[TABLE]

where the overhead tilde $\widetilde{\phantom{f}}$ denotes dimensionful quantities. In the above, $\widetilde{f}$ is the Boltzmann distribution function, $\widetilde{m}$ and $\widetilde{{\bm{p}}}$ are the particle mass and momentum vector, respectively, while $\widetilde{\tau}$ is the relaxation time. The Maxwell-Boltzmann equilibrium distribution function is given by:

[TABLE]

where $\widetilde{n}$, $\widetilde{T}$ and $\widetilde{u}$ are the macroscopic density, temperature and velocity of the fluid. The Shakhov term $\mathbb{S}$ is given by:

[TABLE]

where ${\rm Pr}$ gives the Prandtl number (see below), $\widetilde{{\bm{\xi}}}=\widetilde{{\bm{p}}}-\widetilde{m}\widetilde{{\bm{u}}}$ is the peculiar momentum, while $\widetilde{{\bm{q}}}$ is the heat flux.

The macroscopic quantities $\widetilde{n}$, $\widetilde{{\bm{u}}}$, pressure tensor $\widetilde{T}_{ij}$ and $\widetilde{{\bm{q}}}$ can be obtained as moments of $\widetilde{f}$:

[TABLE]

while the pressure is obtained as $\widetilde{P}=\frac{1}{3}(\widetilde{T}_{xx}+\widetilde{T}_{yy}+\widetilde{T}_{zz})=\widetilde{n}\widetilde{K}_{B}\widetilde{T}$.

By employing the Chapman-Enskog expansion, it can be seen that, for three-dimensional flows, the Shakhov collision term gives rise to the following expressions for the transport coefficients, namely the dynamic (shear) viscosity $\widetilde{\mu}$ and the heat conductivity $\widetilde{\kappa}_{T}$ Ambru\cbs and Sofonea (2012):

[TABLE]

where the Prandtl number ${\rm Pr}$, calculated as

[TABLE]

is adjustable according to Eq. (3).

II.2 Non-dimensionalization convention

Our non-dimensionalization convention follows the one employed in Ref. Sofonea and Sekerka (2005), being based on the following reference quantities: the reference length $\widetilde{l}_{\rm ref}=\widetilde{L}$, the reference temperature $\widetilde{T}_{\rm ref}$, the reference mass $\widetilde{m}_{\rm ref}$, and the reference density $\widetilde{n}_{\rm ref}$. In the context of the Couette flow, the reference length is taken as the distance between the parallel plates, $\widetilde{l}_{\rm ref}=\widetilde{L}$, and the reference temperature is the wall temperature, $\widetilde{T}_{\rm ref}=\widetilde{T}_{w}$. The reference mass is taken equal to the particle mass $\widetilde{m}_{\rm ref}$, while the reference speed is:

[TABLE]

where $\widetilde{c}_{s;{\rm ref}}$ is the sound speed at the reference temperature and $\gamma$ is the adiabatic index. Since we only consider monatomic ideal gases ($\gamma=5/3$), the Mach number corresponding to the non-dimensionalized velocity $u$ is:

[TABLE]

The non-dimensionalized distribution function $f$ is defined as:

[TABLE]

such that the non-dimensional form of the Maxwell-Boltzmann distribution (2) is:

[TABLE]

where ${\bm{p}}\equiv\widetilde{{\bm{p}}}/\widetilde{m}_{\rm ref}\widetilde{c}_{\rm ref}$ and the non-dimensional mass $m=\widetilde{m}/\widetilde{m}_{\rm ref}=1$ is kept explicitly for the sake of clarity of the mathematical relations presented in what follows. Finally, the Boltzmann equation (1) can be non-dimensionalized after multiplying both sides by $\widetilde{t}_{\rm ref}(\widetilde{m}_{\rm ref}\widetilde{k}_{B}\widetilde{T}_{\rm ref})^{3/2}/\widetilde{n}_{\rm ref}$, yielding:

[TABLE]

The exact expression for the nondimensionalized relaxation time $\tau$ depends on the collision term model and will be discussed separately in Secs. IV, V and VI.

III Lattice Boltzmann models for the Couette flow

In the remainder of this paper, we focus on the Couette flow between parallel plates. The coordinate system is chosen such that the $x$ axis is perpendicular to the plates, which are located at $x_{\rm left}=-L/2$ and $x_{\rm right}=L/2$, where the non-dimensionalized channel length $L=1$ is kept explicitly in order to facilitate the physical interpretation of the mathematical expressions appearing below. The plates are set in motion along the $y$ direction with constant velocites $u_{\rm left}=-u_{w}$ and $u_{\rm right}=u_{w}$, while their temperature is kept constant ($T_{\rm left}=T_{\rm right}=T_{w}$).

In order to simulate the Couette flow, we employ the LB models presented in Ref. Ambru\cbs and Sofonea (2016a). While these models were used in Ref. Ambru\cbs and Sofonea (2016a) for the simulation of the $2D$ Couette flow in the Boltzmann-BGK model, their applicability to the present problem is immediate since the extra momentum space degree of freedom $p_{z}$, which is perpendicular to the flow direction and to the normal to the walls, can be eliminated by introducing the reduced distributions $\phi$ and $\chi$, as discussed below. It was shown in Ref. Ambru\cbs and Sofonea (2016a) that, for the simulation of flows at non-negligible values of ${\rm Kn}$, the half-range Gauss-Hermite quadrature should be used on the axis perpendicular to the walls (the $x$ axis), while a relatively low-order full-range Gauss-Hermite quadrature is adequate for the direction parallel to the walls. The resulting models are denoted by ${\rm HHLB}(N_{x};Q_{x})\times{\rm HLB}(N_{y};Q_{y})$, where $N_{x}$ and $N_{y}$ are the orders of the polynomial expansions of the equilibrium distribution along the $x$ and $y$ axes, respectively. $Q_{x}$ and $Q_{y}$ denote the quadrature orders employed on the $x$ semiaxes and on the full $y$ axis of the momentum space. The resulting velocity set comprises $2Q_{x}\times Q_{y}$ vectors.

In Ref. Ambru\cbs and Sofonea (2016a) it was shown, for a two-dimensional mixed quadrature LB model, that the value of $Q_{x}$ required to obtain accurate simulation results depends on ${\rm Kn}$ and on the wall velocity $u_{w}$. Furthermore, it was also shown that, in the particular case of the BGK implementation of the collision term, the simulations employing $Q_{y}\geq 4$ at fixed $Q_{x}$ produced identical results. In this section, a similar inequality will be derived in the case of the Shakhov model.

Since the Couette flow is completely homogeneous along the directions which are parallel to the walls (i.e. the $y$ and $z$ directions), Eq. (11) reduces to:

[TABLE]

The $p_{z}$ degree of freedom of the momentum space can be eliminated by integrating Eq. (12) with respect to $p_{z}$. Defining the reduced distributions $\phi\equiv\phi(x,p_{x},p_{y},t)$ and $\chi\equiv\chi(x,p_{x},p_{y},t)$ via:

[TABLE]

the following two equations are obtained:

[TABLE]

where $\chi^{\rm(eq)}=T\phi^{\rm(eq)}$ and $\phi^{\rm(eq)}$ is just the 2D Maxwell-Boltzmann distribution ($\alpha\in\{x,y\}$):

[TABLE]

The terms $\mathbb{S}_{\phi}$ and $\mathbb{S}_{\chi}$ are given by:

[TABLE]

Let us now expand $\phi$ and $\phi^{\rm(eq)}(1+\mathbb{S}_{\phi})$ with respect to the full-range Hermite polynomials on the $y$ axis, as follows Hildebrand (1987):

[TABLE]

where ${\overline{p}}_{y}\equiv p_{y}/p_{0,y}$ represents the particle momentum along the $y$ axis in units of an arbitrary momentum scale $p_{0,y}$ (we only consider $p_{0,y}=1$ in this paper), while the weight function $\omega({\overline{p}}_{y})$ for the full-range Hermite polynomials is

[TABLE]

Furthermore, the coefficients $\mathcal{F}_{\ell}$ and $\mathcal{F}^{\mathbb{S}_{\phi}}_{\ell}$ can be computed as:

[TABLE]

The compatibility between Eqs. (17) and (19) is ensured by the orthogonality relation obeyed by the Hermite polynomials Olver et al. (2010):

[TABLE]

Substituting (17) into the Boltzmann equation (12) and projecting on the space of Hermite polynomials gives:

[TABLE]

where, as before, $\mathcal{F}_{\ell}$ and $\mathcal{F}^{\mathbb{S}_{\phi}}_{\ell}$ depend on $x$, $p_{x}$ and $t$. Since $\phi^{\rm(eq)}$ depends on $n$, ${\bm{u}}$ and $T$, while $\mathbb{S}_{\phi}$ depends on ${\bm{q}}$, the coefficients $\mathcal{F}^{\mathbb{S}_{\phi}}_{\ell}$ always involve the coefficients $\mathcal{F}_{\ell^{\prime}}$ of orders $0\leq\ell^{\prime}\leq 3$. In this paper, only the moments of $f$ up to $\bm{q}$ are tracked. These moments can be expressed with respect to the coefficients $\mathcal{F}_{\ell}$ with $0\leq\ell\leq 3$. In our implementation, we consider the expansion (15) of $\phi^{\rm(eq)}$ through a direct product procedure, by separately expanding the factors $g_{x}$ and $g_{y}$. More precisely, $g_{x}$ is expanded with respect to the half-range Hermite polynomials and $g_{y}$ is expanded with respect to the full-range Hermite polynomials up to order $N_{y}=Q_{y}-1$ Ambru\cbs and Sofonea (2016a, b). Since $\mathbb{S}_{\phi}$ is a polynomial of the third order in ${\bm{p}}$, the evolution of ${\bm{q}}$ requires the recovery of the moments of $\phi^{\rm(eq)}$ of order $3+3=6$ on each axis. Such moments can be exactly recovered when $N_{y}\geq 6$.

The above analysis shows that, when the Shakhov collision term is employed, the numerical results for the $3D$ Couette flow considered in this paper when $Q_{y}>7$ must coincide with those obtained with $Q_{y}=7$. Furthermore, as discussed in Ref. Ambru\cbs and Sofonea (2016a), the accurate simulation of flows at large values of ${\rm Kn}$ requires the increase of the quadrature $Q_{x}$. It is worth mentioning that if ${\rm Pr}=1$, the Shakhov term $\mathbb{S}$ in Eq. (3) vanishes and the LB models with $Q_{y}\geq 4$ give identical results, as discussed in Ref. Ambru\cbs and Sofonea (2016a).

IV Low Mach number validation

For low Mach number flows, the Boltzmann equation can be approximated via its linearized form, which was tackled by many authors using numerical simulations Cercignani (1969); Kim et al. (2008); Sharipov (2009, 2016); Wu et al. (2017) or semi-analytic methods Sone (2007); Yap and Sader (2012); Li et al. (2015); Jiang and Luo (2016). In order to approach the assumptions of the linearized regime, the simulations presented in this section are performed with $u_{w}=10^{-5}$, such that the temperature and density profiles remain nearly constant throughout the channel. At higher values of $u_{w}$, deviations from the benchmark results can be expected, since the profile of $u_{y}/u_{w}$ does not scale perfectly with $u_{w}$, as also demonstrated in Ref. Rogozin (2016). At $u_{w}=10^{-5}$, a quadrature order $Q_{y}=2$ is sufficient for the direction pointing along the flow.

Since the value of ${\rm Pr}$ does not influence the results in this regime, the simulations are performed with ${\rm Pr}=1$ (the BGK approximation). Furthermore, since the linearized regime analysis is oblivious of the number of degrees of freedom of the particle constituents, our simulations are performed only at the level of the $\phi$ function, as follows:

[TABLE]

where $\phi^{\rm(eq)}$ is given in Eq. (15) with the temperature obtained as $T=\frac{1}{2n}(T_{xx}+T_{yy})$. The relaxation time is implemented through:

[TABLE]

where $k$ is interpreted as the Knudsen number Sone (2007); Li et al. (2015); Jiang and Luo (2016).

The numerical results presented in this section were obtained with the ${\rm HHLB}(N_{x};Q_{x})\times{\rm HLB}(1;2)$ models, for various values of $Q_{x}$. For consistency with the subsequent numerical results sections, the expansion order $N_{x}$ was set to $N_{x}={\rm min}(6,Q_{x}-1)$. The spatial domain was discretized using $S$ nodes, stretched according to Eq. (47) with $A=0.98$. More details on the numerical scheme are provided in the Appendix.

In the first part, we consider a comparison with the benchmark data obtained by Jiang and Luo in Ref. Jiang and Luo (2016) for the following four quantities: the slip velocity $u_{\rm slip}$, the velocity derivative at the channel center $u^{\prime}(0)$, the mass flow rate $\dot{m}$, and the non-diagonal stress $T_{xy}$, which are introduced below.

The slip velocity is obtained by subtracting the fluid velocity at the wall $u_{y}(L/2)$ from the wall velocity:

[TABLE]

where $u_{y}(L/2)$ is obtained by quadratic extrapolation from the fluid nodes with indices $S$, $S-1$ and $S-2$ ($S$ is the index of the last node inside the fluid domain and $u_{y;i}$ is the $y$ component of the fluid velocity in node $i,\,1\leq i\leq S$):

[TABLE]

This expression of $u_{y}(L/2)$ is third order accurate with respect to the spacing $\delta\eta$ of the parameter of the stretched grid, which is defined in the Appendix.

The velocity derivative at the channel center is obtained using the following three-point formula:

[TABLE]

which is sixth order accurate if we take into account the antisymmetry of the velocity profile with respect to the channel center.

The half-channel mass flow rate is obtained using the rectangle integration method:

[TABLE]

which is second order accurate with respect to $\delta\eta$. In the hydrodynamic regime, when $u_{y}^{\rm hydro}=2u_{w}x/L$ Kundu et al. (2015); Rieutord (2015), the half-channel mass flow rate is given by:

[TABLE]

Finally, the non-diagonal stress $T_{xy}$ is obtained by averaging $T_{xy}$ over the half-channel, using the equivalent of Eq. (27):

[TABLE]

In the ballistic regime, $T_{xy}$ is given by Sone (2007); Ambru\cbs and Sofonea (2016a):

[TABLE]

In Fig. 1, the results obtained with our mixed quadrature LB models are compared with those reported by Jiang and Luo in Table 2 of Ref. Jiang and Luo (2016), for (a) the relative slip velocity $u_{\rm slip}/u_{w}$, (b) the normalized velocity derivative at the channel center $Lu^{\prime}_{y}(0)/2u_{w}$, (c) the normalized mass flow rate $\dot{m}/\dot{m}_{\rm hydro}$ and (d) the normalized non-diagonal stress $T_{xy}/T_{xy}^{\rm bal}$. It can be seen that very reasonable agreement is obtained with the ${\rm HHLB}(3;4)\times{\rm HLB}(1;2)$ model on a grid with $S=16$ nodes. The small discrepancies seen in Fig. 1(b) at small values of $k$ are removed by doubling the grid points, while the discrepancies observed at large values of $k$ in Figs. 1(a-c) can be removed by increasing the quadrature order $Q_{x}$ of the half-range Gauss-Hermite quadrature.

The comparison shown in Fig. 1 allows a qualitative assessment to be made at the level of absolute differences between our results and the benchmark results. We now discuss the relative error $\varepsilon(A)$, which is defined for a quantity $A$ as

[TABLE]

where $A_{\rm LB}$ and $A_{\rm lin}$ are the values of $A$ obtained using our models and those reported in Ref. Jiang and Luo (2016), respectively. The use of $\varepsilon(A)$ augments the differences between our LB results and the benchmark data in the regions where $A$ is small. In Fig. 2(a-d), the relative errors $\varepsilon(A)$, computed for $A\in\{u_{y}(L/2),u^{\prime}_{y}(0),\dot{m},T_{xy}\}$, are represented with respect to $k$ (23) for various quadrature orders. The dashed lines in panels (a–c) and (d) indicate the $1\%$ and $0.1\%$ relative error thresholds, respectively. It can be seen that for the quantities (a) $u_{y}(L/2)$, (b) $u^{\prime}_{y}(0)$, and (c) $\dot{m}$, the relative error at high values of $k$ can be decreased below $1\%$ only when high quadrature orders $Q_{x}$ are employed. This is because these three quantities go to [math] as $k\rightarrow\infty$ and hence the absolute error must decrease significantly in order to achieve the $1\%$ threshold for the relative error. By contrast, the relative error in $T_{xy}$, shown in Fig. 2(d), is well below $1\%$ even when $Q_{x}=4$. For this quantity, the relative error becomes less that $0.1\%$ when $Q_{x}\geq 7$. All results presented in Fig. 2 were obtained on a $1D$ grid with $S=16$ points, stretched according to Eq. (46) with $A=0.98$, using the models ${\rm HHLB}(N_{x};Q_{x})\times{\rm HLB}(1;2)$ of various quadrature orders $Q_{x}$ and $N_{x}={\rm min}(6,Q_{x}-1)$.

In the second and final part, a comparison between our results and those reported in Ref. Li et al. (2015) for the nonlinear part of the velocity profile $u_{y}^{\rm nl}$ is considered. This nonlinear part refers to the departure of the solution of the kinetic equation from the straight line profile predicted via the Navier-Stokes equations. The construction of $u_{y}^{\rm nl}$ is made by first obtaining a linear velocity profile which vanishes at the channel center and which attains the value predicted by the kinetic equation on the channel wall. This linear profile can be regarded as the solution of the Navier-Stokes equations with the correct velocity slip taken into account. Subtracting the velocity profile obtained by solving the kinetic equation gives a profile which vanishes, by construction, at the channel center and at the wall. This nonlinear velocity profile $u_{y}^{\rm nl}$ is defined as:

[TABLE]

It can be seen in Fig. 3(a) that the resulting profiles cancel at $x=0$ and $x=L/2$, while reaching a maximum value inside the channel. The value of this maximum increases with $k$ up to a maximum value, and decreases afterwards as $k\rightarrow\infty$, when the velocity profile is trivially a straight line: $u_{y}^{\rm bal}(x)=0$. Our results are in excellent agreement with those reported in Ref. Li et al. (2015) and were obtained using the models ${\rm HHLB}(N_{x};Q_{x})\times{\rm HLB}(1;2)$, with $N_{x}={\rm min}(Q_{x}-1,6)$, while $Q_{x}=4,7,21$ and $80$ for $k=0.03$, $0.1$, $1$ and $10$, respectively. In order to maintain a good accuracy, the number of grid points had to be increased to $S=64$ for $k=0.03$ and $S=32$ for $k\in\{0.1,\,1,\,10\}$. These simulation parameters are summarized in Table 1. For $k=0.03$, the approach of our LB results towards the benchmark solution as the grid is successively refined can be seen in Fig. 3(b).

V Comparison with DSMC: Maxwell molecules

We now benchmark our LB results against the direct simulation Monte Carlo (DSMC) results for Maxwell molecules reported by Struchtrup et al in Refs. Struchtrup and Torrilhon (2007, 2008); Schuetze . The working gas in these simulations is Argon, such that the reference mass is taken to be $\widetilde{m}_{\rm ref}=\widetilde{m}_{\rm Ar}\simeq 6.63\times 10^{-26}\,{\rm kg}$. Taking the reference (and hence, the wall) temperature to be $\widetilde{T}_{\rm ref}=\widetilde{T}_{w}=273\,{\rm K}$, the reference velocity (7) is $\widetilde{c}_{\rm ref}\simeq 238.35\,{\rm m}/{\rm s}$ and the sound speed is $\widetilde{c}_{s}\simeq 307.71\,{\rm m}/{\rm s}$ at $\widetilde{T}=\widetilde{T}_{w}$. The average particle number density was taken to be $\widetilde{n}=1.4\times 10^{20}\,{\rm molecules}/{\rm m}^{3}$ and the mean free path is $\widetilde{\lambda}=0.008833\,{\rm m}$. The Knudsen number ${\rm Kn}$ is thus controlled by varying the domain size, such that the reference length and reference time depend on ${\rm Kn}$, as shown in Table 2.

In the Maxwell molecules model, the viscosity coefficient has a linear temperature dependence:

[TABLE]

where $\widetilde{\mu}_{\rm ref}$ is the viscosity at the reference temperature. This expression for the viscosity can be achieved within the single relaxation time approximation by setting the non-dimensionalized relaxation time to Bird (1994):

[TABLE]

This expression ensures that the viscosity obtained via the Chapman-Enskog expansion is linear with respect to the temperature:

[TABLE]

The Prandtl number is fixed at ${\rm Pr}=2/3$ using the Shakhov model, as discussed in Sec. II.

In this section, we compare the results obtained using our LB models and the DSMC results at the level of the profiles of the density $n$, pressure $P$, temperature $T$, viscous stress $\Pi_{ij}=T_{ij}-nT\delta_{ij}$, velocity $u_{y}$ and heat fluxes $q_{x}$ and $q_{y}$. We considered three batches of simulations, which are discussed below. For all simulations, a grid with $S=16$ nodes, stretched according to Eq. (46) with $A=0.98$, was employed. The time step was set to $\delta t=5\times 10^{-4}$. The model used was ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ for all simulations, except at ${\rm Kn}=1$. Since at ${\rm Kn}=1$, the flow enters the transition regime, the quadrature order had to be increased to $Q_{x}=11$ and the simulations in this regime were performed with the model ${\rm HHLB}(6;11)\times{\rm HLB}(6;7)$. We note that increasing the expansion order of $g_{x}$ with respect to the half-range Hermite polynomials from $N_{x}=6$ to higher values does not have any visible influence on the results.

In the first batch of simulations, the relative wall velocity difference is fixed at $2\widetilde{u}_{w}=300\ {\rm m}/{\rm s}$ ($u_{w}\simeq 0.63$) and ${\rm Kn}$ is varied from $0.01$ up to $1$. Our simulation results are shown in Fig. 4 alongside the DSMC results. An excellent agreement can be seen for the velocity (a), the heat fluxes [(b) and (c)], the density (d) and the non-diagonal component $\Pi_{xy}$ of the viscous stress tensor (i). The temperature, the pressure and the diagonal components of the stress tensor present visible deviations when ${\rm Kn}\gtrsim 0.25$.

The second and third simulation batches are performed at ${\rm Kn}=0.1$ and ${\rm Kn}=0.5$, respectively, for values of the relative wall velocity difference $2\widetilde{u}_{w}$ between $200\ {\rm m}/{\rm s}$ (${\rm Ma}\simeq 0.65$) and $1000\ {\rm m}/{\rm s}$ (${\rm Ma}\simeq 3.25$), corresponding to $u_{w}\simeq 0.42$ and $u_{w}\simeq 2.1$, respectively. The simulation results are shown in Figs. 5 and 6. In these figures, a very good agreement can be seen between the LB and DSMC results. It is interesting to note that there is some disagreement in the results for $\Pi_{ij}$, $T$ and $P$ at high wall velocities ($u_{w}\gtrsim 1.68$), even at ${\rm Kn}=0.1$. This disagreement seems to indicate that the relaxation time model becomes inaccurate at high shearing rates.

The comparisons presented in this section validate the LB models with mixed Gauss-Hermite quadratures, for a wide range of the Knudsen number, as well as of the plate velocities. The simulations were performed using a discretization of the velocity space employing $2Q_{x}Q_{y}=98$ distinct vectors ($Q_{x}=Q_{y}=7$) for ${\rm Kn}<1$ and $154$ distinct velocities ($Q_{x}=11$, $Q_{y}=7$) at ${\rm Kn}=1$. This makes our proposed method highly efficient for the study of channel flows in the rarefied regime.

VI Comparison with DSMC: Realistic potentials

We now consider the validation of our LB models in the case of the Couette flow of two noble gases, namely Helium (He) and Argon (Ar). The transport coefficients for these gases were measured experimentally and the experimental data can be found in Ref. Kestin et al. (1984). More recently, these transport coefficients were computed using ab initio potentials and the results were reported in Refs. Cencek et al. (2012) and Vogel et al. (2010) for dilute Helium and Argon gases. In this subsection, we compare our LB simulation results with the results reported in Ref. Sharipov and Strapasson (2013), calculated using the direct simulation Monte Carlo (DSMC) method based on ab initio potentials over a wide range of the gas rarefaction parameter $\delta$. The results reported in Ref. Sharipov and Strapasson (2013) concern mixtures of He and Ar, including the limiting cases of pure He and pure Ar. The treatment of gas mixtures requires more elaborate models, such as the McCormack model, which was introduced in Ref. McCormack (1973) for the linearized Boltzmann equation. Outside the linear regime, a relaxation time model for isothermal binary fluids was proposed in Ref. Sofonea and Sekerka (2001). To the best of our knowledge, there is no established relaxation time model which can be used to simulate non-isothermal gas mixtures. For simplicity, in this subsection we only consider the case of pure monatomic gases (He or Ar), which can be easily treated in the framework of the Shakhov model.

The Maxwell molecules model can be considered as a particular case of interaction model, for which the dynamic viscosity $\widetilde{\mu}$ varies with temperature according to:

[TABLE]

where $\widetilde{\mu}_{\rm ref}$ is the value of the viscosity at the reference temperature $\widetilde{T}_{\rm ref}$ and the viscosity index $\omega$ takes the value $\omega_{\rm Maxwell}=1$ in the case of Maxwell molecules. In a more general formulation, the variable hard sphere (VHS) model gives rise to values of the viscosity index of the form $\omega=\frac{1}{2}+\nu$, where $\nu$ is a constant controlling the dependence of the collisional cross section on the relative speed of the interacting particles Bird (1994). The particular case of the hard sphere (HS) model is recovered by setting $\nu=0$. The VHS paradigm works remarkably well for quasi-isothermal flows. However, when the temperature variations in the flow are large, the data tabulated in Refs. Cencek et al. (2012); Vogel et al. (2010) for the transport coefficients of various gases indicate that the viscosity index $\omega$ is a slowly-varying function of the temperature Sharipov and Strapasson (2013).

In order to obtain reasonable agreement between the collisional model employed in DSMC and the realistic data available for the transport coefficients, the generalized hard sphere (GHS) model was introduced in Ref. Hassan and Hash (1993). This model can reproduce with remarkable accuracy the experimental data in Ref. Kestin et al. (1984) to temperatures as high as $1500\ {\rm K}$, where the deviation from the experimental data is about $5\%$ Macrossan (2003).

In this section, we consider the viscosity law which emerges from the model proposed by Sutherland Sutherland (1893). In this model, the interaction potential is considered to be the superposition between a hard-sphere-like infinite potential barrier around the repulsive core, followed by an attractive tail. This potential can be used to derive the scattering cross-sections of two-particle collisions and the value of the viscosity coefficient emerging from this model is given by Sutherland (1893):

[TABLE]

where $\widetilde{\mu}_{\rm ref}$ is the fluid viscosity at the temperature $\widetilde{T}=\widetilde{T}_{\rm ref}$ and $\widetilde{S}$ is Sutherland’s constant, which has the dimension of temperature. According to Ref. Macrossan (2003), the Sutherland model can be fitted in order to reproduce with reasonable accuracy the experimental data reported in Ref. Kestin et al. (1984).

In order to validate our simulation results against the results reported in Ref. Sharipov and Strapasson (2013), we take the reference temperature to be equal to the wall temperature, which was set to $\widetilde{T}_{\rm ref}=300\ {\rm K}$ therein. At this temperature, the values of the viscosity for ${\rm He}$ and ${\rm Ar}$ reported in Ref. Kestin et al. (1984) are $\widetilde{\mu}_{\rm ref}^{\rm He}\simeq 20.04\ {\rm\mu Pa}\cdot{\rm s}$ and $\widetilde{\mu}_{\rm ref}^{\rm Ar}\simeq 22.83\ {\rm\mu Pa}\cdot{\rm s}$. It can be seen that the viscosity of He given above does not coincide with the value obtained via the ab initio formulation, namely $\widetilde{\mu}_{\rm ref;ab\ initio}^{\rm He}\simeq 19.91\ {\rm\mu Pa}\cdot{\rm s}$ Cencek et al. (2012); Sharipov and Strapasson (2013). Thus, in deriving the value of the Sutherland constant $\widetilde{S}$, we consider the values of the transport coefficients obtained in the framework of the ab initio calculations, as presented in Refs. Cencek et al. (2012) and Vogel et al. (2010) for Helium and Argon, respectively.

Since the maximum value of the temperature attained in Ref. Sharipov and Strapasson (2013) is $\lesssim 2\widetilde{T}_{\rm ref}=600\ {\rm K}$, we seek the values of $\widetilde{S}$ in Eq. (37) which best reproduce the reference data over the temperature range $300\ {\rm K}$ - $600\ {\rm K}$. A non-linear fit gives the following values for $\widetilde{S}$:

[TABLE]

With the above choice of parameters, the maximum relative deviation of the viscosity from the tabulated data is less than $0.8\%$ and $0.2\%$ for ${\rm He}$ and ${\rm Ar}$, respectively. The result of the fit is shown in Fig. 7.

The Sutherland model can be implemented by setting the non-dimensional relaxation time $\tau$ within the Shakhov model to:

[TABLE]

where $\delta$ is the rarefaction parameter employed in Ref. Sharipov and Strapasson (2013), while $S_{\rm He}\simeq 0.3101$ and $S_{\rm Ar}\simeq 0.5239$ after non-dimensionalization.

The profiles of $n$, $T$ and $u_{y}$ for $\delta=0.1$, $1$ and $10$ at wall velocity difference $2u_{w}=2\sqrt{2}$ (${\rm Ma}\simeq 2.2$) are shown in Fig. 8. Good agreement can be observed in general, with the largest discrepancies occuring in the temperature profile for $\delta=1$ and in the velocity profile at $\delta=0.1$. For $\delta=10$ and $1$, the simulations results were obtained with the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ model using $S=16$ lattice nodes stretched according to Eq. (46) with $A=0.98$ and the time step $\delta t=5\times 10^{-4}$. At $\delta=0.1$, the flow is within the transition regime and the simulation results obtained with $Q_{x}=7$ are no longer accurate within the S-model. We thus obtained the results at $\delta=0.1$ using the ${\rm HHLB}(6;30)\times{\rm HLB}(6;7)$ model with $S=16$, $A=0.98$ and $\delta t=2.5\times 10^{-4}$ (increasing $N_{x}$ to values higher than $6$ did not make any visible differences to the results). The simulation parameters are summarized in Table 3. We note that increasing $Q_{x}$ did not bring the LB simulation results closer to the DSMC profiles, as indicated in Fig. 9. This seems to indicate the fundamental limitation of the relaxation time approach, which fails to provide accurate results far within the transition regime.

In order to validate our model with the results reported in Ref. Sharipov and Strapasson (2013), we perform simulations at wall velocity differences $2u_{w}=2\sqrt{2}$ (${\rm Ma}\simeq 2.2$) and $2u_{w}=0.2\sqrt{2}$ (${\rm Ma}\simeq 0.22$), corresponding to $U=2v_{0}$ and $U=0.2v_{0}$ in Ref. Sharipov and Strapasson (2013), for various values $\delta\in\{0.01,0.1,1,10,20,40\}$ of the rarefaction parameter, covering the slip, transition and free molecular flow regimes. The validation is performed at a quantitative level based on the numerical results for the gradient $\nu$ of the velocity at the center of the channel (Table 4), the value $\Pi$ of the shear stress (Table 5) and the value $T_{0}$ of the temperature measured in the center of the channel (Table 6). In all cases, the LB results obtained using the above mentioned values of the simulation parameters were compared with the results from Ref. Sharipov and Strapasson (2013). The discrepancy is quantified by the relative error:

[TABLE]

where $A_{\rm LB}$ represents our simulation result and $A_{\rm ref}$ is the reference value from Ref. Sharipov and Strapasson (2013). The comparison is performed for $A\in\{\nu,\Pi,T_{0}\}$, where the values of the velocity gradient $\nu$, of the temperature $T_{0}$ (calculated in the channel center), as well as of the quantity $\Pi$ (derived from the non-diagonal stress $T_{xy}$), are introduced below.

We begin by considering the dimensionless velocity gradient defined in Ref. Sharipov and Strapasson (2013), which in our non-dimensionalization convention reads:

[TABLE]

The derivative appearing above is computed using Eq. (26). As can be seen in Table 4, our results are in very good agreement with those reported in Ref. Sharipov and Strapasson (2013) for both He and Ar for $\delta\geq 1$. These results were obtained using the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ model and have the absolute error bounded by $\pm 0.001$. At $\delta=0.1$ and $0.01$, in order to obtain LB results with the same $\pm 0.001$ absolute error, the quadrature order was raised up to $Q_{x}=30$ and $100$, respectively. For the convenience of our readers, these simulation parameters are summarized in Table 3. The effect of increasing $Q_{x}$ at small values of $\delta$ on the absolute error $\nu(Q_{x})-\nu_{\rm conv}$ computed with respect to the value $\nu_{\rm conv}$ obtained using the ${\rm HHLB}(6;100)\times{\rm HLB}(6;7)$ model on a grid with $S=32$ points ($A=0.98$) is shown in Fig. 10(a). The necessity to increase $Q_{x}$ as $\delta$ is decreased was also demonstrated in Fig. 2(b) in the context of the analysis of the linearized limit of the Boltzmann equation discussed in Sec. IV (a decreasing value of $\delta$ corresponds to an increasing value of $k$ in Fig. 2).

Next, we consider the non-dimensional quantity $\Pi$, defined as Sharipov and Strapasson (2013):

[TABLE]

where $\widetilde{n}_{\rm avg}$ represents the average density inside the channel. Although $\Pi$ should be constant in the stationary state of the Couette flow, small fluctuations of this quantity are always present across the channel. For this reason, $\Pi$ is computed using the average value of $T_{xy}$, obtained according to Eq. (29). Table 5 summarizes our results for $u_{w}\in\{\sqrt{2},0.1\sqrt{2}\}$. It can be seen that the relative error between the LB and DSMC results is less than $2\%$ for all tested values of the parameters. For consistency with the results reported for $\nu$ in Table 4, the models summarized in Table 3 were employed. However, we note that the values of $\Pi$ obtained using the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ model are within $0.2\%$ error with respect to the LB values inscribed in the table, as can be seen from Fig. 10(b). This observation is consistent with the results presented in Fig. 2(d) for the error of $T_{xy}$ as compared with the solution of the linearized Boltzmann equation.

Finally, the temperature $T_{0}$ in the channel center is obtained using the following formula:

[TABLE]

which is sixth order accurate with respect to the spacing $\delta\eta$ for even functions of $x$. In the above, $T_{s}$ corresponds to the temperature $T(x_{s})$ at point $x_{s}$, which is given by Eq. (47). Table 6 shows a comparison between the LB and the DSMC results for the temperature $T_{0}$ in the center of the channel, obtained for $u_{w}=\sqrt{2}$ and various values of $\delta$. As was the case for $\Pi$, the results shown in the table were obtained using the models summarized in Table 3, however, the results obtained using the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ model are within less than $0.1\%$ relative error with respect to the LB results obtained with the model ${\rm HHLB}(6;100)\times{\rm HLB}(6;7)$ on a finer grid ($S=32$, $A=0.98$), as can be seen in Fig. 10(c).

We end this section with a comment on the accuracy of the simulation results presented herein. The results reported in Fig. 8 and in Tables 4–6 were obtained on a grid with $S=16$ nodes, stretched according to Eq. (47) with $A=0.98$. The time step was set to $\delta t=5\times 10^{-4}$ for $Q_{x}=7$ and $\delta=2.5\times 10^{-4}$ for $Q_{x}=30$ and $100$, as summarized in Table 3. The results shown in Table 4 for $\nu$ have a maximum absolute error of $0.001$, while the results for $\Pi$ and $T_{0}$ shown in Tables 5 and 6 have relative errors of less than $0.1\%$. We checked that the effects of halving $\delta t$ or doubling $S$ were within these error bounds.

A comparison of the LB and the DSMC data for $\Pi$ and $T_{0}$ shown in Tables 5 and 6 together with the convergence analysis shown in Fig. 10, panels (b) and (c), reveals that the Shakhov model can provide accurate estimates of the DSMC results (relative errors below $2\%$) with a very modest cost, since the results obtained with the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ on a grid with $16$ nodes have less than $0.2\%$ error with respect to the results obtained using the ${\rm HHLB}(6;100)\times{\rm HLB}(6;7)$ model on a grid with $32$ nodes. In terms of absolute error, the values of the velocity gradient $\nu$, obtained using our LB models are within less than $5\%$ of the DSMC results. However, because the value of $\nu$ decreases as $\delta$ is decreased, the relative error becomes very large when $\delta\lesssim 0.1$. Furthermore, accurate estimates of $\nu$ using our LB models require the quadrature order $Q_{x}$ to be increased to large values as $\delta$ is decreased ($Q_{x}=100$ was employed at $\delta=0.01$).

VII Conclusion

In this paper, we studied the 3D Couette flow using the mixed quadrature lattice Boltzmann models introduced in Ref. Ambru\cbs and Sofonea (2016a), which employ the half-range Gauss-Hermite quadrature of order $Q_{x}$ on the axis perpendicular to the walls (the $x$ axis) and the full-range Gauss-Hermite quadrature of order $Q_{y}$ on the axis parallel to the flow direction (the $y$ axis). The third degree of freedom in the momentum space was removed through the analytic integration of the Boltzmann equation and the subsequent analysis was performed using reduced distribution functions.

We first validated our LB models in the low Mach number regime by comparing our simulation results with the benchmark results obtained in Ref. Jiang and Luo (2016) through a semi-analytic procedure applied to the linearized Boltzmann-BGK equation. To ensure that our simulations remained in the linearized regime, the wall velocity was set to a small value ($u_{w}=10^{-5}$). The validation was performed at the level of the velocity at the wall $u_{y}(L/2)$, the derivative of the velocity at the center of the channel $u^{\prime}_{y}(0)$, the non-diagonal component $T_{xy}$ of the stress tensor, and the half-channel mass flow rate $\dot{m}$. By employing a convergence test, we concluded that the minimum quadrature order required in order to achieve a given accuracy ($1\%$ error tolerance) must be increased as the Knudsen number $k$ is increased. Setting $Q_{x}=4$ ensures that the relative errors in $T_{xy}$ are less than $1\%$ up to $k=100$. For the remaining three quantities, the relative errors are below $1\%$ up to $k\lesssim 0.1$. Setting $Q_{x}=7$ preserves the $1\%$ error threshold up to $k\simeq 1$, while the relative error in $T_{xy}$ is decreased below $0.1\%$ up to $k=100$. Between $1\lesssim k\lesssim 100$, the quadrature order has to be increased dramatically in order to ensure that the relative errors in $u_{y}(L/2)$, $u^{\prime}_{y}(0)$ and $\dot{m}$ remain below $1\%$. This is partly due to the fact that the absolute values of these quantities decrease as $k$ is increased, such that maintaining a $1\%$ relative error entails an effective increase of the simulation accuracy. These convergence tests were performed employing a grid with $S=16$ points spanning half of the flow channel, stretched towards the bounding wall. When comparing the nonlinear part of the velocity profile with the data reported in Ref. Li et al. (2015), we found very good agreement after refining the grid to $S=64$ points for $k=0.03$ ($Q_{x}=4$ was sufficient here) and $S=32$ points for $k=0.1$ ($Q_{x}=7$), $k=1$ ($Q_{x}=21$) and $k=10$ ($Q_{x}=80$).

Next, we compared the LB profiles of the macroscopic quantities (particle number density, velocity, pressure tensor and heat flux) with the Direct Simulation Monte Carlo (DSMC) results for Maxwell molecules reported in Refs. Struchtrup and Torrilhon (2007, 2008); Torrilhon and Struchtrup (2008); Taheri et al. (2009); Schuetze . The LB profiles were obtained for ${\rm Kn}<1$ with the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ model, employing $Q_{x}=Q_{y}=7$ and an expansion up to $N_{x}=N_{y}=6$ of the equilibrium distribution with respect to the half-range (on the $x$ axis) and full range (on the $y$ axis) Hermite polynomials. At ${\rm Kn}=1$, the quadrature order on the $x$ axis was raised to $Q_{x}=11$. We found that the DSMC results for Knudsen numbers between $0.01\leq{\rm Kn}\leq 1$ and for wall velocities between $0.42\leq u_{w}\leq 2.1$ could be reasonably well recovered by employing the Shakhov collision term. We found deviations between our LB results and the DSMC data at ${\rm Kn}\gtrsim 0.25$, as well as when the wall velocity exceeded $u_{w}\gtrsim 1.68$.

Finally, we performed a comparison with the results for pure Helium and Argon obtained in Ref. Sharipov and Strapasson (2013) using an interaction model based on ab initio potentials. In order to match the ab initio transport coefficients, we implemented the Sutherland model and obtained the Sutherland constant by fitting the analytic expression for the viscosity to the tabulated data reported in Refs. Cencek et al. (2012) and Vogel et al. (2010). The relative errors of the viscosity obtained in the frame of the Sutherland model compared to the tabulated data are below $1\%$ for the temperature range relevant for the simulations considered in this paper ($300\ {\rm K}$ – $600\ {\rm K}$). At the level of the profiles for the density $n$, velocity $u_{y}$ and temperature $T$, the LB results are in very good agreement with the DSMC data when the rarefaction parameter $\delta$ satisfies $\delta>1$. At $\delta\lesssim 1$, visible deviations occur in the profile of $T$. Moreover, the velocity profile also deviates from the DSMC data at $\delta=0.1$. This observed discrepancy is due to the limitations of the relaxation time approximation of the collision integral, since increasing the quadrature order does not bring our FDLB results closer to the DSMC data. A quantitative analysis at the level of the temperature $T_{0}$ at the channel center and of the shear pressure $\Pi$ shows that the deviations of our FDLB models from the DSMC results are within a few percent. The velocity gradient $\nu$ at the channel center presents an increasing relative error as $\delta$ is decreased, which may also be due to the fact that $\nu$ decreases towards [math] as $\delta$ is decreased. In terms of an absolute error, the LB results for $\nu$ still remain within a few percent of the DSMC data (for the Argon gas at $\delta=0.01$ and $u_{w}=\sqrt{2}$, $\nu_{\rm DSMC}=0.048$, while the FDLB result is $\nu_{\rm FDLB}=0.027$).

The analysis presented in this paper indicates that the solution of the S-model equation seems to be within a few percent of the DSMC results in the context of the Couette flow for velocities up to $2.1$ and for values of the rarefaction parameter $\delta$ down to $0.01$ (values of the Knudsen number ${\rm Kn}$ up to $\simeq 70$).

We finish this paper by noting that the ${\rm HHLB}(6;7)\times{\rm HLB}(6;7)$ model, which employs only $2Q_{x}Q_{y}=98$ distinct velocities, can be used to obtain a very good estimate (within $1\%$ relative error) of the solution of the Shakhov model equation at the level of the temperature in the channel center and non-diagonal component of the stress tensor. We thus conclude that the use of half-range quadratures is an essential ingredient when considering the channel flow of rarefied gases.

Acknowledgements

This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, Project No. PN-II-ID-PCE-2011-3-0516. The authors express their gratitude to Professor Henning Struchtrup (Department of Mechanical Engineering, University of Victoria, Canada) for kindly providing the DSMC results for Maxwell molecules used for the validation of the models introduced in this paper. The authors are grateful to Prof. Li-Shi Luo (Department of Mathematics and Statistics, Old Dominion University, Norfolk, VA) for suggesting the validation of their results by considering the non-linear part of the velocity profile in the linearized regime. The authors are indebted to Prof. Felix Sharipov (Universidade Federal do Paraná, Curitiba, Brazil) for suggesting the validation tests in the context of realistic, ab-initio potentials and for kindly providing DSMC results for comparison.

Appendix A Numerical scheme

The simulation results presented in this paper were obtained using an explicit third order total variation diminishing (TVD) Runge-Kutta (RK-3) time marching procedure Gottlieb and Shu (1998); Henrick et al. (2005); Shu and Osher (1988); Trangenstein (2007), as described in Subsec. A.1. In order to increase the simulation efficiency, we follow Refs. Guo and Zhao (2003); Mei and Shyy (1998) and employ a grid stretching algorithm to increase the grid resolution in the vicinity of the solid wall, as described in Subsec. A.2. The fifth-order weighted essentially non-oscillatory (WENO-5) scheme Gan et al. (2011); Jiang and Shu (1996) employed for the advection is presented in Subsec. A.3. Finally, the implementation of the diffuse reflection and bounce-back boundary conditions is discussed in Subsec. A.4.

A.1 Time stepping

In order to implement the time stepping algorithm, it is convenient to cast the Boltzmann equation (11) in the following form:

[TABLE]

where $F\in\{\phi,\chi\}$ is any of the two reduced distributions introduced in Sec. III, while $\mathbb{S}_{F}$ is given in Eq. (16).

Let us consider a discretization of the time coordinate using equal time steps $\delta t$, such that $t_{\ell}=\ell\,\delta t$. The distribution functions at time step $\ell$ can be written as $F_{\ell}\equiv F(t_{\ell})$. The value of $F_{\ell+1}$ can be obtained using the third-order Runge-Kutta TVD method introduced in Ref. Shu and Osher (1988), using two intermediate steps, as follows:

[TABLE]

The Butcher tableau Butcher (2008) for the above scheme is summarized in Table 7.

A.2 Grid stretching

As highlighted in Refs. Guo and Zhao (2003); Mei and Shyy (1998), a finer mesh is needed in the vicinity of solid boundaries as compared to the bulk regions of the flow in order to capture the Knudsen layer effects. This can be achieved by performing a standard grid-stretching procedure and in this paper, we follow Ref. Busuioc and Ambru\cbs and characterize the refined mesh using the non-dimensional parameter $\eta$ as follows:

[TABLE]

where $0\leq\eta\leq{\rm arctanh}\,A$ and $0<A<1$ controls the stretching such that when $A\rightarrow 0$, the grid becomes equidistant with respect to $x$, while as $A\rightarrow 1$, the grid points accumulate towards the right boundary.

The flow domain is discretized using $S$ equidistant values of $\eta$, namely:

[TABLE]

where the points with $1\leq s\leq S$ lie within the flow domain. The value $A=0.98$ was employed for all simulations presented in this paper.

A.3 Advection

The spatial derivative occuring in Eq. (11) can be approximated by considering an equidistant grid with respect to the $\eta$ coordinate (47):

[TABLE]

The flux $\mathcal{F}_{s+1/2}$ corresponds to the interface between the cells centered on $\eta_{s}$ and $\eta_{s+1}$, while the coordinates $x_{s\pm 1/2}$ of these interfaces are obtained by substituting $\eta=\eta_{s}\pm\delta\eta/2$ in Eq. (47). The fluxes are computed using the WENO-5 algorithm Gan et al. (2011); Ambru\cbs and Blaga (2018); Busuioc and Ambru\cbs ; Busuioc et al. . For the case when the advection velocity $p_{x}/m$ is positive, the flux is given by

[TABLE]

where $\mathcal{F}^{q}_{s+1/2}$ ($q=1,2,3$) are interpolating functions, which can be computed as follows:

[TABLE]

The weighting factors $\overline{\omega}_{q}$ are given by:

[TABLE]

where $\delta_{q}\in\{0.1,0.6,0.3\}$ are the ideal weights. The indicators of smoothness $\sigma_{q}$ can be computed using:

[TABLE]

It is customary in numerical algorithms to add a small quantity $\varepsilon\simeq 10^{-6}$ to $\sigma_{q}$ in order to avoid division by zero. This operation can have side effects which depend on the magnitude of the advected quantity $F$, as discussed in Ref. Henrick et al. (2005). In order to avoid such side effects, $\overline{\omega}_{q}$ is computed directly from Table 8 in the limiting cases when one or more of the $\sigma_{q}$ functions vanish.

A.4 Boundary conditions

The Couette flow considered in this paper is symmetric with respect to the channel centerline, thus allowing the simulation domain to be reduced to only the right half of the channel, such that $0\leq x\leq L/2$. The symmetry condition of the Couette flow is immediately achieved when bounce-back boundary conditions are implemented on the centerline at $x=0$. The gas-wall interaction is modeled using diffuse reflection boundary conditions Ambru\cbs and Sofonea (2016a, b); Sofonea (2009), which are implemented at $x=L/2$. In order to apply the fifth order WENO scheme described in the previous subsection, the simulation domain must be extended on both sides through the addition of three ghost nodes. Let the pair of indices $ij$ ( $1\leq i\leq 2Q_{x}$, $1\leq j\leq Q_{y}$) label the momentum vector corresponding to each discrete population.

For the bounce-back condition at $x=0$, the following procedure is performed to define the particle populations in the ghost nodes. Let $F_{1;ij}$, $F_{2;ij}$ and $F_{3;ij}$ be the population of particles of momentum ${\bm{p}}_{ij}=(p_{x,i},\,p_{y,j})$, located in the first three nodes of the simulation domain near the channel centerline. These nodes are counted in the positive (right) direction of the $x$ axis. The first three ghost nodes located at the left of the channel centerline and counted in the negative direction of the $x$ axis, have the populations $F_{0;ij}$, $F_{-1;ij}$ and $F_{-2;ij}$, respectively. To implement the bounce-back condition, these ghost populations are related to the populations in the simulation domain according to

[TABLE]

where the indices $\widetilde{\imath}$ ($\widetilde{\jmath}$) refer to the components $p_{x,\widetilde{\imath}}$ ($p_{y,\widetilde{\jmath}}$) defined through:

[TABLE]

Let $S$ denote the last (rightmost) node located in the flow domain. The first three ghost nodes outside the right boundary will be denoted $S+1$, $S+2$, $S+3$. On the right boundary, the diffuse reflection concept should be imposed. According to this concept, the flux of particles coming from the ghost nodes is Maxwellian and equals $\phi^{\rm(eq)}_{w;ij}p_{x,i}/m$, where $\phi^{\rm(eq)}_{w;ij}$ is defined by Eq. (15). In the frame of the WENO scheme, this can be exactly achieved when Busuioc and Ambru\cbs ; Busuioc et al. :

[TABLE]

where $F^{\rm(eq)}_{w;ij}=\phi^{\rm(eq)}_{w;ij}$ when $F_{s;ij}$ refers to $\phi_{s;ij}$, while $F^{\rm(eq)}_{w;ij}=\chi^{\rm(eq)}_{w;ij}=T_{w}\phi^{\rm(eq)}_{w;ij}$ in the case when $F_{s;ij}$ refers to $\chi_{s;ij}$. Since Eqs. (45) and (48) cannot be used in the nodes $s\in\{S+1,S+2\}$ when $p_{x,i}>0$, the corresponding functions $F_{s;ij}$, which describe the particles travelling rightwards, are extrapolated at every time step by a quadratic procedure:

[TABLE]

The wall density $n_{w}$, which is required in order to construct $F^{\rm(eq)}_{w;ij}\in\{\phi^{\rm(eq)}_{w;ij},\chi^{\rm(eq)}_{w;ij}\}$, is thereafter obtained by imposing mass conservation on the right wall:

[TABLE]

where $\Phi_{S+1/2;ij}$ is the flux (49) corresponding to the reduced distribution $\phi_{ij}$ through the interface between the last fluid cell and the first ghost node.