Self-diffusivity of dense confined fluids

Fluids confined to nanoscale geometries are ubiquitous in biological and engineering applications,^1–3 such as in nanocatalyst construction,⁴ nanostructured water filtration membranes,⁵ and proton transport within fuel cells.⁶ A notable application that motivates this research is the transport of hydrocarbon fluids stored within sedimentary shale rocks,^7–9 where the gas–gas and gas–surface physics through the naturally formed nanopores deviate from our current understanding and modeling capabilities.

The behavior of confined fluids is still poorly understood, for several reasons. The first challenge appears when the mean free path (MFP) λ, which is the distance traveled by fluid molecules between collisions, is comparable to the flow characteristic length H (e.g., the pore height), and so the Knudsen number $Kn=λ/H$ is moderately larger than zero. In this case, the assumption of local thermodynamic equilibrium breaks down,¹⁰ and subsequently the continuum description of fluids fails.¹¹ One of the most recognizable consequences of non-equilibrium effects is the so-called “Knudsen paradox,” which shows up when a constant pressure difference drives a rarefied flow along a narrow straight channel. The mass flow rate displays a characteristic minimum as the inlet pressure is reduced,¹² being in sharp contrast with the monotonic decrease predicted by the Navier–Stokes equations.

A second challenge occurs when the pore characteristic length is similar to the diameter of fluid molecules σ, as typically occurs in shale reservoirs, where pore sizes are very tight with mesopores (⁠$H∼$ 20–200 Å) and micropores (⁠$H<$ 20 Å) dominating the storage of supercritical methane.¹³ In these conditions, molecular ordering happens inside the pores and the effect of dense fluid packing is magnified. Surprisingly, the recent measurement of Poiseuille mass flow rates through ultratight pores, which are characterized by small confinement ratios $R=H/σ$⁠, was found to follow a monotonic increase with Knudsen number.¹⁴ Although the Boltzmann equation predicts the existence of the Knudsen minimum^15–17 and it may explain its disappearance in bent¹⁸ and short¹⁹ channels, it cannot be used to investigate the effects arising due to confinement, which are responsible for the Knudsen minimum disappearance observed in straight channels.

A third challenge arises when one considers the high pressures in these unconventional reservoirs. Pressures can range from 5 to 50 MPa,²⁰ and consequently fluid densities are beyond the dilute-gas Boltzmann limit.²¹ Therefore, the space correlations of molecules, which are disregarded in the derivation of the Boltzmann equation, need to be considered. The Enskog equation approximately considers that the positions of molecules are statistically correlated in dense fluids, and recently it has been shown to predict the Knudsen minimum disappearance within straight channels of molecular dimensions.²² However, a satisfactory explanation of the underlying mechanisms that justify these mass flow rates is missing, particularly in the transition regime (⁠$0.1<Kn<10$⁠).

A possible theoretical explanation of the Knudsen minimum disappearance in straight channels relies on a change of the underpinning physics, where mass transport is no longer occurring because of collective advection. In pressure-driven flows within very tight geometries, there is not enough room for molecules to develop different velocities based on their spatial position, because of the limited number of molecular layers in the confining direction and the large amount of collisions with walls. Therefore, the fluid velocity profile switches from the Poiseuille parabolic shape to plug-like as the characteristic dimension of the channel reduces,^23,24 suggesting that diffusion now dictates the fluid dynamics in this context. Given that the Knudsen minimum only disappears when the confinement is on the molecular scale, it is clear that a deeper understanding of diffusive dynamics is of paramount importance to better explain the fluid behavior, and so this will be the scope of the present work.

The classical modeling approach to diffusive processes is based on Fick's law, where diffusion follows a linear response with the gradient of the species' concentration along the streamwise direction. The simplest diffusive phenomenon, i.e., Brownian motion of identical particles without a net flow induced by a thermodynamic gradient, is referred to as self-diffusion and can be studied by a process of tagging and tracking particles. The constant of proportionality, which is known as the self-diffusivity D, is commonly identified with the mean square displacement (MSD) of molecules in bulk fluids.²⁵ The value of D depends on the physical mechanisms that govern diffusion, e.g., molecular and Knudsen diffusion amongst many others.²⁶ Molecular diffusion dominates in the hydrodynamic limit for Kn <0.01, when interactions are mainly between fluid particles, whereas for Kn >10 the diffusive dynamics can be simplified to Knudsen diffusion, consisting of free flights of particles between collisions with walls.

