2.1. The simplified kinetic model for van der Waals fluids

Following our previous works (Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020; Su et al. Reference Su, Gibelli, Li, Borg and Zhang2023), we expand the collision operator (2.4) into a Taylor series near $\boldsymbol {r}$ and retain up to the second-order terms as shown below:

(2.11)\begin{equation} J_E(f,f) = \chi\,J^{(0)}(f,f) + J^{(1)}(f,f) + J^{(2)}(f,f), \end{equation}

with

(2.12)\begin{equation} \left.\begin{aligned} J^{(0)}(f,f) & = \sigma^2 \iint (f'f_1'-ff_1) \boldsymbol{g}\boldsymbol{\cdot}\boldsymbol{k}\,{\rm{d}}\boldsymbol{k}\, {\rm d}\boldsymbol{\xi}_1, \\ J^{(1)}(f,f) & = \sigma^3\chi \iint \boldsymbol{k} \boldsymbol{\cdot} (f'\,\boldsymbol{\nabla} f_1' + f\,\boldsymbol{\nabla} f_1) \boldsymbol{g} \boldsymbol{\cdot} \boldsymbol{k}\,{\rm{d}} \boldsymbol{k}\,{\rm d}\boldsymbol{\xi}_1\\ & \quad + \frac{\sigma^3}{2}\iint (\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}\chi) (f'f_1' + ff_1) \boldsymbol{g} \boldsymbol{\cdot}\boldsymbol{k}\,{\rm{d}}\boldsymbol{k}\,{\rm d} \boldsymbol{\xi}_1,\\ J^{(2)}(f,f) & = \frac{\sigma^4}{2}\,\chi \iint \boldsymbol{kk}: (f'\,\boldsymbol{\nabla} \boldsymbol{\nabla} f_1' - f\,\boldsymbol{\nabla}\boldsymbol{\nabla} f_1) \boldsymbol{g} \boldsymbol{\cdot} \boldsymbol{k}\,{\rm{d}}\boldsymbol{k}\,{\rm d}\boldsymbol{\xi}_1\\ & \quad+ \frac{\sigma^4}{2}\iint (\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\nabla} \chi) [\boldsymbol{k} \boldsymbol{\cdot}(f'\,\boldsymbol{\nabla} f_1' - f\,\boldsymbol{\nabla} f_1)]\boldsymbol{g}\boldsymbol{\cdot}\boldsymbol{k}\,{\rm{d}}\boldsymbol{k}\,{\rm d}\boldsymbol{\xi}_1\\ & \quad+ \frac{\sigma^4}{8}\iint (\boldsymbol{kk} : \boldsymbol{\nabla}\boldsymbol{\nabla} \chi) (f' f_1' - f f_1)]\boldsymbol{g} \boldsymbol{\cdot}\boldsymbol{k}\,{\rm{d}}\boldsymbol{k}\,{\rm d} \boldsymbol{\xi}_1, \end{aligned}\right\} \end{equation}

where all the quantities are evaluated at the position $\boldsymbol {r}$, and $J^{(0)}(f,f)$ is the Boltzmann collision operator for dilute gases, which can be simplified further by kinetic models. Here, we choose the Shakhov model (Shakhov Reference Shakhov1968), which is written as

(2.13)\begin{equation} J^{(0)}(f,f)\simeq J_s={-}\frac{1}{\tau_s}\left[(f-f^{eq})-f^{eq}(1-Pr)\,\frac{\boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{Q}_k}{5p_0RT}\left(\frac{c^2}{RT}-5\right)\right], \end{equation}

where $\tau _s$ is the relaxation time, $Pr$ is the Prandtl number, $p_0=nk_BT$ is the EoS for ideal gases, $R=k_B/m$ is the specific gas constant, and $f^{eq}$ is the Maxwellian distribution function, which reads as

(2.14)\begin{equation} f^{eq}=n\left(\frac{m}{2{\rm \pi} k_BT}\right)^{{3}/{2}}\exp\left(-\frac{mc^2}{2k_BT}\right). \end{equation}

The terms $J^{(1)}(f,f)$ and $J^{(2)}(f,f)$ describe the dense gas effect arising from increasing density. Considering that (i) for dilute gases far from equilibrium, the density terms $J^{(1)}(f,f)$ and $J^{(2)}(f,f)$ are negligible, and the non-equilibrium effect can be captured by the Shakhov model (2.13), (ii) for the gases of large densities, the density terms $J^{(1)}(f,f)$ and $J^{(2)}(f,f)$ become important, where the gas mean free path should be small as $\lambda \propto 1/n$, implying that the gases are not far from equilibrium, and (iii) for gases not far from equilibrium, the equilibrium distribution function $f^{eq}$ is the leading part of the distribution function $f$, further simplifications can be made on $J^{(1)}(f,f)$ and $J^{(2)}(f,f)$ by replacing the velocity distribution functions therein with their corresponding equilibrium distribution functions, leading to the following two terms (Rangel-Huerta & Velasco Reference Rangel-Huerta and Velasco1996; Kremer Reference Kremer2010; Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020; Su et al. Reference Su, Gibelli, Li, Borg and Zhang2023):

(2.15)\begin{equation} J^{(1)}(f,f)\simeq \mathcal{I}^{(1)}={-}nV_0 \chi f^{eq}\left\{ \begin{array}{l} \boldsymbol{c}\boldsymbol{\cdot} \left[\boldsymbol{\nabla} \ln(n^2\chi T)+\dfrac{3}{5} \left(\mathcal{C}^2-\dfrac{5}{2}\right)\boldsymbol{\nabla} \ln T \right]\\ \quad{}+\dfrac{2}{5}\left[2\mathcal{CC}:\boldsymbol{\nabla} \boldsymbol{u} +\left(\mathcal{C}^2-\dfrac{5}{2}\right) \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \right] \end{array}\right\} \end{equation}

