Mesoscopic method to study water flow in nanochannels with different wettability

Molecular dynamics (MD) simulations is currently the most popular and credible tool to model water flow in nanoscale where the conventional continuum equations break down due to the dominance of fluid-surface interactions. However, current MD simulations are computationally challenging for the water flow in complex tube geometries or a network of nanopores, e.g., membrane, shale matrix, and aquaporins. We present a novel mesoscopic lattice Boltzmann method (LBM) for capturing fluctuated density distribution and a nonparabolic velocity profile of water flow through nanochannels. We incorporated molecular interactions between water and the solid inner wall into LBM formulations. Details of the molecular interactions were translated into true and apparent slippage, which were both correlated to the surface wettability, e.g., contact angle. Our proposed LBM was tested against 47 published cases of water flow through infinite-length nanochannels made of different materials and dimensions―flow rates as high as seven orders of magnitude when compared with predictions of the classical no-slip Hagen-Poiseuille (HP) flow. Using the developed LBM model, we also studied water flow through finite-length nanochannels with tube entrance and exit effects. Results were found to be in good agreement with 44 published finite-length cases in the literature. The proposed LBM model is nearly as accurate as MD simulations for a nanochannel, while being computationally efficient enough to allow implications for much larger and more complex geometrical nanostructures.

Figure 1

Schematic diagram of water flow in a nanochannel [51]. The strength of the water-solid attraction weakens from (a) (hydrophilic) to (c) (hydrophobic). (a) The water-solid interaction is larger than the water-water interaction. The density and viscosity of the blue-square region is higher than the water in the bulk-water region, inducing a negative apparent slip length $l_{s,a}$. (b) The water-solid interaction is comparable to the water-water interaction. The density and viscosity of the green-arrow region is smaller than that of the bulk-water region (green arrows), owing to depletion. Apparent slip $l_{s,a}$ occurs too, but the value is positive. (c) The water-solid interaction is smaller than the water-water interaction. Similar to B, the density and viscosity of the green-arrow region is smaller than that of the bulk-water region (with green arrows.) In addition to the apparent slip $l_{s,a}$ induced by density and viscosity, the true slip $l_{s,t}$ occurs owing to the direct surface sliding of water molecules on the wall. Both of the two slip-length values are positive here.

Figure 2

(a) Typical representation of water molecules near the wall surface in a nanochannel (>1.60 nm) [74]. Water molecules (blue circles) tend to be layered and ordered near the wall surface. The fluctuating green curve is the water-density profile. The first liquid layer starts to develop after the yellow-filled region near the wall surface. The term ${\delta }_1$ is defined as the distance between the solid surface and the location of the first density peak, and ${\delta }_2$ is defined as the distance between the location of the first density peak and the second density peak. The value of stand-off distance δ is generally taken as 0.24 nm for water, and ${\delta }_1$ as well as ${\delta }_2$ are about one molecular diameter, 0.28 nm for water (Appendix pp3). The pink curve is the velocity profile, and no velocity can be observed within the stand-off distance because water density is zero in the region. (b) Determination of the relationship between the first density peak and the contact angle. Data are from MD simulations in the literature [37, 43, 47, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88]. Detailed parameters shown in Appendix pp6. Blue dashed line is optimally fitted curve.

Figure 3

(a) Relationship between the tuning parameters ($G_b$ and $G_w$) and wettability. (b) Density profiles under different contact angles in proximity of the wall surface. For the hydrophilic wall surface (e.g., $\theta =0^{\circ }$), the water density is increased in the near-wall region; for the hydrophobic wall surface (e.g., $\theta ={170}^{\circ }$), the density is decreased in the near-wall region. (c) Velocity profiles under different contact angles in proximity of the wall surface. The colorful solid lines are the velocity profiles. The red and blue dashed arrows are the extrapolated velocity for $\theta =0^{\circ }$ and $\theta ={170}^{\circ }$, respectively. The corresponding apparent slip lengths are $l_{s,a1}$ and $l_{s,a2}$. As shown, $l_{s,a1}$ is a positive slip length, whereas $l_{s,a2}$ is a negative slip length. The wall interface is located in the first grid. Notably, to decouple the contribution of apparent slip and true slip, only the apparent slip after incorporation of Eqs. (10) and (12) is tested here, and the true slip length from Eq. (7) is zero. (d) Density profiles in different radial grid numbers ($\theta ={160}^{\circ }$). At a fixed grid resolution (one molecular diameter), the water confined flow in different channel sizes with the range of 2.72 nm ($N_y=10$) to 83.92 nm ($N_y=300$) is simulated.