Although Fick's law is the cornerstone of modeling diffusion transport in porous media,²⁷ its validity at the nanoscale is not obvious and so it has been subjected to many studies.^28,29 The Fickian approach clearly cannot be applied under extremely tight confinements where single-file diffusion takes place.³⁰ This anomalous transport process occurs when molecules move without being able to overtake each other, abiding to their original order in the row, which differs from the normal diffusion processes described by Fick's law. For larger but still moderately confined channels (⁠$R∼10$⁠), this linear theory is assumed to be valid but some doubts arise with respect to the proportionality factor that has to be used within its formulation. Specifically, the non-bulk structure of the fluid may involve the breakdown of the MSD approach for evaluating the self-diffusion coefficient^31,32 and, indeed, the lower R threshold that dictates the validity of the MSD-based self-diffusivity in the Fickian framework is still not well defined.³³ Local diffusion coefficients have also been introduced to cope with the fluid inhomogeneity,^34–36 but the procedures for averaging the local values to compute global quantities are somehow phenomenological in nature.

Molecular dynamics (MD) computer experiments have served as a useful reference for evaluating the self-diffusion coefficient in both bulk^37–40 and confined fluids.^41–44 However, a satisfactory overall picture of the diffusion dynamics of fluids under confinement is still lacking and there are no analytical derivations to predict self-diffusivities accurately in the entire range of Knudsen numbers. Importantly, the crossover between molecular and Knudsen diffusive mechanisms has not been generalized beyond the rarefied description of non-confined gases,^45–47 and as we explain above, this transition flow regime under tight confinement and high fluid packing is an area with critical implications to shale gas media.

The aim of this paper is to systematically investigate the effect of rarefaction, confinement, and fluid-wall friction on the self-diffusivity. The novelty of this study is presented next. First, we measure the self-diffusion coefficient in a wide range of fluid densities (from dense to rarefied), channel heights (from tight to quasi-bulk), and microscopic wall properties (from rough to frictionless) that has not been subjected to investigation before. Here, we use event-driven molecular dynamics (EDMD), which, compared to the Enskog theory, provides the exact description of the monatomic hard-sphere dynamics in the whole range of reduced fluid densities. Second, we assess the validity of using the well-known MSD procedure for studying diffusion in strongly inhomogeneous fluids. Third, we use a splitting procedure of the colliding molecules to understand the interplay between molecular and Knudsen mechanisms in the transition regime, and how these are influenced by confinement and rarefaction. This calls into question previously observed diffusive scaling laws. Fourth, we develop a semi-analytical theory that predicts the self-diffusion coefficient in the wide parametric space of consideration, which has so far been missing in the literature.

The rest of the paper is organized as follows. In Sec. II we outline the EDMD methodology and the procedures for calculating self-diffusivities, from either following particles' trajectories through their random walk or according to Fickian theory. Results for self-diffusivity are presented in Sec. III, including the validation of MSD and Fickian approaches and splitting results of diffusive mechanisms and demonstration of the proposed theory. Section IV concludes with the main findings and future work.

A system composed of N hard-sphere particles with diameter σ is studied in a slit geometry, defined by two infinite parallel plates at a distance $H=h+σ$ apart, where h is the effective transversal space accessible to the center of spherical molecules (see Fig. 1). The walls are assumed to be structureless flat surfaces, and the fluid–wall interactions are described by the Maxwell scattering kernel with tangential momentum accommodation coefficient (TMAC). Namely, a fraction of the molecules impinging on the wall (given by TMAC) is diffusely re-emitted after being thermalized with the wall, while the remaining partition (⁠$1−TMAC$⁠) is specularly reflected. This is the simplest possible simulation setup that permits one to capture the key features of the diffusion process.⁴⁸ 

The self-diffusion coefficient is evaluated by two different methods, as illustrated in Fig. 1: (a) Under equilibrium conditions, the self-diffusivity is computed using the Einstein relation,⁴⁹ measuring the squared deviation of each particle position $ri$⁠, with respect to a reference position over time t, as follows:

