Direct matching between the flow factor approach model and molecular dynamics simulation for nanochannel flows

Mathematically formulating nanochannel flows is challenging. Here, the values of the characteristic parameters were extracted from molecular dynamics simulation (MDS), and directly input to the closed-form explicit flow factor approach model (FFAM) for nanochannel flows. By this way, the physical nature of the simulated system in FFAM is the same with that in MDS. Two nano slit channel heights respectively with two different liquid-channel wall interactions were addressed. The flow velocity profiles across the channel height respectively calculated from MDS and FFAM were compared. By introducing the equivalent value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\Delta_{im} } \mathord{\left/ {\vphantom {{\Delta_{im} } D}} \right. \kern-\nulldelimiterspace} D}$$\end{document}Δim/D, FFAM fairly agrees with MDS for all the cases. The study values FFAM in simulating nanochannel flows.

Nanofluidics have been developed fast. They have the important applications in super purification, hemofiltration, DNA analysis, biosensors, drug delivery, heat and mass transfer, and micromachines etc. [1][2][3][4][5][6][7][8] . The studies on them have been mainly experimental or by molecular dynamics simulation (MDS) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] . In these devices, the nanochannel flow has been shown by MDS to be distinctly different from the classical recognition [9][10][11][12][13][14][15] . It may be much slower or much faster than described by the classical hydrodynamic flow theory, depending on the fluid-channel coupling. The rheological properties including the viscosity and density of the confined fluid in a nanochannel are usually evolved due to the fluid-channel wall interaction; they are largely different from the fluid bulk values [16][17][18] . Also, the wall slippage easily occurs in a nanochannel 19,20 . It is believed that these new phenomena should be the important factors causing the very special nanochannel flow. It was also recognized that the effects of the discontinuity and inhomogeneity of the fluid across the channel height contribute very significantly to the nanochannel flow 21 . It was however felt a pity that no powerful constitutive models were available for efficiently analyzing the nanoscale flows in the realistic nanofluidics with the length on the micrometer or millimeter scales in one or two dimensions. This resulted in the lacking of the development of the systematic design theory for a lot of nanofluidic devices such as nanoporous filtration membranes and micro/nano bearings. The researches on these devices were mainly focused on the experimental manufacturing and tests but lacked intensive theoretical understanding [2][3][4]8 .
Molecular dynamics simulation has been popularly used in the study of nanochannel flows for capturing the local physical property evolution and the flow velocity profile of the confined fluid [9][10][11][12][13][14][15][22][23][24][25][26][27][28][29][30][31][32][33][34][35] . It can reveal the flow phenomena in a nanoscale zone; but it is hardly applicable for an engineering flow with macroscale sizes, because of the cost of both the over long computational time and the over large computer storage. It is mandatory to develop the feasible flow models for engineering nanochannel flows; they can be used in the theoretical research and in the design of relevant nanofluidics.
Zhang 36 proposed the flow factor approach model (FFAM), an equivalently continuum model for nanochannel flows, for effectively analyzing the nanoscale flows in very small engineering surface separations. He considered both the discontinuity and inhomogeneity of the fluid across the channel height. He used the concept of the local viscosity, which is dependent on the local density; these two parameters both are dependent on the separation between the neighboring fluid molecules across the channel height. The coarse-graining or multiscale simulation methods were also proposed for saving the computational time and computer storage when simulating the nanochannel flows. Bhadauria and Aluru 37 proposed a quasi-continuum hydrodynamic model for the nano slit pore flow; they computed the density profile and the viscosity change across the channel height; they also incorporated the wall slippage effect by using the Dirichlet type slip boundary condition. Ghorbanian [42][43][44][45][46] have examined this model by comparing it with MDS by using the characteristic parameter values which were ever suspected by some readers to be fitting, although satisfactory agreements between these two approaches were existing. As a new corroboration work, the present paper presents the direct comparison between FFAM and MDS in the calculated flow velocity profiles, by directly inputting the values of the characteristic parameters extracted from MDS to FFAM. By this way, the averaged local density and local viscosity profiles across the channel height in FFAM are the same with those in MDS; in the comparison the physical structure of the simulated system in FFAM is forehand set as the same with that in MDS. The comparison shows that good agreements between these two approaches are existing. The present study not only directly substantiates FFAM but also confirms the earlier work 36,[42][43][44][45][46] .
The importance of the present work is to provide the base stone for the theoretical derivation of the flow equation for the nanoscale flow which is based on FFAM, though this equation itself matched well with the MDS results 21,43 .

