2.1. Fluid–fluid interactions

According to the kinetic theory for dense fluids and mean-field theories, the fluid–fluid interactions $Q_{\textit{ff}}$ can be decomposed into short-range repulsive interactions $Q_{\textit{ff}}^{\textit{sg}}$ for $|\boldsymbol{r}'|=|\boldsymbol{r}_2-\boldsymbol{r}_1|\lt \sigma _f$ and long-range attractive interactions $Q_{\textit{ff}}^{lg}$ for $|\boldsymbol{r}'| \geqslant \sigma _f$ . This decomposition (Shan et al. Reference Shan, Su, Gibelli and Zhang2023b ) is expressed as

(2.2a) \begin{align} Q_{\textit{ff}} & = Q_{\textit{ff}}^{\textit{sg}} + Q_{\textit{ff}}^{lg}, \end{align}

(2.2b) \begin{align} Q_{\textit{ff}}^{\textit{sg}} & = \chi (\bar {n}) Q_{\textit{ff}}^0 + Q_{\textit{ff}}^1, \end{align}

(2.2c) \begin{align} Q_{\textit{ff}}^{lg} & = m^{-1}\nabla _{\boldsymbol{\xi }}f\cdot \nabla \phi _{\textit{ff}}^{lg}, \end{align}

where $m$ is the mass of a fluid molecule, $\chi$ is the radial distribution function (Carnahan & Starling Reference Carnahan and Starling1969), which quantifies the increase in collision frequency between fluid molecules due to their finite size. It can be calculated from the local average density $\bar {n}=\int w(\boldsymbol{r}')n(\boldsymbol{r}+\boldsymbol{r}')\mathrm{d}\boldsymbol{r}'$ (Tarazona Reference Tarazona1985; Tarazona, Marconi & Evans Reference Tarazona, Marconi and Evans1987) as

(2.3) \begin{equation} \chi (\bar {n})=\frac {1}{2}\frac {2-\eta (\bar {n})}{[1-\eta (\bar {n})]^3}, \quad \eta (\bar {n}) = \frac {\bar {n}V_0}{4},\quad V_0 = \frac {2\pi \sigma _f^3}{3}, \end{equation}