Here, d is the dimensionality of the system, and $⟨Δr2(t)⟩$ is double averaged over the number of particles N and multiple time origins τ, where $τ>tj$ for all j. Unlike the bulk case, diffusion is not an isotropic process in confined geometries,⁵⁰ and the self-diffusivity in the $x−$direction takes a lower value than in the yz plane, with d = 2 to compute the mean D^E in these two symmetric dimensions. (b) In the presence of a concentration difference in the streamwise direction, the self-diffusivity is computed by using the phenomenological Fick's first law,⁵¹ i.e., as the proportionality factor between the concentration gradient dC/dy and the flux of particles J that arises as a result of this driving force

Two remarks are in order: first, the Fick's law refers to the linear relation between mass flux and concentration gradient, regardless of the fluid rarefaction state and not only in the continuum regime.⁵² Second, the Fickian self-diffusivity is different from the transport diffusivity that is defined in the presence of a convective flow.⁵³ 

Three dimensionless groups can be identified to systematically describe the different diffusive processes that may take place in the current problem formulation: the Knudsen number $Kn=λ/h$⁠, the confinement ratio $R=h/σ$⁠, and the reduced density $η=nπσ3/6$⁠, where n is the number density. In addition to differentiating between flow regimes, the Knudsen number, combined with the confinement ratio, determines the diffusive mechanism that prevails at a given rarefaction and confinement degrees. The reduced density represents the packing fraction of the fluid, e.g., larger η implies less free space for the particles to move. By using the Enskog theory, these groups are interrelated through

where χ is the contact value of the pair correlation function in a hard-sphere fluid in uniform equilibrium.²¹ An approximate but accurate expression for χ can be obtained from the hard-sphere fluid equation of state proposed by Carnahan and Starling⁵⁴ 

The time evolution of the hard-sphere dynamics is simulated using event-driven molecular dynamics (EDMD),⁵⁵ which is an adaptation of MD simulations to discrete potential systems. This simulation technique is event-driven in the sense that the state of the system jumps from one time to another corresponding to the earliest collision event. The time step is not constant throughout the simulation run, like in MD simulations, as it depends on the spatial coordinates and velocities of all molecules in the system. The algorithm consists of three basic steps: (a) evaluating the time of the earliest collision event, (b) moving ballistically all particles for that time interval, and (c) updating the velocity of the particles that have collided with another particle or the wall, according to hard-sphere dynamics or the Maxwell scattering kernel, respectively.

The main advantage of EDMD relies on the computational savings with respect to simulations that consider soft potentials, like Lennard–Jones, because it avoids the expensive calculation of multibody intermolecular forces. To simplify the notation and data analysis in the rest of the paper, all the physical quantities are made dimensionless by considering the molecular diameter σ as the reference length, the particle mass m as the reference mass, and $σ/kT/m$ as the reference time, where the denominator is the reference thermal velocity, with k being the Boltzmann constant and T the temperature. Accordingly, self-diffusivities will be expressed in dimensionless $σkT/m$ units.

The computational domain is an orthogonal box of dimensions $[H,ly,lz]$⁠, applying periodic boundary conditions in the $y$- and $z$-directions. The cross section of the simulation box $H×lz$ is chosen to match the target reduced density η. Note that a careful choice of the streamwise length l_y is required when evaluating the self-diffusivity through Eq. (2), as will be discussed in Sec. III A. The number of particles was set over $5×103$ in order to reduce the system-size dependence of the computed self-diffusivities.⁵⁶ Initially, molecules are randomly placed across the simulation box ensuring that they do not overlap, and velocities are sampled from the Maxwell–Boltzmann distribution by using the Box–Muller algorithm. The reduced density is varied in the range $η=[0.000 05, 0.35]$⁠, with confinement ratios within $R=[2, 100]$⁠. Consequently, the Knudsen number spans all flow regimes,⁵⁷ i.e., $Kn ∼[0.001, 1000]$ as given by Eq. (3).

As detailed in Sec. II A, the self-diffusion coefficient has been computed using two different approaches. The computation of self-diffusivities by Eq. (1) is performed based on equilibrium simulations. Time averaging starts when the initial ballistic regime dies out and the simulation time is set long enough to ensure that the ergodic hypothesis is satisfied. Note that there are two main sources of uncertainties in evaluating self-diffusion coefficients. First, accurate self-diffusivity results require an unbiased MSD, and therefore uncorrelated sampling⁵⁸ needs to be considered when maximizing the number of samples per run, e.g., multiple time origins t_j besides the initial simulation time t₀. Second, using the Einstein relation is associated with relevant uncertainties, because of the unattainable search of perfectly linear regions between MSD and time.⁵⁹ 