and

(2.16)\begin{equation} J^{(2)}(f,f)\simeq\mathcal{I}^{(2)}=\boldsymbol{\nabla}\boldsymbol{\cdot}\left[f^{eq}\,\frac{\mu_B}{p_0}\,(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u})\left(\mathcal{C}^2-\frac{3}{2}\right)\boldsymbol{c}\right]+\mathcal{R}, \end{equation}

where $\mathcal {C}=(m/2k_BT)^{1/2}\boldsymbol c$ is the non-dimensional peculiar velocity, and $\mu _B$ is the bulk viscosity. Here, $\mathcal {R}$ is a second-order quantity that has no contribution to the transfer of mass, momentum and energy, so it can be ignored hereafter in the kinetic model. Since $n \to 0$ leads to $\mu _B\to 0$ for monatomic gases, the dense gas terms $\mathcal {I}^{(1)}$ and $\mathcal {I}^{(2)}$ become negligible in the dilute gas limit. In this case, the present kinetic model reduces to the Shakhov model.

It should be noted that the bulk viscosity $\mu _B$ appears when calculating the pressure tensor and heat flux from the Enskog equation using the Chapman–Enskog analysis (Ferziger & Kaper Reference Ferziger and Kaper1972; Chapman & Cowling Reference Chapman and Cowling1990), but is absent in the previous simplified kinetic models (Luo Reference Luo1998; He & Doolen Reference He and Doolen2002; Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020; Takata et al. Reference Takata, Matsumoto and Hattori2021). Although it is a small quantity involved in the second-order term of the Taylor series (Rangel-Huerta & Velasco Reference Rangel-Huerta and Velasco1996; Kremer Reference Kremer2010), it is important in many applications (Jaeger, Matar & Müller Reference Jaeger, Matar and Müller2018), such as sound attenuation and shock wave propagation, where gases undergo strong compression or expansion (Hoover et al. Reference Hoover, Evans, Hickman, Ladd, Ashurst and Moran1980a,Reference Hoover, Ladd, Hickman and Holianb).

For simplicity, the pair correlation function $\chi$ in (2.11) can be absorbed into the relaxation time $\tau _s$ in (2.13). The final evolution equation of the kinetic model for the van der Waals fluids can be written as

(2.17)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol{\xi} \boldsymbol{\cdot}\boldsymbol{\nabla} f + \frac{\boldsymbol{F}_{ext}+\boldsymbol{F}_{att}}{m} \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{\xi}}f = J^{(0)}_s + \mathcal{I}^{(1)} + \mathcal{I}^{(2)}, \end{equation}

where

(2.18)\begin{equation} J^{(0)}_s={-}\frac{1}{\tau}\left[(f-f^{eq})-f^{eq}(1-Pr)\,\frac{\boldsymbol{c}\boldsymbol{\cdot}\boldsymbol{Q}_k}{5p_0RT}\left(\frac{c^2}{RT}-5\right)\right], \end{equation}

with the relaxation time $\tau =\tau _s/\chi$. The attractive part of the LJ potential, i.e. $\phi _{att}=-4\epsilon (\sigma /r)^{6}$, is chosen to simulate the molecular attraction in the mean-field force term (2.6). It should be emphasised that (2.17) is second-order accurate in the Taylor series of the Enskog collision operator (2.4), with omitted second-order quantities that have no contribution to mass, momentum and energy transfer. It should be noted that although the collisional terms $J^{(0)}_s$, $\mathcal {I}^{(1)}$ and $\mathcal {I}^{(2)}$ are derived from the Enskog collision operator (2.4), the kinetic model (2.17) is not restricted to hard-sphere molecules as the transport coefficients can be corrected to account for the influence of intermolecular potentials. In the following subsections, we will demonstrate the thermodynamic consistency of our kinetic model and how to obtain correct transport coefficients.

2.2. The hydrodynamic equations and relaxation time

Using the Chapman–Enskog expansion (see Appendix A for the details), the following hydrodynamic equations can be obtained:

(2.19)\begin{align} \left.\begin{gathered} \frac{\partial \rho}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho \boldsymbol{u})=0,\\ \frac{\partial (\rho \boldsymbol{u})}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} \{\rho \boldsymbol{uu}+[p-\mu_B(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})] \boldsymbol{\mathsf{U}}-2\mu_s\mathring{\boldsymbol{\mathsf{S}}}-\boldsymbol{\mathsf{K}}\}-n\boldsymbol{F}_{ext}=0,\\ \frac{\partial (\rho E)}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} \{\rho E \boldsymbol{u}-\kappa\,\boldsymbol{\nabla} T+ [p-\mu_B(\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u})]\boldsymbol{\mathsf{U}}\boldsymbol{\cdot}\boldsymbol{u} -(2\mu_s\mathring{\boldsymbol{\mathsf{S}}}+\boldsymbol{\mathsf{K}})\boldsymbol{\cdot}\boldsymbol{u}\} - n\boldsymbol{F}_{ext}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gathered}\right\} \end{align}

where $E=(3RT+\boldsymbol {u}^2)/2$ is the combination of internal and kinetic energies per unit mass of monatomic gases, and $\mathring {\boldsymbol{\mathsf{S}}}$ and $\boldsymbol{\mathsf{K}}$ are the rate-of-shear and capillary tensors, respectively, which can be found in Appendix A. According to (A20), the relaxation time $\tau$ and the Prandtl number $Pr$ can be obtained as

(2.20)\begin{equation} \left.\begin{gathered} \tau= \dfrac{\mu_s}{nk_BT}\,\dfrac{1}{1+\dfrac{2}{5} nV_0\chi}, \\ Pr = \dfrac{5k_B}{2m}\,\dfrac{{1+\dfrac{3}{5} nV_0\chi}}{{1+\dfrac{2}{5} nV_0\chi}}\,\dfrac{\mu_s}{\kappa}. \end{gathered}\right\} \end{equation}