Figure 4

Confined water flow in an infinite-length nanochannel. (a) Schematic of simulated domain with periodic boundary conditions (BC) on both sides. (b) Example of LBM simulation results for confined water flow in an infinite-length nanochannel ($\theta ={160}^{\circ }$). Fluid density (lattice unit) near the wall surface is decreased in the hydrophobic tube by molecular interactions. (c) Comparison of experiments or MD simulation-enhancement factor ε from literature and simulation-enhancement factor ε by our proposed LBM under identical conditions (infinite-length nanochannel). Comparison includes 47 cases from literature, in which 23 cases are experiments [35, 42, 99, 100, 101, 102] and 24 cases are MD simulations [98, 103, 104, 105, 106, 107]. The comparison shows most cases from different research can be recovered well. The black diagonal line guides the reader's eye to the comparison. Error bars in enhancement factor ε originate from uncertainties in determination of contact angles or tube size.

Figure 5

Comparison of the density profiles from MD simulations (gray dashed line) in the literature and the one by the proposed LBM (red dots) under the same conditions. The nine cases cover the wettability from hydrophilic to superhydrophobic channels with different dimensions and materials [43, 75, 76, 77, 78, 79, 80], and all the comparisons show the proposed mesoscopic LBM can reproduce the continuous density profile in a molecularly coarse-grained way. The contact angle, channel size, and LBM grids of the simulated cases are (a) $\theta ={13}^{\circ }$, $d=2.5\mathrm{nm}$, $N_y=9$ [75]; (b) $\theta ={30}^{\circ }$, $d=10\mathrm{nm}$, $N_y=36$ [80]; (c) $\theta ={52}^{\circ }$, $d=4\mathrm{nm}$, $N_y=15$ [43]; (d) $\theta ={71}^{\circ }$, $d=5\mathrm{nm}$, $N_y=18$ [76]; (e) $\theta =86.1^{\circ }$, $d=7\mathrm{nm}$, $N_y=25$ [79]; (f) $\theta ={120}^{\circ }$, $d=2.5\mathrm{nm}$, $N_y=9$ [75]; (g) $\theta ={154}^{\circ }$, $d=10\mathrm{nm}$, $N_y=36$ [78]; (h) $\theta ={163}^{\circ }$, $d=10\mathrm{nm}$, $N_y=36$ [78]; (i) $\theta ={178}^{\circ }$, $d=10\mathrm{nm}$, $N_y=36$ [78].

Figure 6

Comparison of the velocity profiles from MD simulations (gray hollow dots) in the literature and the one by the proposed LBM (colored dots) under the same conditions. The five cases cover the wettability from hydrophilic to hydrophobic channels with different dimensions and materials [25, 81, 84, 108], and all the comparisons show the proposed mesoscopic LBM can reproduce the velocity profile in a molecularly coarse-grained way. The external forces in these cases are shown in a form of MD in the original literature. To avoid the unit conversion and achieve the fitting, we assume the external force is 1 MPa/nm in LBM simulations, and the simulated LBM velocity profiles are scaled by the center velocity of MD to obtain the final values. The contact angles of the simulated cases are (a) $\theta ={43}^{\circ }$ [81]; (b) $\theta =105.5^{\circ }$ [25]; (c) $\theta =130.5^{\circ }$ [25]; (d) $\theta ={150}^{\circ }$ [84]; (e) $\theta =161.5^{\circ }$ [108].

Figure 7