The calculation of the self-diffusivity by Eq. (2) uses the same setup, but particles are tagged as either type A (tracer) or type B (solvent), with all fluid molecules being mechanically indistinguishable. The overall density, accounting all tagged and non-tagged particles, remains constant throughout the domain. Initially, all molecules are assigned to type B. During the simulation, those which cross the right-hand and left-hand side boundaries of the simulation box are assigned to type A and type B, respectively, independently of their initial tag (see Fig. 1). This mimics the presence of two reservoirs, virtually full of tracers and of solvent in both ends of the channel. Molecular motion and collisions are performed in the same manner for both types, but macroscopic quantities such as number density and mean velocities are collected separately for each of them. The concentration profile of type A particles is then allowed to relax, and the time averaging starts at the steady state when the profile becomes linear in the middle of the simulation box, a few molecular diameters away from the boundaries to avoid non-diffusive end effects. The concentration gradient is then computed by minimizing the mean square error of the linear fit to the concentration field of type A particles in this region.

The Einstein and Fickian self-diffusivities match up in the continuum regime,⁶⁰ but their agreement has been questioned in tightly confined geometries, as introduced in Sec. I. Preferential fluid packing next to the walls extends to a larger relative portion of the channel when R values are low. In this condition, most of the MSD trajectories are made up of contributions from particles crossing inhomogeneous areas with different local densities. Accordingly, it is not obvious whether the microscopic Einstein self-diffusivity still agrees with the coefficient coming from the macroscopic Fickian theory. In order to shed light on this question, we compute the self-diffusion coefficients based on Eqs. (1) and (2) for a wide range of η and R, with emphasis on large reduced densities that lead to non-homogeneous density profiles.

A preliminary observation in our results must be made before moving into a detailed comparison. As shown in Table I, the Einstein self-diffusivity D^E does not significantly depend on the length of the domain in the streamwise direction, which is expected since the MSD measurements of particles continue seamlessly across the periodic boundaries, as if to represent an infinitely long slit. On the other hand, this is not the case for the Fickian self-diffusivity, in which D^F is seen to increase steadily with l_y and eventually levels off for $ly≳2λ$⁠, regardless of the confinement ratio and the fluid density. This occurs because D^F is now influenced by the inlet- and outlet-type boundary conditions, and a critical length is required to offset the non-equilibrium perturbation of particles changing tag at the boundaries. We find this requirement on the streamwise length to be particularly severe in the rarefied regime, due to the inverse proportionality between η and λ.

Figure 2 shows the Einstein and Fickian self-diffusivities as a function of the confinement ratio R for different large packing fractions, when walls are fully diffuse (⁠$TMAC=1$⁠). For the aforementioned reasons, the Fickian self-diffusivities are obtained by running simulations where the streamwise length of the computational box is larger than two molecular mean free paths. Our results reveal that the predictions given by Eqs. (1) and (2) are almost indistinguishable in the whole parametric space. Particularly, Einstein and Fickian self-diffusivities agree under tight confinement (⁠$R≲20$⁠) despite the strong fluid inhomogeneities that exist (see the inset of Fig. 2). It is worth stressing that this conclusion may be a consequence of the slit geometry considered in this study, which allows particles to swap places in the spanwise z-direction even in molecular-like confinements (⁠$R≤3$⁠). Hence, diffusion is not anomalous, although there is not enough room for particles to overtake their counterparts through the confining x-dimension. This 2D planar channel feature ensures that the linearity between MSD and time assumed by Eq. (1) always holds, and anomalous diffusion never takes place.

At a given reduced density η, the molecular mobility decreases with narrower channels because, on average, particles collide more frequently with the confining walls. Throughout these collisions, particles experience a change in momentum that hinders their displacement in the two non-constrained $y, z$-directions, parallel to the walls. For this reason, the self-diffusivity shows a decreasing tendency with increasing confinement in Fig. 2, which is logarithmic-like as predicted theoretically by some studies^62,63 (dashed lines are included as a guide to the eye). If the confinement ratio is large (⁠$R≥50$⁠), an excellent agreement is found between our simulations of a dense fluid and the bulk self-diffusion coefficients at the same reduced density, with all deviations below 2.5%. Table I and Fig. 2 represent the validation of the EDMD code, which was developed in house, to accurately compute the self-diffusion coefficient for both confined and bulk-like fluids.