The proper assignment of the relaxation time $\tau$ and Prandtl number $Pr$ is a key requirement for the thermodynamic consistency of the kinetic model, which depends on the appropriate determination of the shear viscosity $\mu _s$ and thermal conductivity $\kappa$ of the fluids, to be discussed in § 2.3.

It should be noted that the hydrostatic pressure $p$ in (2.19) satisfies the van-der-Waals-type EoS, where both the volume exclusion and the intermolecular attraction are considered. The specific form of the EoS depends on the choice of the pair correlation function $\chi$. If we choose $\chi =1/(1-nV_0)$, then the hydrostatic pressure (A22) recovers the exact van der Waals EoS (1.1). However, the shielding effect of the gas molecules is not taken into account by this choice. As the density inhomogeneity may become significant in a nano-confined system, we use the Fischer–Methfessel method to consider the density dependence of the pair correlation function (Fischer & Methfessel Reference Fischer and Methfessel1980), and the Carnahan–Starling EoS to consider the shielding effect (Carnahan & Starling Reference Carnahan and Starling1969). Therefore, the pair correlation function can be written as

(2.21)\begin{equation} \chi = \chi[\boldsymbol{r}_1, \boldsymbol{r}_2\mid n(\boldsymbol{r}_1 \to \boldsymbol{r}_2)] \simeq \chi(\bar{n})=\frac{1-0.5\eta}{(1-\eta)^3},\quad \eta=0.25\bar{n}V_0, \end{equation}