The flow factor approach model for nanochannel flow
For a better understanding, FFAM is introduced in this section. This model has been shown in details in Refs. 21,36,42 . Here, it is only briefly repeated. Figure 1 shows the flow factor approach model for the flow in a nano slit pore, where the molecules of the fluid with the equivalent number n across the channel height are ordered normal to the channel wall; because of the interaction between the fluid and the channel wall, the separation between the neighboring fluid molecules across the channel height is varied; this results in the density and viscosity variations across the channel height.
According to the model, in the pressure driven flow shown in Fig. 1, the velocity of the fluid molecule is 21,36,42 : Here, i is the order number of the fluid molecule across the channel height, u a is the velocity of the (n − 1)th fluid molecule (on the upper solid wall), u b is the velocity of the 0th fluid molecule (on the lower solid wall), p is the driving pressure, D is the fluid molecule diameter, x is the coordinate shown in Fig. 1, and η line,i−1 and i−1 are respectively the local viscosity and the separation between the ith and (i − 1)th fluid molecules across the channel height. It is formulated that: and In this model, it is assumed that the distribution of the separation between the neighboring fluid molecules across the channel height is symmetrical with respect to the median plane of the channel; it is taken that www.nature.com/scientificreports/ Here, q 0 and m are respectively constant, and n is an odd number ( n ≥ 5 ). The values of � i+1 /� i and η line,i /η line,i+1 both may be varied across the channel height; the model takes q 0 as the average value of � i+1 /� i (for i = 0, 1,…, (n − 1)/2-2).
As the flow in FFAM is essentially molecular-scale and non-continuum, the dimension of the channel height in FFAM should be no more than the total thickness of the physical adsorbed layers respectively freely formed on both the upper and lower channel wall surfaces. Equation (1) shows that for no wall slippage i.e. u a = u b = 0 , once the values of D, m, n, q 0 , � im /D , η line,im and ∂p/∂x are known, the flow velocity u i of each fluid molecule across the channel height can be directly calculated 36,42 . These input parameters are the characteristic parameters in FFAM. The values of D, n and ∂p/∂x can be determined beforehand. While, the values of m, q 0 , � im /D and η line,im can be calibrated by MDS or found/deduced from other scientific data sources. The natures of both the flowing liquid and the channel wall determine the values of all these characteristic parameters except ∂p/∂x . By setting the values of these characteristic parameters in FFAM as the same with those in MDS, the physical nature of the simulated system in FFAM is the same with that in MDS.