The self-diffusivity results that are presented in the remaining part of the paper are obtained using the Einstein relation. This choice stems from its overall lower computational cost over the Fickian approach, and also because of the convenience of producing self-diffusivities independent of l_y, as shown in Table I. To provide a more comprehensive picture, besides the self-diffusivities reported in Fig. 2, we also show results for the rarefied end of reduced densities (lower η) in Fig. 3. Note that $η=0.01$ can be approximately considered as the threshold between dense and dilute gas, since for lower values of the reduced density the compressibility factor of the fluid deviates from unity (ideal gas behavior) less than 5%. Self-diffusion coefficients are presented as a function of Kn in Fig. 3(a), while the analysis is performed with respect to η in Fig. 3(b).

As expected, self-diffusivities increase (decrease) with Knudsen number (packing fraction) because particles have more room to move freely before colliding with another entity in the system. The Knudsen diffusion mechanism, where Kn goes to infinity because of the zero density limit, is associated with the larger mobility for a given R. The lower self-diffusivity threshold is related to larger reduced densities through the molecular diffusion mechanism, where the fluid enters the metastable regime after freezing at $η=0.494$⁠.⁴⁰ 

Collisions are the main driving mechanism for diffusion, and as such, making a distinction between fluid–fluid (molecular) and fluid-wall (Knudsen) collisions would give us a better understanding of how these are influenced by rarefaction and confinement. Particles are scattered isotropically when they collide with other particles, which is in clear contrast to when they collide with the wall, where particle re-emission occurs toward the half-space occupied by the fluid.

Figure 4 shows the different collision frequencies, defined per unit molecule, between fluid–fluid (⁠$υFF$⁠), fluid-wall (⁠$υFW$⁠), and the total amount (⁠$υFF+υFW$⁠) measured directly from our EDMD simulations in two slits of R = 2 and R = 100. For smaller Knudsen numbers (Kn < 2), fluid–fluid collisions dominate as expected. At $Kn∼2$⁠, we observe a crossover point for both slits at which fluid-wall and fluid–fluid collisions are similar. For larger Knudsen numbers (Kn > 2), fluid-wall prevails with respect to intermolecular collisions. A high confinement has a role in both increasing the number of collisions (⁠$υFW$ and $υFF$⁠), which is attributed to dense effects becoming more prominent, and in narrowing the difference between $υFW$ and $υFF$ in the earlier Knudsen numbers, which in this case is at least one order of magnitude, i.e., $υFW$ becomes more important for small R.

The fluid–fluid collision frequency $υFF$ can be estimated from elementary arguments of kinetic theory of a dense gas as the mean thermal speed $v¯$ over the mean free path λ,

while a different theoretical derivation is involved when evaluating the fluid-wall collision frequency $υFW$⁠, where the characteristic length is now the particle-center accessible region h,

Equation (5) for fluid–fluid collisions agrees very well with our EDMD results, except for slight deviations at low Knudsen numbers in the tighter channel (R = 2), which occur because small density inhomogeneities are apparent in the middle of the channel for ultratight confinements.

Fluid-wall interactions deviate significantly from Eq. (6) at high fluid packing. This occurs because the ordering of fluid density near the wall can be several times higher than the nominal density, on which the theoretical equation has been derived, which leads to an increase in the collision frequency with the wall. For R = 100, this effect of fluid-wall collisions due to layering is negligible on the total number of collisions, as can be seen in Fig. 4, where $υFF≫υFW$⁠.

These findings emphasize the importance of the transition regime (⁠$0.1<Kn<10$⁠), in which there is an interplay of both fluid–fluid and fluid-wall collisions impacting the self-diffusivity, which are also influenced by ordering at low Kn and high confinement, as can be seen in Fig. 4.

A more convenient parameter that we will use to distinguish the interplay between the different collision frequencies is the collision frequency ratio β, defined as

A very good agreement was found between this theoretical approximation and our EDMD results for channels with large R, whereas, as shown in Fig. 5, Eq. (7) overpredicts the actual β with increasing magnitude as R gets smaller. The largest discrepancies are associated with the dense η limit, which produces the lower Knudsen numbers according to Eq. (3). As before, the underpinning reason for this discrepancy is that Eq. (7) assumes that the fluid is homogeneous, which is not the case for high confinements.