where $\bar {n}=\int w(\boldsymbol {r}')\,n(\boldsymbol {r}+\boldsymbol {r'})\,\mathrm {d}\boldsymbol {r'}$ is the local average density, with $w(\boldsymbol {r})$ being a weight function (Tarazona Reference Tarazona1985; Vanderlick, Scriven & Davis Reference Vanderlick, Scriven and Davis1989). Substituting (2.21) into (A22), we can get the hydrostatic pressure $p$ satisfying the EoS

(2.22)\begin{equation} p=nk_BT\,\frac{1+\eta+\eta^2-\eta^3}{(1-\eta)^3}-an^2. \end{equation}

Other cubic equations of state, such as the Soave–Redlich–Kwong and Peng–Robinson types, can also be recovered by choosing appropriately the pair correlation function $\chi$ using the above approach. This hydrostatic pressure, which is recovered from the non-equilibrium transport equation, is consistent with the thermodynamic theory for the equilibrium state, demonstrating that our kinetic model (2.17) is thermodynamically consistent.

2.3. Transport coefficients for van der Waals fluids

The transport coefficients in (2.8) and (2.9) are obtained through the first-order Chapman–Enskog expansion of the Enskog equation, and include both the kinetic and collisional contributions. For simplicity, the derivation of (2.8) and (2.9) was based on the hard-sphere molecules, i.e. all intermolecular collisions are rigid and elastic. To improve accuracy of the predictions for real gases, the molecular dimensions are assumed to change with temperature, i.e. a higher temperature leads to a smaller molecular diameter, which has been adopted widely in MET (Hanley et al. Reference Hanley, McCarty and Cohen1972), kinetic reference theory (Karkheck & Stell Reference Karkheck and Stell1981) and other models (Guo, Zhao & Shi Reference Guo, Zhao and Shi2005; Guo et al. Reference Guo, Zhao and Shi2006; Shan et al. Reference Shan, Wang, Zhang and Guo2020). This modification accounts for the softness of molecules during the collision, but the effect of the gas molecular attraction on transport coefficients can still not be considered properly. In contrast to the Enskog equation, which describes the dynamics of hard-sphere gases, no molecular potential model appears explicitly in our kinetic model (2.17). Instead, the intermolecular potential, including molecular attraction for real gases, is included in the transport coefficients.

For different molecular potential models (Chapman & Cowling Reference Chapman and Cowling1990), the shear viscosity of real dilute gases can be written as

(2.23a,b)\begin{equation} \mu_0=\frac{\mu_0^{hs}}{\varOmega^{(2,2)}},\quad \mu_0^{hs}=\frac{5}{16\sigma^2}\sqrt{\frac{mk_BT}{\rm \pi}}, \end{equation}

where $\mu _0^{hs}$ is the viscosity of dilute gases of hard-sphere molecules, and $\varOmega ^{(2,2)}$ is the reduced collision integral depending on the intermolecular potential, which accounts for the effect of gas molecular attraction on viscosity and is difficult to obtain theoretically. Neufeld, Janzen & Aziz (Reference Neufeld, Janzen and Aziz1972) proposed an empirical form of the integral that performs well (with error less than 0.1 %) in the temperature range $0.3\leq \hat {T}\leq 100$, with $\hat {T}=k_BT/\epsilon$, which can be written as

(2.24)\begin{equation} \varOmega^{(2,2)}=\frac{c_1}{\hat{T}^{c_2}} + c_3\exp(c_4\hat{T})+c_5\exp(c_6\hat{T})+c_7\hat{T}^{c_8}\sin(c_9\hat{T}^{c_{10}}+c_{11}), \end{equation}

with corresponding coefficients given in table 1.

Table 1. Coefficients for calculating the transport integral in (2.24).

To obtain the shear viscosity and thermal conductivity of the van der Waals fluids, we use the method proposed by Chung, Lee & Starling (Reference Chung, Lee and Starling1984) and Chung et al. (Reference Chung, Ajlan, Lee and Starling1988), which is based on the kinetic theory and experimental correlation. For convenience, we convert the original expression to the following form where the shear viscosity can be calculated as

(2.25)\begin{equation} \mu_s = \mu_0^{hs}\left(\frac{F_AF_B}{\varOmega^{(2,2)}}+F_C\right), \end{equation}

with

(2.26a)\begin{gather} F_A=1-0.2756\omega, \end{gather}

(2.26b)\begin{gather}F_B=\frac{1}{G_v}+A_6\eta, \end{gather}

(2.26c)\begin{gather}F_C=\frac{1}{\hat{T}^{{1}/{2}}}\,A_7\eta^2G_v\exp\left(A_8+\frac{A_9}{\hat{T}}+\frac{A_{10}}{\hat{T}^2}\right), \end{gather}

(2.26d)\begin{gather}G_v=\frac{(A_1/\eta)[1-\exp({-}A_4\eta)]+A_2\chi\exp(A_5\eta)+A_3\chi}{A_1A_4+A_2+A_3}, \end{gather}

where $F_A$ accounts for the effect of the acentric of molecules, with $\omega$ being the acentric factor, and $F_B$ and $F_C$ account for the dependence of viscosity on gas density. For monatomic gases, the acentric factor is $\omega =0$, so $F_A=1$. The coefficients $A_1$–$A_9$ can be calculated by

(2.27)\begin{equation} A_i=a_0(i)+a_1(i)\,\omega, \end{equation}

with the corresponding constants listed in table 2.

Table 2. Coefficients for calculating the viscosity of van der Waals fluids in (2.27).

Similarly, the thermal conductivity of dilute gases can be calculated as

(2.28a,b)\begin{equation} \kappa_0=\frac{\kappa_0^{hs}}{\varOmega^{(2,2)}},\quad \kappa_0^{hs}=\frac{75k_B}{64m\sigma^2}\sqrt{\frac{mk_BT}{\rm \pi}}. \end{equation}

The thermal conductivity of van der Waals fluids at large densities can be calculated as

(2.29)\begin{equation} \kappa=\kappa_0^{hs}\left(\frac{F_P F_A F_D}{\varOmega^{(2,2)}}+F_E\right), \end{equation}

with

(2.30a)\begin{gather} F_D=\frac{1}{G_t}+B_6\eta, \end{gather}

(2.30b)\begin{gather}F_E=0.8906B_7\eta^2G_t, \end{gather}

(2.30c)\begin{gather}G_t=\frac{(B_1/\eta)[1-\exp({-}B_4\eta)]+B_2\chi\exp(B_5\eta)+B_3\chi}{B_1B_4+B_2+B_3}, \end{gather}

where $F_P$ accounts for the polyatomic effect on thermal conductivity, which is unity for monatomic gases, $F_D$ and $F_E$ account for the dependence of thermal conductivity on density, and the coefficients $B_1$–$B_7$ can be calculated from

(2.31)\begin{equation} B_i=b_0(i)+b_1(i)\,\omega \end{equation}

using the constants listed in table 3.

Table 3. Coefficients to calculate the thermal conductivity of van der Waals fluids in (2.31).

Since the correlated density-dependent functions are introduced to extend the Enskog model (2.8) and (2.9) for real gases by taking gas molecular attraction into account, the modified Enskog model (2.25) and (2.29) is referred to as the correlated Enskog model in this study. Once the shear viscosity $\mu _s$ and thermal conductivity $\kappa$ are calculated from (2.25) and (2.29), respectively, the relaxation time $\tau$ and Prandtl number $Pr$ can be determined through (2.20).

One last parameter that needs to be determined is the bulk viscosity $\mu _B$, which was derived for the hard-sphere fluids as

(2.32)\begin{equation} \mu_B^{hs}=(nV_0)^2\chi \mu_0^{hs}. \end{equation}

This equation overestimates the bulk viscosity of dense LJ fluids according to Hoover et al. (Reference Hoover, Evans, Hickman, Ladd, Ashurst and Moran1980a) and Borgelt, Hoheisel & Stell (Reference Borgelt, Hoheisel and Stell1990). The overestimation is inherent in the calculation of the shear viscosity and thermal conductivity of the Enskog model given by (2.8) and (2.9), as these two transport coefficients are a combination of kinetic and collisional contributions. Taking the shear viscosity (2.8) as an example, (2.8) can be rewritten as

(2.33)\begin{equation} \mu^{hs} = \mathop{ \underbrace{\frac{\mu_0^{hs}}{\chi}\left(1+\frac{2}{5}\,nV_0\chi\right)}}\limits_{\mu_k} + \mathop{\underbrace{\frac{\mu_0^{hs}}{\chi}\left(1+\frac{2}{5}\,nV_0\chi\right)\frac{2}{5}\,nV_0\chi+\frac{3}{5}\,\mu_B^{hs}}}\limits_{ \mu_c}, \end{equation}

where $\mu _k$ and $\mu _c$ are the kinetic and collisional contributions to the shear viscosity, respectively. An overestimation of the bulk viscosity in $\mu _c$ will naturally lead to an overestimation of the shear viscosity, especially at large densities where $\mu _c$ dominates. This explains that the transport coefficients of the Enskog model at large densities are not accurate.

Gray & Rice (Reference Gray and Rice1964) proposed an explicit formula for the bulk viscosity, suggesting that the bulk viscosity consists of three parts: the hard-core collision part $\mu _B^{hs}$, the long-range attractive part $\mu _B^{att}$, and the cross (intermediate) part $\mu _B^{crs}$ between hard-core collision and long-range attraction, namely $\mu _B=\mu _B^{hs}+\mu _B^{att}+\mu _B^{crs}$. There are conflicting explanations for this formula. Madigosky (Reference Madigosky1967) stated that the cross part $\mu _B^{crs}$ is negligible when $\hat {T}>1$ and the long-range attractive part satisfies $\mu _B^{att} \propto \rho ^2$, and $\mu _B^{att}$ is always positive. On the contrary, Collings & Hain (Reference Collings and Hain1976) found that the cross part $\mu _B^{crs}$ cannot be neglected, and the long-range attractive part $\mu _B^{att}$ can be negative at large densities, which is consistent with the fact that the Enskog prediction of the transport coefficients is much larger than the experimental values at large densities, where the contribution of the long-range molecular attraction to the bulk viscosity is ignored.

A two-parametric function has recently been proposed by Chatwell & Vrabec (Reference Chatwell and Vrabec2020) to calculate the bulk viscosity, which is in good agreement with the experimental data and the MD simulation results at ultra-low temperature and ultra-high density conditions. However, it may become problematic when the density reduces or temperature increases, as non-physical bulk viscosity would appear. Overall, the bulk viscosity for dense monatomic gases needs further investigation. Here, we adopt an empirical approach (Hoover et al. Reference Hoover, Ladd, Hickman and Holian1980b) to calculating the bulk viscosity, which considers the effect of attraction between gas molecules as

(2.34)\begin{equation} \mu_B =nV_0y\left(\frac{\epsilon}{k_BT}\right)^{{1}/{12}}\mu_0^{hs}, \end{equation}

with

(2.35a)\begin{gather} y=2.722x+3.791x^2 +2.495x^3-1.131x^5, \end{gather}

(2.35b)\begin{gather}x=0.477465nV_0\left(\frac{\epsilon}{k_BT}\right)^{{1}/{4}}. \end{gather}

To be consistent with the shear viscosity and thermal conductivity, we refer to (2.34) as the correlated Enskog model since the effect of gas molecular attraction is included.

3.1. Model analysis and comparison

The pair correlation function plays an essential role in the MET. A key requirement for determining the pair correlation function is that $\chi \to 1$ as $n \to 0$, so that the Enskog-type equation for dense gases reduces to the Boltzmann-type equation in the dilute gas limit. However, the MET does not satisfy this requirement when the van der Waals pressure is chosen, as shown by (1.4), which makes the MET inaccurate in capturing the effect of the long-range molecular attraction. A temperature-dependent diameter (Hanley et al. Reference Hanley, McCarty and Cohen1972) can be employed to correct this error, which relates the co-volume $V_0$ with the second virial coefficient $B$ through $V_0 =B+T\,\mathrm {d}B/\mathrm {d}T$, and leads directly to the following EoS for real gases:

(3.2a,b)\begin{equation} p=Znk_BT, \quad Z=1+n\chi\left(B+T\,\frac{\mathrm{d}B}{\mathrm{d}T}\right), \end{equation}

where $Z$ is the compressibility factor. Clearly, the real gas EoS recovers the ideal gas EoS as the compressibility factor $Z \to 1$ when $n \to 0$. However, the compressibility factor $Z$ may be less than unity near the critical temperature (Mahmoud Reference Mahmoud2014), which means that the pair correlation function $\chi$ may be negative in (3.2a,b) as both $n>0$ and $V_0 =B+T\,\mathrm {d}B/\mathrm {d}T>0$, thus leading to negative shear viscosity and thermal conductivity, as can be seen from (2.8) and (2.9), which is physically inappropriate. Therefore, the MET is not suitable for modelling gas dynamics of the van der Waals fluids.

As shown in figure 2, the Enskog prediction overestimates the shear viscosity and thermal conductivity at large densities. The assumption of a state-dependent diameter (2.7) attenuates this overestimation at low temperatures (see figure 2a), but leads to an overestimation of the shear viscosity at low densities (see figure 2b). Overall, the correlated Enskog model agrees well with the experimental data for a wide range of temperatures and densities, particularly for shear viscosity and thermal conductivity, which indicates the importance of molecular attraction.

Figure 2. Comparison of transport coefficients of argon: (a,b) for the shear viscosity at $T= 173.0$ K and 298.0 K, respectively, with the experimental data from Haynes (Reference Haynes1973); (c,d) for the thermal conductivity at $T= 298.15$ K and 348.15 K, respectively, with the experimental data from Michels, Sengers & Van de Klundert (Reference Michels, Sengers and Van de Klundert1963); and (e,f) for the bulk viscosity with the experimental data from Malbrunot et al. (Reference Malbrunot, Boyer, Charles and Abachi1983) and Madigosky (Reference Madigosky1967), respectively. The correlated Enskog model considers the effect of gas attraction on shear viscosity and thermal conductivity using the approach of Chung et al. (Reference Chung, Ajlan, Lee and Starling1988), and on the bulk viscosity using the approach of Hoover et al. (Reference Hoover, Ladd, Hickman and Holian1980b). The softened Enskog uses a state-dependent molecule diameter (2.7) in the Enskog prediction of transport coefficients (2.32). The Enskog (corrected $\mu _B$) uses the corrected bulk viscosity in the Enskog prediction of shear viscosity (2.33).

Similar to shear viscosity and thermal conductivity, molecular attraction is important for bulk viscosity. However, the bulk viscosity is predicted more accurately by the Enskog formula (2.32) with a state-dependent diameter (2.7), i.e. the softened Enskog prediction, as shown in figures 2(e) and 2(f). Consequently, the Enskog prediction of shear viscosity and thermal conductivity is in better agreement with the experimental data at large densities if we use this corrected bulk viscosity in (2.33), replacing the original hard-sphere bulk viscosity $\mu _B^{hs}$ (see figures 2a,b,c,d). So the overestimation of bulk viscosity from the Enskog theory leads to the overestimation of shear viscosity and thermal conductivity (see figures 2a,b,c,d) at large densities.

To evaluate the effect of viscosity models on gas dynamics, a Poiseuille-type flow is investigated with the bottom and top wall temperatures $T_c=173$ K and $T_h=373$ K, respectively, the averaged density $\rho _{avg}=150\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width $H=15$ nm, and the external force $F_{ext}=0.0003\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$. The profiles of gas density, viscosity and velocity distributions are shown in figure 3. Although the density distribution is barely changed, the molecular attraction affects the velocity and viscosity profiles obviously.

Figure 3. The effect of viscosity models on (a) density and viscosity, and (b) velocity profiles, where ‘viscosity without attraction’ refers to the hard-sphere model (2.8), and ‘viscosity with attraction’ refers to the Chung model (2.25).

The present kinetic model (2.17) will then be evaluated by comparison with the simulation results of MD, the Shakhov–Enskog model (Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020), and the NS equations. For incompressible, steady-state and laminar flows, the NS equations reduce to

(3.3a)\begin{gather} \mu\,\frac{\partial^2 u}{\partial y^2}+F_{ext}n=0, \end{gather}

(3.3b)\begin{gather}\kappa\,\frac{\partial^2~T}{\partial y^2}+\mu \left(\frac{\partial u}{\partial y}\right)^2=0, \end{gather}

with the second-order boundary condition for velocity slip and the first-order boundary condition for temperature jump, namely

(3.4a)\begin{gather} u_s=\left.\pm A_1\lambda\,\frac{\partial u}{\partial y}\right|_{y=0}-\left.A_2\lambda^2\,\frac{\partial^2 u}{\partial y^2}\right|_{y=0}, \end{gather}

(3.4b)\begin{gather}T_j=\left.\beta\,\frac{2\gamma}{\gamma+1}\,\frac{\lambda}{Pr}\,\frac{\partial T}{\partial y}\right|_{y=0}, \end{gather}

where the slip coefficients $A_1=1.0$ and $A_2=0.5$ are chosen (Chapman & Cowling Reference Chapman and Cowling1990), and $\beta =(2-\sigma _T)/\sigma _T$ with the chosen thermal accommodation coefficient $\sigma _T=1.0$.

3.2. Poiseuille flows

In Poiseuille flows, an external body force acts on all the fluid molecules in the $x$ direction with the wall temperature $T_w=273$ K. By solving (3.3) and (3.4), the velocity and temperature distribution across the channel can be obtained as

(3.5a)\begin{gather} u(y)={-}\frac{F_{ext}n}{2\mu}\,[y^2-yH-H^2(A_1\,{Kn} +2A_2\,{Kn}^2)], \end{gather}

(3.5b)\begin{gather}T(y)={-}\frac{(F_{ext}n)^2}{24\mu \kappa}\,(2y^4-4Hy^3+3H^2y^2-H^3y-L_TH^4)+T_w, \end{gather}

where $Kn$ is the Knudsen number defined as

(3.6)\begin{equation} {Kn}=\frac{1}{\sqrt{2}\,n{\rm \pi} \sigma^2\chi H}, \end{equation}

and $L_T$ is the thermal jump length expressed as

(3.7)\begin{equation} L_T=\beta\,\frac{2\gamma}{\gamma+1}\,\frac{{Kn}}{Pr}. \end{equation}

Figure 4 shows the density and velocity profiles of the Poiseuille flows under a small external body force at different densities, i.e. different degrees of non-equilibrium (rarefaction) effect. As shown, the results of our kinetic model are in good agreement with the MD data for a broad range of densities (the reduced number density $\eta$ ranges from 0.00031 to 0.14). By contrast, the Shakhov–Enskog model (Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020), which neglects the gas molecular attraction, overestimates the density near the wall and underestimates the overall velocity profiles, particularly at large densities. The reduced fluid density at the wall when the molecular attraction is considered is due to the ‘pull’ of the molecules in the bulk. The NS prediction, on the other hand, is better at large densities where the non-equilibrium effect is not significant.

Figure 4. Density and velocity profiles: (a,b) for $\rho _{avg}=10\ \textrm {kg}\ \textrm {m}^{-3}$ ($\eta = 0.00031$, $Kn = 2.56$); (c,d) for $\rho _{avg}=150\ \textrm {kg}\ \textrm {m}^{-3}$ ($\eta = 0.047$, $Kn = 0.15$); and (e,f) for $\rho _{avg}=450\ \textrm {kg}\ \textrm {m}^{-3}$ ($\eta = 0.14$, $Kn =0.039$). The external force $F_{ext} = 0.001\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$ is small, so the viscous dissipation is negligible, the channel width is $H=5$ nm, and the wall temperature is $T_w = 273$ K.

For high-speed flows, the viscous dissipation plays an important role, which is investigated in figure 5 with the average density $\rho _{avg}=350\ \textrm {kg}\ \textrm {m}^{-3}$, channel width $H=5$ nm, and wall temperature $T_w=273$ K. Two large external forces are considered, namely $F_{ext}=0.01$ and $0.02\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$, respectively. Again, the density oscillation, parabolic velocity and quartic temperature profiles are well captured by the current kinetic model, while the Shakhov–Enskog model and the NS equation show large errors. The discrepancy in the results between the current kinetic model and the Shakhov–Enskog model suggests the important role of the long-range molecular attraction in gas dynamics, leading to reduced density near the wall, and enhanced velocity slip and temperature jump.

Figure 5. The density, velocity and temperature profiles at different external forces: (a–c) for $F_{ext} = 0.01\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$; and (d–f) for $F_{ext} = 0.02\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$. The average density is $\rho _{avg}= 350\ \textrm {kg}\ \textrm {m}^{-3}$ ($\eta =0.11$), the channel width is $H=5$ nm, and the wall temperature is $T_w= 273$ K. The resulting Knudsen number is $Kn=0.055$.

The effect of temperature on density and velocity profiles is shown in figure 6, with the average density $\rho _{avg}=350\ \textrm {kg}\ \textrm {m}^{-3}$, channel width $H=5$ nm, and external force $F_{ext}=0.001\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$. It is very clear that the density and velocity profiles predicted by our kinetic model agree with the MD data, while the NS equation fails to predict density variation, and the Shakhov–Enskog model overpredicts the density near the wall. The main difference between our model and the Shakhov–Enskog model is that we include the gas molecular attraction, which can pull the gas molecules to the bulk. As a result, our prediction of the density at the wall is significantly smaller than the Shakhov–Enskog model for hard-sphere molecules, which ignores the molecular attraction. As shown in figure 6(e), even at high temperatures, the gas density of van der Waals fluids is still significantly affected by the long-range molecular attraction, i.e. the gas molecular attraction is not negligible even at high temperatures, which has largely been ignored in the previous studies.

Figure 6. The density and velocity profiles at different temperatures: (a,b) for $T=253$ K; (c,d) for $T=313$ K; and (e,f) for $T=373$ K. The average density is $\rho _{avg}= 350\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width is $H=5$ nm, and the external force is $F_{ext} = 0.001\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$.

The velocity decreases with the temperature, as shown in figures 6(b), 6(d) and 6(f), which is caused by the larger near-wall density and viscosity at high temperatures. The larger near-wall density means more efficient momentum transfer between the fluid and the wall, leading to smaller velocity slips at the solid surfaces, while the larger viscosity means more flow resistance for bulk fluid in the channel. This can be seen more clearly by normalising the slip velocity $u_s$ by $F_{ext}H/(mu_m)$, namely

(3.8a,b)\begin{equation} \hat{u}_s=\frac{u_smu_m}{F_{ext}H},\quad u_m=\sqrt{\frac{2k_BT}{m}}, \end{equation}

where $u_m$ is the most probable velocity. The variation of the normalised slip velocity with temperature is shown in figure 7. The normalised slip velocity of hard-sphere gases predicted by the Shakhov–Enskog model is nearly constant as the temperature changes, while the present kinetic model and MD simulation predict a decreasing slip velocity with temperature. For the hard-sphere gases, the density distribution is not affected by the temperature and the velocity distribution $u(y)\propto 1/\mu _s^{hs}$. As shown by (2.8), the shear viscosity of hard-sphere gases satisfies $\mu _s^{hs} \propto \sqrt {T}$, so the normalised velocity is temperature-independent as $u_m \propto \sqrt {T}$. However, the relationship between viscosity and temperature $\mu _s \propto \sqrt {T}$ no longer holds for the van der Waals fluids as the collision integral $\varOmega ^{(2,2)}$, which is temperature-dependent, comes into play. Furthermore, a new dimensionless number, namely the reduced temperature $\hat {T}=k_BT/\epsilon$, is introduced to signify the competition between gas molecular attraction and kinetic energy for the van der Waals fluids. As the temperature increases, the gas molecules gain more kinetic energy to overcome the attractive forces holding them to the bulk. This results in greater accumulation of gas molecules near the walls, leading to the increased momentum transfer between the solid and the gas, thus reducing the slip velocity.

Figure 7. The variation of the normalised slip velocity with temperature.

As the viscous dissipation is non-negligible for fluids under high shear rates, we investigate its effect on the normalised mass flow rate $Q_n$, which is defined as

(3.9)\begin{equation} Q_n= \frac{\displaystyle\int\nolimits_{0}^{H} n(y)\,u(y)\,\mathrm{d}y}{n_{avg}F_{ext}H^2/(mu_m)}. \end{equation}

As shown in figure 8, increased viscous dissipation reduces the mass flow rate in all the flow regimes. This is because the viscous dissipation leads to a smaller slip velocity and a larger flow resistance, as shown by the results in figure 5. The effect of viscous dissipation on mass flow rate is similar to that of confinement, which is also included in figure 8 for comparison. However, the confinement reduces the mass flow rate more significantly for small-${Kn}$ flows, resulting in the disappearance of the Knudsen minimum. On the contrary, the viscous dissipation flattens the variation curve ($Q_n \sim Kn$), but no Knudsen minimum disappearance is observed.

Figure 8. The variation of normalised mass flow rate with the Knudsen number at different external forces and confinements.

3.3. Couette flows

In the Couette flows, the top and bottom walls move at speed $u_w$ in opposite directions, as shown in figure 1(b). No external force is exerted on fluid molecules, so $F_{ext}=0$ and the wall temperature is set to be $T_w=273$ K. The velocity and temperature profiles can also be obtained by solving (3.3) and (3.4), which are written as

(3.10a)\begin{gather} u(y)=\frac{2u_w}{H}\,y-u_w, \end{gather}

(3.10b)\begin{gather}T(y)={-}\frac{2\mu u_w^2}{\kappa H^2}\left(y^2-Hy-L_T H^2\right)+T_w. \end{gather}

Figure 9 shows the density, velocity and temperature profiles of the Couette flows at different shear rates, with the average density $\rho _{avg}=350\ \textrm {kg}\ \textrm {m}^{-3}$, and the channel width $H=5$ nm. The present kinetic model captures the density oscillation, linear velocity distribution and parabolic temperature distribution, which are in good agreement with the MD data. Similar to the Poiseuille flow, the Shakhov–Enskog model, which neglects the long-range attraction between gas molecules, overestimates the density near the wall, and underestimates both the velocity slip and the temperature jump. As the shear rate increases, gas molecules are more likely to accumulate near the wall, as the long-range molecular attraction may not be sufficient to pull the gas molecules to the bulk, also shown in figure 10. In contrast to the density peak near the wall, the viscosity is the lowest in this region; see figure 10(d). This is because the bulk gas has higher temperatures due to viscous heating. Figure 10(b) shows that stronger viscous dissipation causes a reduction in velocity slip as a combined consequence of a higher density peak and a greater viscosity near the wall.

Figure 9. The density, velocity and temperature profiles: (a–c) $u_w=200\ \textrm {m}\ \textrm {s}^{-1}$; and (d–f) $u_w=400\ \textrm {m}\ \textrm {s}^{-1}$. The average density is $\rho _{avg}= 350\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width is $H=5$ nm, and the wall temperature is $T_w= 273$ K.

Figure 10. Distribution of (a) density, (b) velocity, (c) temperature, and (d) viscosity across the channel at different wall velocities. The average density is $\rho _{avg}= 350\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width is $H=5$ nm, and the wall temperature is $T_w= 273$ K.

The viscous dissipation effect on Couette flows under tighter confinement is also investigated, where the channel width shrinks from 5 nm to 2 nm, as shown in figure 11. For such a case, both the non-equilibrium and confinement effects become stronger. The results from the Shakhov–Enskog model and the NS equations exhibit larger discrepancies compared to the MD simulation results, while our kinetic model can still capture accurately the density, velocity and temperature profiles.

Figure 11. The density, velocity and temperature profiles: (a–c) $u_w=200\ \textrm {m}\ \textrm {s}^{-1}$; and (d–f) $u_w=400\ \textrm {m}\ \textrm {s}^{-1}$. The average density is $\rho _{avg} =350\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width is $H=2$ nm, and the wall temperature is $T_w =273$ K. The resulting Knudsen number is $Kn=0.14$.

3.4. Couette–Fourier flows

The Couette–Fourier flow differs from the Couette flow only in the different wall temperatures, with the top wall temperature at $T_h=373$ K and the bottom wall temperature at $T_c=273$ K. By solving (3.3) and (3.4), the velocity and temperature can be obtained as

(3.11a)\begin{gather} u(y)=\frac{2u_w}{H}\,y-u_w, \end{gather}

(3.11b)\begin{gather}\frac{T(y)-T_c}{T_h-T_c}={-}2\,{Br} \left[\left(\frac{y}{H}\right)^2-(1-2L_T)\,\frac{y}{H}-L_T(1-2L_T) \right]+\frac{y}{H}+L_T, \end{gather}

where $Br =\mu u_w^2/[\kappa (T_h-T_c)]$ is the Brinkman number measuring the competition between viscous heating and thermal conduction.

The density, velocity and temperature profiles of the Couette–Fourier flows at two different wall velocities are shown in figure 12, with the channel width $H=5$ nm. At a small wall moving velocity ($u_w=50\ \textrm {m}\ \textrm {s}^{-1}$), the viscous dissipation is negligible, so the density and temperature distributions recover those of the Fourier flows, while the velocity distribution is similar to the Couette flows. When the wall velocity increases to $u_w=300\ \textrm {m}\ \textrm {s}^{-1}$, the temperature profile becomes a combination of the linear and parabola distributions resulting from the Fourier and Couette flows, respectively. Again, the present kinetic model predicts these profiles accurately, while the results of the Shakhov–Enskog model deviate significantly from the MD data, particularly for the density and temperature profiles.

Figure 12. The density, velocity and temperature profiles of the Couette–Fourier flows at different wall velocities: (a–c) $u_w=50\ \textrm {m}\ \textrm {s}^{-1}$; and (d–f) $u_w=300\ \textrm {m}\ \textrm {s}^{-1}$. The average density is $\rho _{avg}=350\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width is $H=5$ nm, and the top and bottom wall temperatures are $T_h=373$ K and $T_c=273$ K, respectively.

With increased viscous dissipation at large wall velocities, the heat generated in the gases leads to higher gas temperatures, as shown in figure 13(a). The viscous dissipation increases the heat transfer rate between the gas and the cold (bottom) wall, while it limits the heat transfer rate between the gas and the hot (top) wall, which can be seen clearly from the heat flux variation in figure 13(b). When the wall velocity is sufficiently large, the hot wall can also be heated by the gases due to the large amount of heat generated by viscous heating. Thus our kinetic model may provide a design simulation tool to develop next-generation technologies such as nanoscale evaporative cooling.

Figure 13. (a) Temperature and (b) heat flux distributions of the Couette–Fourier flows at different wall velocities, where the symbols denote the MD data. The averaged density is $\rho _{avg}=350\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width is $H=5$ nm, and the top and bottom wall temperatures are $T_h=373$ K and $T_c=273$ K, respectively.

3.5. Model solution in the dilute and continuum limits

As shown in figures 14(a) and 14(b), the results of the present kinetic model for the van der Waals fluids are in good agreement with the Shakhov–Boltzmann model for dilute gases and the MD simulation when the real gas effects (namely the volume exclusion and the long-range molecular attraction) and the confinement effect are negligible. This is because the density terms $\mathcal {I}^{(1)}$ (2.15) and $\mathcal {I}^{(2)}$ (2.16) and the mean-field force (2.6) become negligible, and the kinetic model (2.17) reduces to the Shakhov model for the hard-sphere molecules in the dilute limit. Meanwhile, the results of our kinetic model, NS equations and MD simulations are very close in the continuum limit where the non-equilibrium and confinement effects are sufficiently small, as shown in figures 14(c) and 14(d). This is also expected because the NS equations are recovered from the kinetic model (2.17) in the small-${Kn}$ limit, as shown in Appendix A. Therefore, the present kinetic model, which is an extension of the EV model for the hard-sphere molecules to consider real gas effects, is capable of simulating non-equilibrium flows of the confined van der Waals fluids.

Figure 14. The present kinetic model agrees with (a,b) the Boltzmann model in the dilute limit, and (c,d) the NS equations in the continuum limit. The average density $\rho _{avg}= 2\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width $H=100$ nm and the external force $F_{ext}=0.00003\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$ in (a,b) correspond to $\eta =0.00062$ and $Kn = 0.64$; the average density $\rho _{avg}= 450\ \textrm {kg}\ \textrm {m}^{-3}$, the channel width $H=50$ nm and the external force $F_{ext}=0.00005\ \textrm {kcal}\ \textrm {mol}^{-1}\ \textrm {\AA }^{-1}$ in (c,d) correspond to $\eta =0.14$ and $Kn = 0.0039$.