with $w(\boldsymbol{r}')$ being a weight function, $\eta$ the reduced number density and $V_0$ the covolume of fluid molecules. This approach has successfully reproduced the fluid structure under strong confinements (Vanderlick et al. Reference Vanderlick, Scriven and Davis1989; Guo et al. Reference Guo, Zhao and Shi2005; Frezzotti & Gibelli Reference Frezzotti and Gibelli2008; Bhadauria et al. Reference Bhadauria, Sanghi and Aluru2015; Bhadauria & Aluru Reference Bhadauria and Aluru2017). As $\chi$ varies with the local average density, it is space-dependent for inhomogeneous fluids. Besides, $\phi_{ff}^{lg} = \int_{|\boldsymbol{r}'|\gt \sigma _f} n(\boldsymbol{r}+\boldsymbol{r}')\phi _{att}(|\boldsymbol{r}'|)\mathrm{d}\boldsymbol{r}'$ is the mean-field attractive potential with $\phi _{att}(|\boldsymbol{r}'|)=-4\epsilon _{\textit{ff}}(\sigma _{f}/|\boldsymbol{r}'|)^6$ being the attractive part of the Lennard–Jones (LJ) potential and $\epsilon _{\textit{ff}}$ the fluid–fluid interaction strength.

The short-range repulsion $Q_{\textit{ff}}^{\textit{sg}}$ consists of a collision term ( $Q_{\textit{ff}}^0$ ) and a finite-size correction ( $Q_{\textit{ff}}^1$ ). The collision operator $Q_{\textit{ff}}^{0}$ can be simplified using the classic relaxation approach. Here, we adopt the Shakhov model (Shakhov Reference Shakhov1968), i.e.

(2.4a) \begin{align} Q_{\textit{ff}}^{0} & = -\tau _f^{-1}\left (f-f^{Sh}\right ), \end{align}

(2.4b) \begin{align} f^{Sh} & = f^{eq}\left [1+(1-Pr )\frac {\boldsymbol{c}\cdot \boldsymbol{Q}_k}{5p_0RT}\left (\frac {|\boldsymbol{c}|^2}{RT}-5\right )\right ], \end{align}

(2.4c) \begin{align} f^{eq} & = n(2\pi RT)^{-\frac {3}{2}}\exp \left (-\frac {|\boldsymbol{c}|^2}{2RT}\right ), \end{align}

where $\tau _f$ is the relaxation time, $f^{Sh}$ is the reference VDF, $p_0=nk_B T$ is the ideal gas pressure with $k_B$ being the Boltzmann constant and $T$ the temperature, $Pr$ is the Prandtl number, $\boldsymbol{Q}_k$ is the kinetic heat flux, $R=k_B/m$ is the specific gas constant, $f^{eq}$ is the local Maxwellian distribution function and $\boldsymbol{c}=\boldsymbol{\xi }-\boldsymbol{u}$ is the peculiar velocity with $\boldsymbol u$ being the flow velocity. This model improves upon the Bhatnagar–Gross–Krook collision operator (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954) to recover the correct Prandtl number. This technique has been proved as an effective procedure to capture the complex dynamics of dense fluids (Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020; Su et al. Reference Su, Gibelli, Li, Borg and Zhang2023; Li et al. Reference Li, Su, Shan, Li, Gibelli and Zhang2024a ).

The first-order term $Q_{\textit{ff}}^{1}$ accounts for the finite size effect (Chapman & Cowling Reference Chapman and Cowling1990), which prevents fluid molecules from overlapping. It approximates the repulsive forces among inhomogeneous fluid molecules (Shan et al. Reference Shan, Ju, Su, Guo and Zhang2023a ) and can be expressed as

(2.5) \begin{align} Q_{\textit{ff}}^{1}&=-\bar {n}V_0\chi (\bar {n}) f^{eq}\left \{\frac {2}{5}\left [2\boldsymbol{\mathcal{C}}\boldsymbol{\mathcal{C}}:\nabla \boldsymbol{u} +\left (|\boldsymbol{\mathcal{C}}|^2-\frac {5}{2}\right )\nabla \cdot \boldsymbol{u} \right ]\right . \notag \\ &\quad +\left .\boldsymbol{c}\cdot \left [ \frac {2\boldsymbol{A}}{\bar {n}}+\frac {\boldsymbol{B}}{\chi (\bar {n})} +\frac {3}{5}\left (|\boldsymbol{\mathcal{C}}|^2-\frac {5}{6}\right ) \nabla \ln T\right ]\right \}, \end{align}

where $\boldsymbol{\mathcal{C}}=\boldsymbol{c}/u_m$ is the non-dimensional peculiar velocity with $u_m=\sqrt {2RT}$ being the most probable speed, and $\boldsymbol{A}$ and $\boldsymbol{B}$ are vector functions (Guo et al. Reference Guo, Zhao and Shi2005), expressed as

(2.6a) \begin{align} \boldsymbol A(\boldsymbol r) & = \frac {120}{\pi \sigma _f^5}\int _{|\boldsymbol r'|\lt \sigma _f/2} \boldsymbol r' \bar {n}(\boldsymbol r+\boldsymbol r')\:\mathrm{d}\boldsymbol r', \end{align}

(2.6b) \begin{align} \boldsymbol B(\boldsymbol r) & = \frac {120}{\pi \sigma _f^5}\int _{|\boldsymbol r'|\lt \sigma _f/2} \boldsymbol r' \chi [\bar {n}(\boldsymbol r+\boldsymbol r')]\:\mathrm{d}\boldsymbol r'. \end{align}

2.2. Fluid–surface interactions

Similar to fluid–fluid interactions, we decompose fluid–surface interactions into short-range ( $Q_{\textit{fs}}^{\textit{sg}}$ for $h\lt H_{ex}$ ) and long-range ( $Q_{\textit{fs}}^{lg}$ for $h \geqslant H_{ex}$ ) contributions:

(2.7a) \begin{align} Q_{\textit{fs}} & = Q_{\textit{fs}}^{\textit{sg}}+ Q_{\textit{fs}}^{lg}, \end{align}

(2.7b) \begin{align} Q_{\textit{fs}}^{\textit{sg}} & = Q_{\textit{fs}}^0 +Q_{\textit{fs}}^{ex}, \end{align}

(2.7c) \begin{align} Q_{\textit{fs}}^{lg} & = m^{-1}\nabla _{\boldsymbol{\xi }}f\cdot \nabla \phi _{\textit{fs}}^{lg}. \end{align}

Here, we utilise the 10-4-3 potential (Steele Reference Steele1973) to characterise long-range fluid–surface interactions, i.e.

(2.8) \begin{equation} \phi _{\textit{fs}}^{lg}(y) =2n_D\pi \epsilon _{\textit{fs}}\left [\frac {2}{5}\left (\frac {\sigma _{\textit{fs}}}{y}\right )^{10} -\left (\frac {\sigma _{\textit{fs}}}{y}\right )^{4}-\frac {\sigma _{\textit{fs}}^4}{3\Delta (y+0.61\Delta )^3} \right ], \end{equation}

where $y$ is the position in terms of the physical wall, $n_D$ is the non-dimensional solid number density, $\Delta$ is the average distance between two adjacent layers in solids, see figure 1(a), and $\sigma _{\textit{fs}}$ and $\epsilon _{\textit{fs}}$ are the length and energy parameters for fluid–surface interactions, which can be calculated using the Lorentz–Berthelot combination rules (Lorentz Reference Lorentz1881; Berthelot Reference Berthelot1898), i.e. $\epsilon _{\textit{fs}}=\sqrt {\epsilon _{\textit{ff}}\epsilon _{ss}},\,\sigma _{\textit{fs}}=(\sigma _{f}+\sigma _{s})/2$ with $\epsilon _{ss}$ representing the solid–solid interaction strength.

The short-range interaction ( $Q_{\textit{fs}}^{\textit{sg}}$ ) includes both information (momentum and energy) exchanges between the fluid and the solid ( $Q_{\textit{fs}}^{0}$ ) and the fluid–solid volume exclusion effect ( $Q_{\textit{fs}}^{ex}$ ). Specifically, the term $Q_{\textit{fs}}^{ex}$ accounts for the finite size of fluid and solid molecules, as well as the lattice structure of the solid substrate, which creates a depletion region near the solid surface (Joseph & Aluru Reference Joseph and Aluru2008; Hadjiconstantinou Reference Hadjiconstantinou2024). This region defines the hydrodynamic wall location, typically on the fluid side of the depletion area (Bocquet & Barrat Reference Bocquet and Barrat1994; Bhadauria & Aluru Reference Bhadauria and Aluru2017). The exact location depends on multiple factors, such as the fluid–solid interaction strength and surface lattice structure (Bhadauria et al. Reference Bhadauria, Sanghi and Aluru2015; Hadjiconstantinou Reference Hadjiconstantinou2024) and plays a crucial role in modelling NCFs (Kobryn & Kovalenko Reference Kobryn and Kovalenko2008; De La Torre et al. Reference de la Torre, Duque-Zumajo, Camargo and Español2019). We will subsequently introduce how to model $Q_{\textit{fs}}^{ex}$ and $Q_{\textit{fs}}^{0}$ , respectively.

2.2.1. Fluid–solid volume exclusion

As there are no fluid molecules in the fluid–solid volume exclusion region ( $h\lt H_{ex}\!$ ), the modelling of $Q_{\textit{fs}}^{ex}$ is essentially equivalent to the determination of the hydrodynamic wall. Following the approach used in the kinetic theory for dense fluids (Chapman & Cowling Reference Chapman and Cowling1990), we begin by considering the molecules of the hard-sphere fluid and solid with diameters of $\sigma _f$ and $\sigma _s$ , respectively. The solid molecules are bound by spring-like forces, as illustrated in figure 2(a). The distance between adjacent lattice sites is $2a$ , where $2a$ can be assumed no larger than $\sigma _{\textit{fs}}$ , a condition satisfied in most physical scenarios. The amplitude of a solid molecule’s vibration is $d$ , as illustrated in figures 2(a) and 2(b). Using the Debye model (Born & Huang Reference Born and Huang1996), the root-mean-square amplitude of the vibration can be estimated as

(2.9) \begin{equation} d \simeq \sqrt {\frac {3\hbar ^2}{Mk_B\Theta _D} \left (\frac {T_w}{\Theta _D} + \frac {1}{4}\right )}, \end{equation}

where $\hbar$ is the reduced Planck constant, $M$ is the mass of the vibrating molecules, $\Theta _D$ is the Debye temperature and $T_w$ is the solid temperature. Take gold molecules to represent the solid substrate as an example. Within a temperature range of 100–600 K, the calculated vibration amplitude of solid molecules is approximately $0.06{-}0.12$ Å, significantly smaller than the lattice spacing ( $4.07$ Å) and molecular diameter ( $3.29$ Å) of the solid substrate. Therefore, it can be assumed that solid molecules remain effectively at their equilibrium positions due to minimal vibration amplitudes.

Figure 2.

(a) Lattice vibration around the equilibrium position; (b) the potential a solid molecule feels during vibration; (c) the position of the hydrodynamic wall as an average of two extreme contacts.

As shown in figure 2(c), the physical and hydrodynamic walls do not coincide in NCFs. The physical wall is defined at the mass centres of the innermost solid molecules that contact fluids. Given the negligible vibration amplitude, the physical wall can be approximated as the equilibrium position, as illustrated in figure 2(a). To determine the hydrodynamic wall location, we consider the interaction between fluid molecules and a solid substrate, taking into account the non-negligible molecular sizes and the lattice structure of the solid. The extreme distances between the physical wall and the incident fluid molecules are denoted as $H_1$ and $H_2$ , respectively, as shown in figure 2(c), which can be calculated as

(2.10a) \begin{align} H_1 & = \frac {\sigma _f+\sigma _s}{2}=\sigma _{\textit{fs}}, \end{align}

(2.10b) \begin{align} H_2 & = \sqrt {\sigma _{\textit{fs}}^2-a^2}. \end{align}

Following the Lorentz combination rule (Lorentz Reference Lorentz1881) commonly used in MD simulations, the hydrodynamic wall location is similarly determined as the arithmetic average of these two distances, namely $H_{ex}=(H_1+H_2)/2$ .

To test our molecular kinetic model with the above hydrodynamic wall, we calculate the solvation force for various confinement space $H$ . The solvation force, defined as the average normal force per unit area exerted by fluid molecules on the surface (with dimensions equivalent to pressure), exhibits oscillatory behaviour as the confinement space changes. It can be estimated from our model as

(2.11) \begin{equation} F_s \approx-\int _{H_{ex}}^{H-H_{ex}} n(y) \frac {\textrm {d} \phi _{\textit{fs}}^{lg} (y) } {\textrm {d} y} {\textrm {d} y}. \end{equation}

where $F_s$ is the solvation force related to the disjoining pressure $p_d$ and the bulk pressure $p_b$ through the relationship $F_s=p_d +p_b$ .

As illustrated in figure 3, the solvation force oscillates with increasing confinement and converges to the bulk pressure $p_b$ as the confinement space $H$ becomes sufficiently large. Our molecular kinetic model reproduces this oscillatory behaviour, as observed in the Monte Carlo simulations (Snook & van Megen Reference Snook and van Megen1980), confirming the appropriateness of the present molecular kinetic model and the way the hydrodynamic wall is determined. Therefore, the hydrodynamic wall location can be approximately chosen at $h=H_{ex}$ , which means that the term $Q_{\textit{fs}}^{ex}$ can be neglected in the hydrodynamic bulk flow region ( $h\gt H_{ex}$ ) when modelling fluid–surface interactions.

Figure 3.

Oscillation of the solvation force compared with the Monte Carlo simulations (Snook & van Megen Reference Snook and van Megen1980). The following parameters are used: average number density $n=0.5925\sigma _f^{-3}$ ; temperature $T=1.2\epsilon _{\textit{ff}}/k_B$ ; $\epsilon _{ss}=\epsilon _{\textit{ff}}$ .

2.2.2. Fluid–surface collisions

To model the collision term $Q_{\textit{fs}}^0$ , which describes the momentum and energy exchanges between fluids and solids during molecular collisions, we assume that solid molecules remain in Maxwellian equilibrium (Frezzotti & Gibelli Reference Frezzotti and Gibelli2008; Aoki, Giovangigli & Kosuge Reference Aoki, Giovangigli and Kosuge2022), unaffected by fluid–solid collisions due to their strong rigidity, as depicted in figure 1(c). Given that the rates of momentum and energy transfer between fluids and solid surfaces differ (Cercignani Reference Cercignani1988; Woods Reference Woods1993; Struchtrup Reference Struchtrup2013), we introduce a total energy distribution function $g=0.5|\boldsymbol{\xi }|^2f$ . Thus, $Q_{\textit{fs}}^0$ for momentum ( $Q_{\textit{fs}}^{m}$ ) and energy ( $Q_{\textit{fs}}^{e}$ ) can be modelled separately, each by a relaxation process (Beenakker, Borman & Krylov Reference Beenakker, Borman and Krylov1994; Aoki et al. Reference Aoki, Giovangigli and Kosuge2022),

(2.12a) \begin{align} Q_{\textit{fs}}^{m} & = -\tau _{m}^{-1} \left (f-f_w^{eq}\right ), \end{align}

(2.12b) \begin{align} Q_{\textit{fs}}^{e} & = -\tau _{e}^{-1} \left (g-g_w^{eq}\right ), \end{align}

where $\tau _m$ and $\tau _e$ represent the relaxation times of momentum and energy, which quantify the exchange rates of momentum and energy between the fluid and the solid. These exchanges occur only when fluid molecules are in close contact with the solid surface, so both $\tau _m$ and $\tau _e$ take finite values at $h=H_{ex}$ . Away from the surface, where there is no momentum or energy exchange between the fluid and the solid, these relaxation times effectively tend to infinity. Given this physical interpretation, $\tau _m$ and $\tau _e$ only play roles at $h=H_{ex}$ , as the infinite relaxation at other places does not change the fluid dynamics. In this regard, the collision terms $Q_{\textit{fs}}^m$ and $Q_{\textit{fs}}^e$ act as effective boundary conditions at $h=H_{ex}$ and their influence is localised to the solid surface. Meanwhile, the solid Maxwellian can be written as

(2.13a) \begin{align} f_w^{eq} & = n_w (2\pi RT_w)^{-\frac {3}{2}}\exp \left (-\frac {|\boldsymbol{\xi } - \boldsymbol{u}_w|^2}{2RT_w}\right ), \end{align}

(2.13b) \begin{align} g_w^{eq} & = 0.5|\boldsymbol{\xi }|^2f_w^{eq}, \end{align}

where $n_w$ is the fluid density at the hydrodynamic wall, and $\boldsymbol{u}_w$ is the moving velocity of the surface. Since no fluid penetrates the solid, zero mass flux is achieved at $h=H_{ex}$ .

4.1. Numerical procedures

Both the proposed molecular kinetic model and MD simulations are employed to investigate NCFs, aiming to demonstrate the accuracy and predictive capability of our model. In the MD simulations, the lowest fluid density considered is 60 kg m $^{-3}$ , as lower densities introduce significant uncertainty and statistical noise, particularly under strong confinements, e.g. the flows in a 2 nm channel.

We simulate argon confined between two planar solid surfaces with a face-centred-cubic crystal structure. Argon is selected as the working fluid due to its simplicity and well-characterised intermolecular potential (Rahman Reference Rahman1964; Fan et al. Reference Fan, Wu and Wang2020b ). Following established methodologies (Bhadauria et al. Reference Bhadauria, Sanghi and Aluru2015; Bhadauria & Aluru Reference Bhadauria and Aluru2017), our model can potentially be extended to accommodate other fluids and electroosmotic flows. To simulate fluid flows, argon atoms are subjected to an external acceleration rate $\boldsymbol{a}_e$ .

For the thermal cases presented in figure 4, the fluid is not thermostatted, allowing for a direct examination of energy evolution in NCFs. In the investigation of the Knudsen minimum phenomenon (figures 5 and 6), an isothermal condition is maintained to better distinguish the effects of fluid–surface interactions and real fluid behaviours on fluid transport properties, where the fluid is governed by the Nosé–Hoover thermostat for stable temperature control, avoiding the need for small external forces that could introduce substantial statistical noise.

The solid atoms are tethered to their lattice sites via a linear spring model, allowing for vibration while maintaining sufficient rigidity to contain the argon atoms. The Nosé–Hoover thermostat is employed to control the temperature of the solid substrate, ensuring thermal stability throughout the simulations. The rough surface in figure 5(a) is achieved by randomly displacing the solid atoms in the $y$ direction, with the maximum displacement $\Delta y_{max}=0.44$ Å. The confinement space, i.e. the channel height $H$ , is defined as the distance between the innermost layer of solid atoms at their equilibrium positions, consistent with the physical wall shown in figure 1(b).

The interactions between argon and solid atoms are modelled using the 12-6 LJ potential with a cutoff distance of 1.3 nm, which is the standard practice in MD simulations. The relevant parameters are $\epsilon _{ss}=0.0390$ kcal mol $^{-1}$ , $\epsilon _{\textit{ff}}=0.2376$ kcal mol $^{-1}$ , $\sigma _{s}=3.2930$ Å and $\sigma _f=3.4039$ Å. For figures 6(c) and 6(d), wettability strengths, quantified by the parameter $\epsilon _D = \epsilon _{\textit{fs}}/\epsilon _{\textit{ff}}$ , are varied by multiplying the solid energy parameter $\epsilon _{ss}$ by a factor of 2.5 and 0.1, resulting in $\epsilon _D=0.641$ and $0.128$ , respectively. After achieving equilibrium, the simulations are run with a time step of 0.5 fs for a minimum duration of 5 ns to ensure accurate statistical sampling. The periodic boundary condition is imposed in all the directions. All MD simulations are conducted using the large-scale atomic/molecular massively parallel simulator (Thompson et al. Reference Thompson2022).

Figure 4.

The NCFs behaviours in comparison with MD simulations. Panels (a) and (b) are for $H=2.5$ nm, $T_w=233$ K, $|\boldsymbol{a}_e|=2.095\times 10^{12}$ m s $^{-2}$ and $\rho _{avg}=100$ kg m $^{-3}$ ;panels (c) and (d) are for $H=5.0$ nm, $T_w=233$ K, $|\boldsymbol{a}_e|=8.382\times 10^{11}$ m s $^{-2}$ and $\rho _{avg}=500$ kg m $^{-3}$ .

Figure 5.

(a) Mass flow rates in atomically smooth ( $\omega _m=0.30$ ) and rough ( $\omega _m=0.82$ ) channels. (b) Contributions of real fluid effects and fluid–surface interactions for flows in perfectly rough channels, i.e. $\omega _m=1.0$ .

The proposed molecular kinetic model is solved using the discrete unified gas kinetic scheme (Guo et al. Reference Guo, Xu and Wang2013, Reference Guo, Wang and Xu2015; Zhang et al. Reference Zhang, Zhu, Wang and Guo2018), which can simulate strongly inhomogeneous fluid flows (Shan et al. Reference Shan, Wang, Zhang and Guo2020). All the length and energy parameters used in the kinetic model, including $\sigma _f$ , $\sigma _s$ , $\epsilon _{\textit{ff}}$ and $\epsilon _{ss}$ , are consistent with those employed in MD simulations for meaningful comparisons.

The accommodation coefficients in our molecular kinetic model are flow system dependent, meaning that once the types of solid and fluid are specified, these coefficients remain constant for a broad range of density, velocity and confinement, if the temperature variation is small. To determine their values, we perform MD simulations and extract velocity and temperature profiles. We then fit these profiles using our molecular kinetic model to obtain the accommodation coefficients. Once determined, these values can be applied to other configurations involving the same fluid–solid system.

Figure 6.

Temperature and wettability effects on the mass flow rate. In (a) and (b), the fluid–surface accommodation coefficients are $\omega _m=0.53$ , $0.30$ and $0.23$ for $T=183$ K, $283$ K and $383$ K, respectively. In (c) and (d), the system temperature $T=283$ K, and the fluid–surface accommodation coefficients are $\omega _m=0.49$ , $0.30$ and $0.18$ for $\epsilon _D=0.641$ , 0.405 and 0.128, respectively.

4.2. Features of NCFs

For force-driven channel flows shown by figure 1(f), figure 4(a,b) presents the density, velocity and temperature profiles for $H=2.5$ nm, $T_w=233$ K, external acceleration $|\boldsymbol{a}_e|=2.095\times 10^{12}$ m s $^{-2}$ in the flow direction and an average density $\rho _{avg}=100$ kg m $^{-3}$ . The fluid–surface accommodation coefficients, $\omega _m=0.32$ and $\omega _e=0.10$ , are determined through independent MD simulations with the same LJ parameters. The results of the present molecular kinetic model align well with the MD simulation data, as shown in figure 4. A pronounced density peak, indicating fluid adsorption due to fluid–surface interactions, is observed near the solid surface and rapidly diminishes with distance. This layering structure of NCFs, though frequently reported in MD simulations, has been largely overlooked by previous kinetic models (Li, Xu & Li Reference Li, Xu and Li2010; Barisik & Beskok Reference Barisik and Beskok2011) due to their oversimplification of fluid–surface interactions.

Momentum and energy transfer between the fluid and solid molecules is modulated by the local density of the adsorbed fluid, which depends on the fluid–surface interaction strength. Notably, the present kinetic model captures non-hydrodynamic effects at the fluid–solid interface, predicting a slip velocity of $u_s\simeq 42$ m s $^{-1}$ and a temperature jump of $T_j\simeq 11.3$ K, both closely matching MD simulation results. The model also accounts for viscous heating effects, as shown in the inset of figure 4(b), though this temperature distribution is less discernible in MD simulations due to thermal fluctuations.

Figure 4(c,d) displays fluid behaviours at $H=5.0$ nm, $T_w=233$ K, $|\boldsymbol{a}_e|=8.382\times 10^{11}$ m s $^{-2}$ and $\rho _{avg}=500$ kg m $^{-3}$ . Once again, the present molecular kinetic model demonstrates excellent agreement with MD data, showcasing its modelling capability for NCFs across a broad range of densities.

Traditional boundary conditions attempt to encapsulate both the intrinsic and apparent non-hydrodynamic effects arising from fluid–solid interactions. However, the thickness of the adsorbed layers in NCFs, which varies with the solid type and thermodynamic conditions, is non-negligible (Hadjiconstantinou Reference Hadjiconstantinou2024). Therefore, the adsorbed layer effects cannot be captured by boundary conditions. Our molecular kinetic model consistently considers molecular interactions between fluid and solid surface molecules, integrating both intrinsic and apparent interfacial dynamics within a unified transport equation and thus more suitable for predicting NCFs behaviours.

4.3. Knudsen minimum

In the Poiseuille flow of a rarefied gas, the mass flow rate, defined as

(4.1) \begin{equation} Q_n=u_m\int \frac {\rho u}{\rho _{avg}|\boldsymbol{a}_e|H^2}\mathrm{d}y, \end{equation}

typically exhibits a characteristic minimum as the Knudsen number varies, a phenomenon known as the Knudsen minimum (Cercignani Reference Cercignani1988). This behaviour arises from the interplay between viscous dissipation within the fluid and velocity slip at the fluid–solid interface (Wu et al. Reference Wu, Liu, Reese and Zhang2016; Sheng et al. Reference Shan, Wang, Zhang and Guo2020). The Knudsen number is defined as

(4.2) \begin{equation} Kn=m[\sqrt {2} \rho _{avg} \pi \sigma _f^2\chi (\rho _{avg}) H]^{-1}. \end{equation}

As the fluid density increases, the Knudsen minimum disappears. This behaviour was first reported by Guo, Zhao & Shi (Reference Guo, Zhao and Shi2006) for confined fluids and later investigated for dense fluids (Wu et al. Reference Wu, Liu, Reese and Zhang2016; Corral-Casas et al. Reference Corral-Casas, Li, Borg and Gibelli2022). However, these studies do not fully account for the complex interactions between fluid and vibrating solid molecules at the interface.

Figure 5(a) illustrates the mass flow rate at $T_w=283$ K, where we focus on isothermal systems with the momentum accommodation coefficient $\omega _m = 0.30$ and $0.82$ for the smooth and rough surfaces, respectively. The predicted mass flow rates show satisfactory agreement with MD simulations across a wide range of Kn, from 0.039–1.025. However, obtaining MD data for Kn ${\gt } 1.025$ is challenging due to significant thermal fluctuations. An increase in the fluid–solid accommodation coefficient causes a transition from specular to diffuse reflections, as shown in figure 1(d,e). Therefore, a larger accommodation coefficient normally corresponds to a more efficient momentum transfer, representing a rougher channel surface. Comparing results obtained for different accommodation coefficients ( $\omega _m=0.30$ and $0.82$ ) highlights the enhancement of the mass flow rate resulting from specular reflections in nanoscale flows confined by atomically smooth surfaces. This enhancement reflects the rapid mass transport enhanced by specular reflection, as observed in MD simulations and experimental studies (Keerthi et al. Reference Keerthi2018; Qian, Wu & Wang Reference Qian, Wu and Wang2023). Note, as the flow velocity in nanochannels is usually small, the thermal dissipation may be neglected, which leads to simplification of the model, i.e. the evolution of energy distribution function in (3.1), can be ignored.

Unlike the classical Knudsen minimum, typically observed at Kn $\sim 1.0$ , our results reveal a terrace regime for $0.2\lt$ Kn ${\lt } 1.0$ in a 2 nm channel. This terrace-like behaviour indicates a balanced competition between the decreasing rate of viscous dissipation and the increasing rate of surface flow resistance within this range.

To elucidate the impact of fluid–surface interactions and real fluid effects on NCFs behaviours, we single out their respective contributions from the present molecular kinetic model. When $H_{att} \ll H$ , the fluid–surface interaction term $Q_{\textit{fs}}^{lg}$ becomes negligible and the hydrodynamic wall coincides with the physical boundary. In this limit, our model simplifies to the Enskog–Vlasov (EV) model (Shan et al. Reference Shan, Su, Gibelli and Zhang2023b ). For consistency, we use $\omega _m=1.0$ in our model and apply the diffuse boundary condition in the reduced EV model. The difference in predictions between the present and EV models highlights the impact of fluid–surface interactions, which uniformly reduce the mass flow rate across all flow regimes by shortening the effective flow path and thereby increasing flow resistance. Although this enhanced resistance is a characteristic of NCFs, the overall flow rate can still be significantly increased by the large velocity slip occurring at smooth solid surfaces, an effect driven by the specular reflection of fluid molecules (Lichter et al. Reference Lichter, Martini, Snurr and Wang2007; Martini et al. Reference Martini, Roxin, Snurr, Wang and Lichter2008a ; Keerthi et al. Reference Keerthi2018), as quantified in figure 5(a).

At low densities, the terms $Q_{\textit{ff}}^{lg}$ and $Q_{\textit{ff}}^1$ describing real fluid effects can be neglected. Therefore, the EV model reduces further to the Shakhov model equation for a perfect gas. The simulation results are in good agreement with the solution of the Boltzmann equation (Ohwada, Sone & Aoki Reference Ohwada, Sone and Aoki1989), as shown in figure 5(b). Real fluid effects, including the finite volume of fluid molecules and intermolecular attractions, reduce the mass flow rate, especially in the cases with Kn ${\lt } 0.1$ and strong confinement (see figure 5). In contrast, when Kn ${\gt } 1.0$ , the flow rate is only marginally affected by real fluid effects due to the insignificance of interactions between fluid molecules.

Figures 6(a) and 6(b) demonstrate that higher temperatures result in more pronounced mass flow rates, a behaviour that contrasts with predictions from the Boltzmann equation, where the temperature has no effect on the flow rate. This discrepancy arises because (i) real fluid viscosity deviates from the relation $\mu \propto T^{\frac {1}{2}}$ for a perfect gas (Shan et al. Reference Shan, Su, Gibelli and Zhang2023b ), and (ii) the fluid–surface accommodation rates vary with temperature (Cercignani Reference Cercignani1988). Notably, while the Knudsen minimum disappears at high temperatures (283 and 383 K) in 2 nm nanochannels, it re-emerges at a lower temperature (183 K), where the convective flow is enhanced and the velocity slip is suppressed. This finding reinforces the explanation of how the Knudsen minimum vanishes for confined dense/real fluids.

Fluid adsorption is macroscopically governed by surface wettability (Fan et al. Reference Fan, Wu and Wang2020a ) and quantified by the contact angle, which is directly related to the fluid–solid interaction strength $\epsilon _D$ (Barrat & Bocquet Reference Barrat and Bocquet1999; Perumanath et al. Reference Perumanath, Chubynsky, Pillai, Borg and Sprittles2023). As shown in figures 6(c) and 6(d), as $\epsilon _D$ decreases, the mass flow rate increases significantly, reflecting the transition from hydrophilic to hydrophobic surface conditions. Under hydrophobic conditions, the lower fluid density near the surface and weaker fluid–surface attraction result in larger velocity slips, thereby leading to more pronounced flow rates. Therefore, our model captures composition-dependent non-hydrodynamic effects observed in MD simulations (Martini et al. Reference Martini, Roxin, Snurr, Wang and Lichter2008b ; Hadjiconstantinou Reference Hadjiconstantinou2024).