In this section, we derive a semi-analytical model that allows us to predict the self-diffusivity results presented in Fig. 3. The diffusion assessment will be based on the β values at different flow regimes which, by definition, must take values in the range $[0,1]$⁠. Small Knudsen numbers, i.e., the continuum regime, are related to β close to unity because of the predominant role of fluid–fluid interactions. On the other hand, large Knudsen numbers in the free molecular regime are associated with small β, which eventually tends to zero when there are just collisions with walls in the ballistic limit. Finally, the analysis for transitional Knudsen numbers considers intermediate β values, accounting for both fluid–fluid and fluid-wall collisions in varying degrees. In each case, we identify the predominant physical mechanism and discuss the appropriate analytical expressions for the self-diffusivity.

For β close to unity, fluid–fluid collisions are the most predominant form of interaction, and the fluid approaches the continuum regime. The classical description of molecular diffusion in the bulk applies,⁶⁴ and the following semi-empirical formula can be used to predict the molecular self-diffusivity of a hard-sphere fluid:³⁹ 

where the fitting parameters are $c1=0.0730, c2=11.6095$⁠, and $c3=−26.9511$⁠, and D_E is the first order approximation of the self-diffusivity according to the Enskog kinetic theory for dense fluids

Note that this semi-empirical approximation accurately reproduces simulation results for $β≥0.9$⁠, i.e., $Kn≤0.1$⁠, showing a maximum 5.7% relative error with respect to our EDMD results when $R≥20$⁠, as emphasized by the dashed lines in Fig. 6. In Eq. (8), the term within brackets is a correction factor that accounts for dense effects not captured by the Enskog theory.

For β close to zero, fluid-wall collisions are the most predominant form of interaction, and the fluid accordingly approaches the free molecular regime. The first order kinetic theory prediction of the Knudsen self-diffusivity in infinite planar channels is given by⁶⁵ 

where $γ=eC=1.781$⁠, with C being the Euler–Mascheroni constant. In Fig. 6, the predictions given by Eq. (10) are the dashed–dotted lines. Note that Eq. (10) underpredicts the self-diffusivity for Kn < 50, with deviations larger than 5%, because intermolecular collisions are not properly considered in this range of Kn while still being relevant, i.e., this can be seen by closer inspection of Fig. 5, where $β∼0.1$⁠.

It is worth noticing that the Knudsen self-diffusivity, as predicted by Eq. (10), scales linearly with the confinement ratio R and shows a logarithmic divergence with respect to the Knudsen number. This self-diffusivity singular behavior can be attributed to the slit geometry.⁶⁶ Indeed, when Kn goes to infinity, there are a growing number of outlier molecules with zero (or close to zero) $x$-velocity components (normal to the walls) that do not suffer any collisions. Therefore, their contribution to the MSD grows linearly in time and eventually causes the self-diffusion coefficient to grow without any bound. The presence of this type of particle behavior is also responsible for the mass flux divergence of a flow driven by a pressure gradient.¹⁵ 

For intermediate β values, both fluid–fluid and fluid-wall collisions are relevant and must be considered. This situation corresponds to the more challenging transition regime, where self-diffusivities take a value in between the molecular and Knudsen results. From Fig. 5, we can roughly estimate that fluid-wall collisions can no longer be neglected at about $Kn∼0.1$⁠, as they represent above 10% of the total number of collisions.

The Bosanquet interpolating formula has been used in the past to deal with bulk and rarefied gases⁴⁵ inside cylindrical capillaries or porous media. In this work, we find that the Bosanquet formula is an acceptable model also for dense, confined fluids, and allows us to predict the self-diffusivity for these intermediate β values,

The rationale underlying the applicability of Eq. (11) is that in the diffusion process, particles perform a random walk composed of steps whose length follow the Poisson distribution, with mean equal to the particle collision frequency. In the transition regime, both fluid–fluid and fluid-wall collisions must be considered, and the corresponding collision frequencies sum up because this constitutes the mean of the sum of two independent Poisson distributions. As the self-diffusivity is proportional to the reciprocal of the collision frequency, Eq. (11) follows this tendency as well.

It is worth noticing that, by construction, the Bosanquet formula agrees exactly with the molecular self-diffusivity in the continuum regime, and with the Knudsen self-diffusivity in the free molecular one. However, using Eq. (10) within Eq. (11) leads to self-diffusivity predictions in poor agreement with numerical results in the transition regime. This can be explained by the fact that Knudsen self-diffusivities given by Eq. (10) encompass two very different effects, namely, the logarithmic divergence due to the slit geometry and the increasing importance of fluid-wall collisions. In order to decouple these two effects and quantify the contribution of fluid-wall collisions (which is the only one that the Bosanquet formula may capture), the Knudsen self-diffusivity has to be defined in a different way, e.g., by phenomenologically extending the expression derived for a cylindrical capillary⁶⁷ to a planar channel configuration

