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 [Formula: see text], FFAM fairly agrees with MDS for all the cases. The study values FFAM in simulating nanochannel flows.

Introduction

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-8]. The studies on them have been mainly experimental or by molecular dynamics simulation (MDS)[1-15]. In these devices, the nanochannel flow has been shown by MDS to be distinctly different from the classical recognition[9-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-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-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-15,22-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 et al.[38] proposed a phenomenological continuum model for force-driven nanochannel flows; they used the effective channel height and the density deficit parameter based on the density layering near the channel walls; their model needs the calibration of the parameters by MDS. Kasiteropoulou et al.[39] used the dissipative particle dynamics method to study the flow properties in a nanochannel. Swift et al.[40] used the lattice Boltzman method to simulate the binary fluid system. Naris and Valougeorgis[41] used the kinetic theory to simulate the non-equilibrium flow in a very narrow channel.

The present study focuses on FFAM for nanochannel flows. The previous studies[36,42-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-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 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.

Figure 1

The flow factor approach model (FFAM) for the flow in a nano slit pore[36,42]. 1-upper plane surface, 2-lower plane surface, 3-ordered fluid molecule.

According to the model, in the pressure driven flow shown in Fig. 1, the velocity of the fluid molecule is[21,36,42]:

Here, is the order number of the fluid molecule across the channel height, is the velocity of the (n − 1)th fluid molecule (on the upper solid wall), is the velocity of the 0th fluid molecule (on the lower solid wall), is the driving pressure, is the fluid molecule diameter, is the coordinate shown in Fig. 1, and and 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 and (for i = 0, 1,…, (n − 1)/2–2). Here, and m are respectively constant, and n is an odd number (). The values of and both may be varied across the channel height; the model takes as the average value of (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. , once the values of D, m, n, , , and are known, the flow velocity 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 can be determined beforehand. While, the values of m, , and 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 . 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 , , , and are respectively the length unit, the depth of the interaction energy well, the permittivity of vacuum, point charge and the distance between the atoms and . The subscripts and respectively represent C and H. and denote the interaction between different methane molecules. Moreover, in the potential model, the bond length is , and (). The intermolecular potential parameters are given in Table 1.

Table 1

MOPLS potential parameters[47].

Type	Parameter	C–C	C–H	H–H
MOPLS	ε/k_BK	46.8	17.17	6.30
	σÅ	3.45	3.06	2.67

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:

Figure 2

Schematic of the methane molecules confined by two parallel silicon atom flats.

The parameters and were obtained by the Lorentz–Berthelot mixing rule:where is the interaction strength factor between the wall atom and the methane molecule. In this study, and were respectively employed. The methane coarse-graining potential parameters and were optimized, by using the relative entropy coarse-graining framework[49,50] in the bulk methane ensemble at the initial condition and . The details of the methane coarse-graining process have been shown in Ref.[49]. The silicon atom intermolecular potential parameters are: and [51]. An ABAB stacking configuration was used to describe the silicon atom flat walls; moreover, every atom was conducted by the harmonic potential site as follows:

Here, is the atom position at time , and was obtained by using the second derivative value of the silicon atom intermolecular potential at . To avoid the effect of the compressibility of the methane flowing through the nanochannel, the external driving body force 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 and [53,54] ( is the mass of methane molecule, is the system temperature, and and are the velocities in the y and z directions respectively). In this simulation, the reduced units are and ; The simulation step is (), and the energy cutoff circle radius is . In order to obtain the accurate results statistically, the results for the initial simulation steps were discarded; the relevant molecule information and parameters for FFAM were calculated in the next 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 bins. Moreover, the nanochannel dimensions are and (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 , , and of FFAM, we supposed the methane fluid layer is five () when the system arrived at the equilibrium state for two cases. These parameters were respectively calculated by the following equations:

Here, is the averaged gap between the ()th fluid layer and the lth fluid layer, the subscripts z indicates the position of the fluid molecule in the z direction, (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 are respectively put below these tables. Tables 6, 7, 8 and 9 respectively show the values of the local viscosity ratio calculated from MDS indicating the significant variation of the local viscosity across the channel height owing to the fluid-channel wall interaction[58]; the calculated corresponding values of are respectively put below these tables.

Table 2

The separations between the neighboring fluid molecules across the channel height calculated from MDS,  = 2.484 nm, ,  = 5.

Separation	Δ_0/D	Δ_1/D	Δ_2/D	Δ_3/D
Value	0.43117	0.51157	0.51193	0.43167

 = 1.1862.

Table 3

The separations between the neighboring fluid molecules across the channel height calculated from MDS,  = 2.484 nm, ,  = 5.

Separation	Δ_0/D	Δ_1/D	Δ_2/D	Δ_3/D
Value	0.46511	0.52826	0.52809	0.46512

 = 1.1356.

Table 4

The separations between the neighboring fluid molecules across the channel height calculated from MDS,  = 2.898 nm, ,  = 5.

Separation	Δ_0/D	Δ_1/D	Δ_2/D	Δ_3/D
Value	0.58531	0.70889	0.70838	0.58549

 = 1.2105.

Table 5

The separations between the neighboring fluid molecules across the channel height calculated from MDS,  = 2.898 nm, ,  = 5.

Separation	Δ_0/D	Δ_1/D	Δ_2/D	Δ_3/D
Value	0.60681	0.72269	0.72285	0.6066

 = 1.1913.

Table 6

The values of calculated from MDS,  = 2.484 nm, ,  = 5.

Ratio of local viscosity	η_line,0/η_line,1	η_line,2/η_line,3
Value	1.3913	0.7182 (= 1/1.3924)

Average value of the local viscosity ratio = 1.3918, m = 1.9362.

Table 7

The values of calculated from MDS,  = 2.484 nm, ,  = 5.

Ratio of local viscosity	η_line,0/η_line,1	η_line,2/η_line,3
Value	1.6248	0.6154 (= 1/1.625)

Average value of the local viscosity ratio = 1.6249, m = 3.818.

Table 8

The values of calculated from MDS,  = 2.898 nm, ,  = 5.

Ratio of local viscosity	η_line,0/η_line,1	η_line,2/η_line,3
Value	1.3011	0.7624 (= 1/1.3116)

Average value of the local viscosity ratio = 1.3064, m = 1.3992.

Table 9

The values of calculated from MDS,  = 2.898 nm, ,  = 5.

Ratio of local viscosity	η_line,0/η_line,1	η_line,2/η_line,3
Value	1.7175	0.5822 (= 1/1.7176)

Average value of the local viscosity ratio = 1.7176, m = 3.09.

The values of (im = (n − 1)/2) calculated from MDS are respectively:

for  = 2.484 nm and , ,

for  = 2.484 nm and , ,

for  = 2.898 nm and , ,

for  = 2.898 nm and , .

The values of calculated from MDS are respectively:

for  = 2.484 nm and , ,

for  = 2.484 nm and , ,

for  = 2.898 nm and , ,

for  = 2.898 nm and , .

In the present study, the above values of m, , and calculated by MDS are directly put into Eq. (1) (FFAM); also the values of D, n and 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  = 5, and the pressure gradient for the Poiseuille flow was taken exactly the same as in MDS i.e. .

For the case  = 2.484 nm

Figure 3a shows the comparison between the flow velocities respectively calculated from FFAM and MDS for the case  = 2.484 nm and . In the figure, the legend “FFAM 1” refers to the calculation from FFAM by using . For this case, the value of calculated from MDS is actually 0.5118 as shown above; the reason for using the corrected value in FFAM is to satisfy the channel height 2.484 nm, because 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 , for giving the channel height 2.484 nm, the equivalent number of the fluid molecules across the channel height should be  = 4.11; if we use n = 5, the value of 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 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 to calculate the flow velocities. In Fig. 3a, we also present the calculation results from FFAM by using , 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 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.

Figure 3

Comparisons between the flow factor approach model (FFAM) and molecular dynamics simulation (MDS) in the flow velocity profiles across the channel height. (a)  = 2.484 nm; FFAM 1: , FFAM 2: ; (b)  = 2.484 nm; FFAM 1: , FFAM 2: ; (c)  = 2.898 nm; (d)  = 2.898 nm.

The same circumstances also exist in Fig. 3b, which is for the case  = 2.484 nm and . In Fig. 3b, the results calculated by “FFAM 1” are considerably smaller than those calculated from MDS, due to the used 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 value 0.355. It was checked that for this case, in “FFAM 2”, the channel height is increased only by 8%.

For the case  = 2.898 nm

Figure 3c shows the comparison between FFAM and MDS for the case  = 2.898 nm and . It is shown that FFAM fairly agrees with MDS. In FFAM, the value of is 0.5476 for giving the channel height 2.898 nm; this value is relatively close to that (0.7086) calculated from MDS.

Figure 3d also shows the good matching between FFAM and MDS for the case  = 2.898 nm and . In FFAM, the value of is 0.5436 for giving the channel height; it is relatively close to the value (0.7228) calculated from MDS.

From Fig. 3a–d, it can be found that with the increase of the fluid-channel wall interaction strength i.e. with the increase of the flow velocities across the channel height are significantly reduced; this shows the varying transport properties of the confined fluid influenced by the fluid-channel wall interaction[58]. Such a flow characteristic is actually reflected by the values of m, , and in FFAM as shown above. Different hydrophilicity or hydrophobicity of the channel wall, different values of m, , and in FFAM; if the wall slippage occurs, and/or are not vanishing and they are the wall slipping velocities required to be solved based on the limiting interfacial shear strength model according to FFAM[60].

Commenting remarks

The above comparisons directly validate FFAM. The subsequently made theoretical derivations[21,43] of the flow equation for the nanoscale flow based on FFAM thus firmly stand. The virtues of FFAM is the flow equation it gives[21,43], which is the governing equation for the nanoscale flow and for the two-dimensional flow reads:where is the mass flow rate per unit channel length, and are respectively the velocities of the fluid molecules on the upper and lower channel walls, is the channel height, is the driving pressure gradient ( is the pressure and is the coordinate), and are respectively the average density and the effective viscosity of the fluid across the channel height, and 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), and can respectively be formulated as the functions of the channel height by the regression equations based on the experimentally measured, MDS or FFAM calculated values of them, can also be regressed out as the function of based on the MDS or FFAM calculation results, and and 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.

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 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 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-15]. FFAM should become a reliable analytical tool for studying the relevant subjects in nanofluidics.