Molecular dynamics simulation
Simulation method and calculated parameters. For comparison with FFAM, molecular dynamics simulation was made for the nano slit pore flow. This section presents the MDS details. Here, the methane nanochannel flow was simulated by the non-equilibrium MDS. The periodic boundary was implemented to the simulation system except for the confined methane molecules in the z-direction. The methane molecules were arranged as a face-centered cubic between two paralled silicon atom flats. The modification of the optimized potential for the liquid simulation (MOPLS) model 47 was used in this work. It is written as: where σ , ε , ε 0 , q and r ij are respectively the length unit, the depth of the interaction energy well, the permittivity of vacuum, point charge and the distance between the atoms i and j . The subscripts a and b respectively represent C and H. α and β denote the interaction between different methane molecules. Moreover, in the potential model, the bond length is l CH = 1.087 Å , and q C = − 4q H = − 0.572 e ( e = 4.803 × 10 −10 esu ). The intermolecular potential parameters are given in Table 1. Figure 2 shows the interactions among different particles when using the non-equilibrium MDS; it also shows the layers used for calculating the relevant parameters of FFAM. In order to solve the interaction problem between the wall atom and the methane molecule, the non-equilibrium multiscale MDS was applied in this work; the intermolecular force between the wall atom and the methane molecule was computed using the methane coarse-graining potential, coupled with the Lorentz-Berthelot mixing rule nearby the wall (see Fig. 1 in Ref. 48 ). The interaction between the wall atom and the methane molecule was modeled by the following equation: where χ is the interaction strength factor between the wall atom and the methane molecule. In this study, χ = 1.0 and χ = 2.0 were respectively employed. The methane coarse-graining potential parameters ε CG = 1.329 kJ mol and σ CG = 3.645 Å were optimized, by using the relative entropy coarse-graining framework 49,50 in the bulk methane ensemble at the initial condition ρ = 377.15 kg/m 3 and T = 140 K . The details of the methane coarse-graining process have been shown in Ref. 49 . The silicon atom intermolecular potential parameters are: ε S = 1.6885 kJ/mol and σ S = 3.826 Å 51 . An ABAB stacking configuration was used to describe the silicon atom flat walls; moreover, every atom was conducted by the harmonic potential site r eq as follows: Here, r(t) is the atom position at time t , and k w = 72ε S / 2 1 / 3 σ 2 S was obtained by using the second derivative value of the silicon atom intermolecular potential at r = r 0 = 2 1 / 6 σ S . To avoid the effect of the compressibility of the methane flowing through the nanochannel, the external driving body force f x = 0.075 (ε/σ ) was employed in the x direction 52 . The quaternion parameters were applied to calculate the orientation coordinates for the methane (tetrahedral) molecule. The fourth Predictor Corrector method was used for solving the motion equations of the methane molecule. The temperature of the simulation system was controlled by the NVT ensemble with a Nosé-Hoover thermostat. In order to check the simulation system temperature, the thermal kinetic energies per molecule were compared in the y and z directions; they are respectively mv 2 y 2 = k B T 2 and mv 2 z 2 = k B T 2 53,54 ( m is the mass of methane molecule, T is the system temperature, and v y and v z are the velocities in the y and z directions respectively). In this simulation, the reduced units are σ = 4.01 Å and ε k B = 142.87 K ; The simulation step is t = 0.001 ( 1.3 × 10 −15 s ), and the energy cutoff circle radius is r cut = 2.5 σ . In order to obtain the accurate results statistically, the results for the initial 1.0 × 10 6 simulation steps were discarded; the relevant molecule information and parameters for FFAM were calculated in the next 5.0 × 10 5 simulation steps. In order to calculate the local viscosity for FFAM, the velocity profiles of the methane fluid were collected by dividing the distance between the two silicon atom flats into 560 bins. Moreover, the nanochannel dimensions are 12 × 12 × 6 and 12 × 12 × 7 (the unit is σ ) respectively, i.e. the distance between the two silicon atom flats is six or seven times of the methane molecule diameter. To obtain the relevant parameters q 0 , η line,j , � im /D and η line,im of FFAM, we supposed the methane fluid layer is five ( n = 5 ) when the system arrived at the equilibrium state for two cases. These parameters were respectively calculated by the following equations: Here, l is the averaged gap between the ( l − 1)th fluid layer and the lth fluid layer, the subscripts z indicates the position of the fluid molecule in the z direction, D = 0.414 nm (the methane molecule diameter), and 〈 〉 denotes the ensemble average. The local viscosities of the methane fluid were calculated by the Poiseuille flow method as follows 55 : Values of the characteristic parameters in FFAM extracted from MDS. For the direct comparison of FFAM with MDS, the input values of the characteristic parameters in FFAM must come from MDS so that these two models are convincing to be initially completely matching. Tables 2, 3, 4 and 5 respectively show the dimensionless separations between the neighboring fluid molecules across the channel height calculated from MDS; the calculated corresponding values of q 0 are respectively put below these tables. Tables 6, 7, 8 and 9 respectively show the values of the local viscosity ratio η line,j η line,j+1 calculated from MDS indicating the sig- u r(t) − r eq = 1 2 k w r(t) − r eq 2 .  www.nature.com/scientificreports/     In the present study, the above values of m, q 0 , � im /D and η line,im calculated by MDS are directly put into Eq. (1) (FFAM); also the values of D, n and ∂p/∂x in FFAM are chosen as the same with those in MDS. By this way, as FFAM describes, the averaged local density and local viscosity profiles across the channel height in FFAM are the same with those in MDS, and the physical nature of the simulated system in FFAM is the same with that in MDS. These are sufficient as the flow velocity profile across the channel height is determined by the averaged density profile in the median section of the channel height, not by the density profile near the solid walls which might be highly oscillatory 59 .

Comparison of FFAM with MDS
The values of the characteristic parameters in FFAM shown above were input to FFAM 21,36,42 to calculate the flow velocities of the fluid molecules across the channel height. These calculated flow velocities were compared with those directly calculated from MDS for the same case. In FFAM, it was taken that n = 5, and the pressure gradient for the Poiseuille flow was taken exactly the same as in MDS i.e. ∂p/∂x = − 3.10621E + 15 Pa/m. Figure 3a shows the comparison between the flow velocities respectively calculated from FFAM and MDS for the case h = 2.484 nm and χ = 1.0 . In the figure, the legend "FFAM 1" refers to the calculation from FFAM by using im D = 0.2713 . For this case, the value of im D calculated from MDS is actually 0.5118 as shown above; the reason for using the corrected value im D = 0.2713 in FFAM is to satisfy the channel height 2.484 nm, because im D = 0.5118 gives the channel height significantly greater than 2.484 nm in FFAM. The reason for this correction is that in actual cases, not completely like as assumed in FFAM, the fluid molecules are often not exactly ordered normal to the channel wall surface. In FFAM, if we use im D = 0.5118 , for giving the channel height 2.484 nm, the equivalent number of the fluid molecules across the channel height should be n eq = 4.11; if we use n = 5, the value of im D must thus be corrected for giving the channel height 2.484 nm. Figure 3a shows that the velocities calculated by "FFAM 1" are considerably smaller than those calculated from MDS especially in the central region of the channel. We checked that the reason for this discrepancy is that the im D value 0.2713 in FFAM 1 is over lower than that (0.5118) calculated from MDS. Thus for this case, in FFAM, it is not rational to use im D = 0.2713 to calculate the flow velocities. In Fig. 3a, we also present the calculation results from FFAM by using im D = 0.36 , as denoted by the legend "FFAM 2". It is shown that the results calculated by "FFAM 2" fairly agree with those calculated from MDS. The reason for this much better matching should be the closer value of im D used in "FFAM 2" than in "FFAM 1". We examined that in "FFAM 2" the channel height is increased only by 8% as compared to the exact value 2.484 nm.