The factor $a1*=4/5π$ is a fitting parameter, which was tuned to match the EDMD data within the early transition regime to capture the initial deviations with respect to the continuum values, and then leads to good predictions of the self-diffusivities across $0≲Kn≲10$ for most of the channel heights considered in our study, when Eq. (12) is included within Eq. (11), as it can be observed in Fig. 6 using solid lines.

More specifically, Eq. (11) provides fairly good predictions of the self-diffusivity except for tight confinement (i.e., $R≤5$⁠), with a relative error below 10%. The deviations for small confinement ratios R can be explained by the poor accuracy of the D_m prediction, as given by Eq. (8). Indeed, the continuum regime is never reached in tight channels because the Kn is still remarkably larger than zero when the freezing reduced density is reached.

From visual inspection of Fig. 3(a), it seems that self-diffusivities scale linearly with R, which would imply that D/R ratios are equivalent at a given Knudsen number, independently on the fluid rarefaction state and channel characteristic length. The collapse of self-diffusivities, which had already been suggested in previous studies,^42,43 may also be inferred from Fig. 7, where we present the EDMD simulation results divided by the confinement ratio at which the simulations were performed. Despite the master curve (representing the scaled results for $R≥5$⁠) looks particularly good on a log –log plot, it is needed to point out that the linear scaling of self-diffusivities with confinement does not actually hold true.

The functional form of Eq. (10) reveals that the linear relationship between D_k and R exists in the free molecular regime, but this is no longer the case in the continuum regime, where fluid–fluid collisions drive the diffusion process. Indeed, the Enskog self-diffusivity D_E, Eq. (9), depends linearly on the confinement ratio R, but the true molecular self-diffusivity D_m, Eq. (8), includes a correction factor that varies with the packing fraction η and, in turn, with the confinement ratio for equal Kn (as emphasized in the inset of Fig. 7). More specifically, at intermediate packing fractions (up to $η∼0.3$⁠), the self-diffusivity is larger than the prediction from Eq. (9) due to hydrodynamic enhancement consequences,⁶⁸ whereas at higher reduced densities, the cage entrapment of particles comes into play causing the actual self-diffusivity to drop.³⁸ Since this property breaks down in the continuum regime, by extension the linear scaling with R cannot be considered as universal throughout the range of Kn.

It is worth noticing that the linear scaling assumption may provide a rough estimate of the self-diffusivity, with deviations within 3% as long as the packing fraction is low enough, e.g., $η≲0.05$⁠. Therefore, it is clear that this simplifying assumption is especially inaccurate for small R, i.e., in shale gas reservoir applications, where the continuum regime can only be reached at a very large fluid density.

The wall roughness is anticipated to affect the self-diffusivity, especially for small β values.⁴⁴ In the Maxwell scattering kernel, the wall roughness is represented by the accommodation coefficient, i.e., from ideal perfectly smooth walls (⁠$TMAC=0$⁠) to more realistic engineering walls (⁠$TMAC=1$⁠). Simulations are performed at different $TMACs=[0,0.05,0.125,0.25,0.5,0.75,1]$⁠, and the resulting self-diffusivities are reported in Fig. 8 as a function of Kn. Results are presented for two different cases, which are representative of loose (R = 100) and tight (R = 5) geometries.

When the fluid is loosely confined, diffusivity results for various TMACs collapse in the continuum regime (⁠$Kn≤0.1$⁠) as expected, whereas this is not the case in the tightly confined scenario. A close inspection of the results for continuum-like Kn reveals an increase in self-diffusivity with decreasing TMAC, which indicates that fluid-wall collisions remain important for low R values in this regime, and cannot be neglected. Note that when molecules are subjected to a specular reflection with the wall, there is no tangential momentum transfer, and so this represents the zero friction limit where the presence of walls is no longer felt by the fluid. Accordingly, the upper limit of self-diffusivities is associated with $TMAC=0$⁠.