Confined water flow in a nanochannel with end effects (finite-length nanochannel). (a) Schematic of simulated domain with two water baths connected on both ends. Rationalization of value $n$ is calibrated in Appendix pp5. (b) Example of LBM simulation results for confined water flow with ends ($\theta ={160}^{\circ }$). Bending of streamlines (yellow arrowed lines) at entrance and exit induces significant additional flow resistance, compared with flow through infinite nanochannels. Density is lattice unit. (c) Comparison of experiments or MD simulation enhancement factor ε from literature and simulation enhancement factor ε from our proposed LBM under identical conditions (finite-length nanochannel). Comparisons, including 50 cases from MD simulations [25, 36, 37, 108, 113]. Comparison shows most cases from different studies can be matched well. The black diagonal line guides the reader's eye to the comparison. Error bars of enhancement factor ε originate from uncertainties in determination of contact angles or external force.

Figure 8

Sketch of the linear (a) (Couette flow) and (b) nonlinear velocity distribution (channel flow) with Navier linear slip assumption in the boundary [57].

Figure 9

Determination of the distance between the solid surface and the location of the first density peak ${\delta }_1$. The data are from MD simulations in the literature, and the simulated wall surfaces are (1) flexible and rigid carbon nanotubes (CNTs) with various surface wettability (Tao et al. [75]); (2) flat walls with two graphene layers (Neek-Amal et al. [76]); (3) Cu and graphene-coated Cu surface wall (Vo et al. [77]); (4) Lennard-Jones hydrophobic substrate (Evans and Wilding [78]); (5) Au surface wall (Vo et al. [79]); (6) double-walled carbon nanotubes (DWCNTs) (Walther et al. [37]); (7) self-assembled monolayers consisting of single membranes of either ${\mathrm{C}}_{20}{\mathrm{H}}_{41}\mathrm{OH}$ or ${\mathrm{C}}_{20}{\mathrm{H}}_4$ (Bonthuis and Netz [128]); (8) sodium-saturated Wyoming montmorillonite surface to represent clay (Botan et al. [81]); (9) hydrophobic hydrogen terminated diamond surface (Sendner et al. [43]); (10) Lennard-Jones substrate (Chen et al. [83]); (11) armchair CNTs (Thomas and McGaughey [84]); (12) two (111) planes with an face-centered cubic (FCC) lattice consisting of 12,288 molecules (11)(Priezjev et al. [129]); (13) zigzag CNTs (Walther et al. [86]); (14) CNTs (Werder et al. [87]); (15) CNTs (Walther et al. [88]). The blue dashed line with the value of 0.28 nm is drawn to guide the reader's eyes.

Figure 10

Determination of the distance between the location of the first density peak and the second density peak ${\delta }_2$. The data are from MD simulations in the literature, and the simulated wall surfaces are (1) flexible and rigid CNTs with various surface wettability (Tao et al. [85]); (2) flat walls with two graphene layers (Neek-Amal et al. [74]); (3) Si(111) plane (Ramos-Alvarado et al. [76]); (4) Cu and graphene-coated Cu surface wall (Vo et al. [77]); (5) Lennard-Jones hydrophobic substrate (Evans and Wilding [78]); (6) Au surface wall (Vo et al. [79]); (7) double-walled carbon nanotubes (DWCNTs) (Walther et al. [37]); (8) self-assembled monolayers consisting of single membranes of either ${\mathrm{C}}_{20}{\mathrm{H}}_{41}\mathrm{OH}$ or ${\mathrm{C}}_{20}{\mathrm{H}}_4$ (Bonthuis and Netz [S10]); (9) sodium-saturated Wyoming montmorillonite surface to represent clay (Botan et al. [81]); (10) hydrophobic hydrogen terminated diamond surface (Sendner et al. [43]); (11) armchair CNTs (Thomas and McGaughey [84]); (12) two (111) planes with an FCC lattice consisting of 12,288 molecules (Priezjev et al. [129]); (13) zigzag CNTs (Walther et al. [86]); (14) CNTs (Werder et al. [87]); (15) CNTs (Walther et al. [88]); (16) Lennard-Jones hydrophobic and hydrophilic walls (Allen et al. [130]). The blue dashed line with the value of 0.28 nm is drawn to guide the reader's eyes.

Figure 11

$P$-ρ plot of the Shan-Chen equations of state. Dashed horizontal line is at $P=0$.

Figure 12

Relationship between enhancement factors and $n$ under different contact angles ($d=2.034\mathrm{nm}$).