For the case h = 2.484 nm.
The same circumstances also exist in Fig. 3b, which is for the case h = 2.484 nm and χ = 2.0 . In Fig. 3b, the results calculated by "FFAM 1" are considerably smaller than those calculated from MDS, due to the used im D value 0.2659 over lower than that (0.5282) calculated from MDS; while, the results calculated by "FFAM 2" fairly agree with those calculated from MDS, due to the used greater but closer im D value 0.355. It was checked that for this case, in "FFAM 2", the channel height is increased only by 8%.
For the case h = 2.898 nm. Figure 3c shows the comparison between FFAM and MDS for the case h = 2.898 nm and χ = 1.0 . It is shown that FFAM fairly agrees with MDS. In FFAM, the value of im D is 0.5476 for giving the channel height 2.898 nm; this im D value is relatively close to that (0.7086) calculated from MDS. where q m is the mass flow rate per unit channel length, u a and u b are respectively the velocities of the fluid molecules on the upper and lower channel walls, h is the channel height, dp/dx is the driving pressure gradient ( p is the pressure and x is the coordinate), ρ eff bf and η eff bf are respectively the average density and the effective viscosity of the fluid across the channel height, and S is the parameter characterizing the non-continuum effect across the channel height.
In the application of FFAM 56,57 , we normally do not need to obtain the input parameter values calibrated from MDS as shown in the above section. Mostly by other way, as shown in Eq. (13), ρ eff bf and η eff bf can respectively be formulated as the functions of the channel height h by the regression equations based on the experimentally measured, MDS or FFAM calculated values of them, S can also be regressed out as the function of h based on the MDS or FFAM calculation results, and u a and u b can respectively be solved based on the no-slip or slip boundary conditions. When the formulations for all these parameters are given, Eq. (13) can be applied for solving the nanoscale flow problem in engineering such as in nanoporous filtration membranes, cellular connexon and membranes and micro/nano bearings etc., which are not solvable by MDS due to the limitation by the computational capacity. The advantage of Eq. (13) (or FFAM) is its convenience for calculating the nanoscale flow problem, just like we solve the flow problem by using the classical continuum hydrodynamic equation. These are obviously substantially progressive towards solving the engineering nanoscale flow problems, not affordable by MDS. www.nature.com/scientificreports/

Conclusions
The flow velocities in the nano slit pores in the Poiseuille flow were calculated from the flow factor approach model (FFAM) by using all the input parameter values obtained from molecular dynamics simulation (MDS). As a new progressive comparison work, in the present study the physical nature of the simulated system including the averaged local density and local viscosity profiles across the channel height in FFAM was forehand set as the same with that in MDS. It was found that FFAM well agrees with MDS because of the very close flow velocities calculated from these two approaches for the same cases. Since FFAM is the equivalent model for a nanochannel flow and the fluid molecules across the channel height in an actual case is often not exactly ordered normal to the channel wall surface, in FFAM the value of im D should often be required to be corrected to give the correct channel height; it can be a little smaller than that calculated from MDS; this im D value should be understood as the equivalent value. As a substantiated model, FFAM can be successfully used in the modeling of nanoscale flows such as occurring in nanoporous filtration membranes 56 , cellular connexons and membranes and micro/nano bearings 57 etc. As a constitutive model, FFAM radically alters the method of the modeling of the molecular scale flow compared to MDS. The advantage of FFAM is its high efficiency in the engineering nanoscale flow modeling just with the modest cost of computational time and computer storage 56,57 , not owned by MDS [9][10][11][12][13][14][15] . FFAM should become a reliable analytical tool for studying the relevant subjects in nanofluidics. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.