As in Sec. III B, here we also derive a semi-analytical formula for retrieving the self-diffusivities as a function of TMAC. Although the Maxwell scattering kernel is used in this study, one cannot expect to use the Smoluchowski prefactor⁶⁹ to rescale the Knudsen self-diffusivity D_k defined through Eq. (12). Indeed, the latter is just an ad hoc definition introduced to permit one to capture the crossover between molecular and Knudsen self-diffusivity using the Bosanquet formula.

A phenomenological approach is used here for extending the self-diffusivity predictions to the case of partial accommodation coefficients. More specifically, numerical experiments suggest that the following functional form for this prefactor is more appropriate to scale the Knudsen self-diffusivity⁷⁰ 

where f = 1.2261 is a fitting parameter and $DkTMAC=1$ is the Knudsen diffusivity obtained for fully diffuse walls, as introduced in Eq. (12).

In Fig. 8, we show that Eq. (13) provides an accurate estimate of the self-diffusivity of a fluid confined within walls of different roughness, as a function of TMAC. For the loose confinement cases, where Eq. (8) accurately reproduces the diffusive dynamics in the continuum regime, the average relative error of the Bosanquet prediction is below 3.2% for the different accommodation coefficients in our assessment. For tighter geometries, where the influence of fluid-wall collisions is more notorious, the transition between molecular and Knudsen self-diffusivities is also successfully predicted at different TMAC values, when defining D_k as in Eq. (13).

It should be mentioned that Eq. (13) predicts an infinite Knudsen self-diffusivity for fully specular walls, i.e., $TMAC=0$⁠. For this case, $1/Dk=0$ and this term can be neglected from Eq. (11), giving D = D_m for all Knudsen numbers and confinement ratios. This prediction is represented by the black solid lines in Fig. 8, showing the maximum self-diffusivity at a given rarefaction state, i.e., bulk behavior of the fluid with no confinement influence.

We have carried out a comprehensive and fundamental study of the self-diffusion process for a hard-sphere fluid confined between two parallel infinite walls, with data measured from EDMD simulations. Three dimensionless groups were adopted to characterize the self-diffusion process, namely the Knudsen number, the confinement ratio and the reduced density. The tangential momentum accommodation coefficient was used to define the roughness of the wall.

We found that the Bosanquet formula was able to provide a very satisfactory prediction of the self-diffusivity from the continuum to the early free molecular regime (⁠$Kn≤10$⁠), beyond moderately low confinements (R > 5) and from smooth to rough surfaces (all TMACs). Importantly, this work provides insights and predictions of self-diffusivity into the challenging transition regime, that was missing in the literature. A splitting procedure of the colliding particles has identified the interplay between the underpinning diffusive mechanisms, namely, molecular and Knudsen diffusion for dense, confined fluids. Some deviations between the theory and the results for extremely tight geometries (⁠$R≤5$⁠) were observed, which can be explained by the inability of the fluid to show a continuum behavior, i.e., fluid-wall collisions cannot be neglected even at densities close to the fluid freezing point that corresponds to the smallest attainable Kn.

In summary, we can now predict self-diffusivity for dense hard sphere fluids confined in two-dimensional slit geometries for any roughnesses, using Eqs. (8), (9), and (11)–(13) for $0<Kn≤10$ and Eq. (10) for Kn >10.

In the validation part of our study, we have analyzed the suitability of using the widely established Einstein relation for describing the diffusion of strongly inhomogeneous fluids. Despite molecules traverse regions with different local densities, the MSD-based approach yields the same values as the Fickian benchmark model, and therefore, the Einstein self-diffusivity is valid for predicting the diffusive flux even at the molecular scale.

Finally, we have also assessed previous work on the scaling of self-diffusivities with confinement. It was concluded that this relationship can only be considered approximately linear when the fluid is sufficiently rarefied (⁠$η≤0.05$⁠) or whether it is confined inside fully specular (frictionless) walls, none of them being appropriate for shale gas flows.

This work lays the foundation to a more precise modeling of dense, confined flows for engineering applications, where there is an interplay between diffusive and advective processes, and to build on this for more realistic surfaces.

This work was financially supported by King Fahd University of Petroleum and Minerals (KFUPM), Saudi Arabia. L.G. and M.K.B. are funded by the Engineering and Physical Sciences Research Council (EPSRC) under Grant Nos. EP/N016602/1, EP/R007438/1, and EP/V012002/1.

The data that support the findings of this study are openly available in Edinburgh DataShare repository at https://doi.org/10.7488/ds/3111, Ref. 71