Next Article in Journal
Conformations of 1-Butyl-3-methylimidazolium Chloride Probed by High Pressure Raman Spectroscopy
Previous Article in Journal
Prediction of Environmental Properties for Chlorophenols with Posetic Quantitative Super-Structure/Property Relationships (QSSPR)
Article

A Velocity Extraction Method in Molecular Dynamic Simulation of Low Speed Nanoscale Flows

School of Civil Engineering and Mechanics, Yanshan University, Qinhuangdao, Hebei Province, P.R. China
Int. J. Mol. Sci. 2006, 7(9), 405-416; https://0-doi-org.brum.beds.ac.uk/10.3390/i7090405
Received: 31 May 2006 / Revised: 14 September 2006 / Accepted: 19 September 2006 / Published: 28 September 2006

Abstract

A new algorithm to extract the velocity caused by the external forces inmolecular dynamic simulation of nanoscale flow problems is proposed. The flowvelocity, an important component in these type of problems, is usually obtained from theaverage value in the time space because the accumulation of the thermal velocity willapproach zero when the time period is large, but this method is not always suitable,especially when the flow velocity is much smaller than the thermal velocity. Based on theidea of the linear accumulation of the flow velocity, in this study a new algorithm isderived to extract the flow velocity. This algorithm can be used to calculate nanoscaleflow problem no matter whether the value of the flow velocity is big or small. Using thisnew algorithm, the 2-D liquid flow of argon in nanochannels was simulated. Thenumerical result demonstrates the effectiveness of the new algorithm.
Keywords: Velocity extraction; molecular dynamic simulation; low speed; nanoscale flow. Velocity extraction; molecular dynamic simulation; low speed; nanoscale flow.

Introduction

Nanoscale flow studies are important in physiology, medicine and for design of nanoscale devices, etc. Molecular dynamic simulation (MDS) is used to calculate nanoscale flow problems because the continuum model is no longer valid when the dimensions of flow systems approach the molecular size. Good reviews regarding the models and solution methods of nanoscale flow problems have been published [1,2]. In many papers, i.e. by Travis et al. [3,4], Thompson and Robbins [5], Heinbuch and Fischer [6], MDS methods were applied to resolve nanoscale flow problems, particularly for the study of the boundary slip phenomena. 2-D and 3-D pressure driven channel flow problems were also investigated using the MDS method. In most papers the fluid was liquid argon and the Lennard-Jones potential was used to describe the intermolecular interactions. Usually the Couette or Poiseuille steady flow was considered to illustrate the effectiveness of the MDS method. Travis et al. [3,4] found that if the pore width was less than five molecular diameters, the Navier-Stokes (NS) equations break down. Travis et al. also reported [3] that when the channel width equals 5.1 molecular diameters, the solution given by the molecular dynamics approach deviates significantly from the NS equations solution, whereas when the channel width equals 10.2 molecular diameters, the molecular dynamics’ solution is similar to that of the Navier-Stokes equations. This means that when the channel width is very small the NS solution is meaningless and we should use molecular dynamics to simulate it. Zhang et al. [7] calculated a confined fluid with a nonequilibrium MDS method and discussed the phenomena of constant pore size and constant load. Thompson and Robbins [5] simulated boundary slips with different solid wall density and interaction strength parameters. They found that the boundary slip could be observed only for very weak fluid-wall interactions. Similar problems were also studied by Koplik et al. [8] and Koplik and Banavar [9]. Todd et al. [10] introduced an effective viscosity method to deal with the inhomogeneous viscosity distribution problem. Instead of solving the difficult inhomogeneous viscosity problem, the effective viscosity method adjusts the real viscosity to match the molecular dynamic simulation result. Moseler and Landman [11] and Eggers [12] have discussed nano jet problems.
Although the MDS method has been used to model some nanoscale flow problems, as mentioned above, there is a crucial problem associated with this approach. In the MDS method, the flow velocity of the molecules caused by the applied external force is mixed with the thermal velocity. If the flow velocity is much smaller than the thermal velocity, which is usually more than 100 m·s-1, it is difficult to extract the flow velocity from the much bigger thermal velocity. As most nanoscale flow problems relevant to real applications deal with low speed flows, the direct application of MDS method is not workable. Until now, in all studies the minimum flow velocity used was several meters per second.
Certainly, MDS is correct for nanoscale flow. The problem is the method used to calculate the flow velocity is not suitable. Thus, how to extract the flow velocity is an important part in the MDS of nanoscale flow. The traditional way to extract the flow velocity is using the idea of time average because the accumulation of the thermal velocity will approach to zero while the time period is large, but this method is not always appropriate, especially when the flow velocity is less than several meters per second. That is, numerically, the accumulation of the thermal velocity will never reach zero and the algorithm based on this idea is not valid for low speed problems, regardless of the elegance of the algorithm.
In this paper we propose a new algorithm to extract the velocity caused by the external forces. The idea of this algorithm is as follows: in each time step, we calculate an increment of the flow velocity, and the total flow velocity is the accumulation of the increments in all time steps. This idea is reasonable because at any given time step, we can get the external force increment from the linear development of the total force function. This algorithm can be used to calculate nanoscale flow problem even the value of the flow velocity is very small. Furthermore, when using MDS method in the surface problems [13,14,15,16], fracture mechanics [17,18,19], heat transfer problems, lubrication problems [20], material science [21] and friction problems [22], we also need to calculate the moving velocity caused by the external force, or other external conditions, so the algorithm presented here is also suitable for these kinds of problems. To demonstrate the effectiveness of this new method, we simulated several examples of 2-D flow in slit nano-channels with different channel heights.

Principle of the new method

Let us consider an isothermal pressure-driven steady state 2-D Poiseuille flow in a straight channel. Let x represent the direction of flow, while h represents the channel height in the y direction, that is the direction perpendicular to the flow. Each channel wall, i.e., the boundary perpendicular to the y direction, is modeled in terms of four rows of particles of solid wall material [5]. Initially, both the fluid particles and the solid wall particles are positioned in fcc lattices [23]. In the x direction we use the periodic boundary conditions [24].
If v represents the total velocity vector of the fluid, it can be expressed by
v = vT + ve .
In a macrosopoic problem, if there is no external force, the fluid will be in an equilibrium state, and the flow velocity v = 0, but when we simulate a nanoscale problem with the MDS method, even if there is no external force and it is an equilibrium molecular dynamics problem, the calculated flow velocity will never be zero. Certainly, at this time, the nonzero flow velocity is due to the modeling and numerical errors. We represent this part of the velocity by vT. When external forces are exerted, a real flow velocity of the fluid will exist. We represent this part of velocity by ve. Theoretically, vT should be zero, but numerically it is not zero. In fact, in MDS vT is not only nonzero, but it is quite big and usually equals several hundred meters per second, so if the real flow velocity is small, it is difficult to know the exact value of ve according to equation (1).
In MDS methods [3,4,5], the flow velocity ve can be calculated by
v e = lim N t 1 N t   k = 1 N t v k = lim N t 1 N t   k = 1 N t ( v k T + v k e ) ,
where N t is the number of time steps. The subscript k refers to the kth time step. Equation (2) is based on the concept that if N t approaches infinity
v T = lim N t 1 N t   k = 1 N t v k T = 0 .
This is the traditional method used in most MDS algorithms. During the calculation we consider a small volume with a specific point inside and vk is the average velocity vector of all particles inside this volume at a particular moment [24]. Many studies verifying this method are available when the scalar velocity ve is not less than several meters per second, but when ve is much smaller than the scalar velocity vT, the above method does not work because, numerically, vT = 0 is not true even when N t is sufficiently big. Therefore we have to devise a new way to calculate the velocity ve.
If , at the k-1th time step, we know the velocity of the ith particle vik-1, then, at the kth time step, the velocity of the ith particle vik can be expressed by
v ik = v ik 1 T + v ik 1 e + Δ v ik T + Δ v ik e .
Let F ik represent the total force exerted on the ith particle at the end of the k-1th time step; F ik T and F ik e represent the force corresponding the equilibrium MDS problem and the force caused by the external loads, respectively. If there is no external load then F ik e = 0 , but F ik T is still nonzero. We have F ik = F ik T + F ik e . In equation (4) the increment of velocity Δ v ik T + Δ v ik e depends on the force F ik . Here Δ v i k T is the velocity increment caused by F ik T at the kth time step; Δ v ik e is the velocity increment caused by F ik e at the kth time step. So the flow velocity of the ith particle at the kth time step v ik e can be calculated by
v ik e = v ik 1 e + Δ v ik e = v ik 1 e + F ik e m i Δ t ,
and we have the flow velocity
v ik e = k = 1 N t Δ v ik e = k = 1 N t F ik e m i Δ t ,
where the initial value of the flow velocity v i 0 e = 0 and m i is the mass of the ith particle.
Figure 1 illustrates the positions of particles i and j at moment t and t + Δ t , corresponding to the k-1th and the kth time steps, respectively. The vector distance between these two particles at moment t is r i j . The vector velocity difference between these two particles is Δ v ij . Because only the relative positions is involved, without lose of generality, we can assume that the position of particle i is not changed. At moment t + Δ t , because of the velocity difference Δ v ij , the particle j moves from point j to point j2. The distance between these two particles is r ij + Δ r ij , where Δ r ij = Δ v ij Δ t . The distance change is Δ r ij = Δ r ij T + Δ r ij e , where Δ r ij T = Δ v ij T Δ t ; Δ r ij e = Δ v ij e Δ t ; Δ v ij T and Δ v ij e are the velocity increments caused by the heat movement and the external forces, respectively.
From Figure 1 we see that Δ r ij = Δ r ij T + Δ r ij e = Δ v ij T Δ t + Δ v ij e Δ t . It means the distance between particles i and j is influenced by both Δ v ij T and Δ v ij e . Obviously, the distances between particles i and j are different for Δ v ij e = 0 and Δ v ij e 0 . Also, the interaction forces are different for Δ v ij e = 0 and Δ v ij e 0 because they are functions of r ij + Δ r i j . At moment t + Δ t , the interaction force is f ( r ij + Δ r ij ) = f ( r ij + Δ r ij T + Δ r ij e ) . The form of f ( r ij + Δ r ij ) depends on the potential function. Thus, we think the interaction force as consisting of two parts. The first part corresponds to Δ v ij e = 0 , that is f ( r ij + Δ r ij T ) , we represent it by F ij T . The second part is the additional force increment because Δ v ij e 0 , that is f ( r ij + Δ r ij T + Δ r ij e ) f ( r ij + Δ r ij T ) , we represent it by F ij eI . If there is no external force, we have Δ v ij e = 0 and then F ij eI = 0 . So F ij eI is caused by external force. Now F ik e consists of two parts, F ik e = F ik eI + F ik eII , where F ik eI represents the interaction force due to the external load caused additional distance change between different particles; F i k e I I represents the external load exerted directly on particle i, usually m i g i . F i k e I I depends on the form of the external load. If the potential function between the ith particle and the jth particle is u ( r i j ) we have the following interaction force function
F ij T + F ij eI = u ( r i j ) = u ( r i j ) r i j r i j r i j .
The interaction force increment between t and t+Δ t is
Δ F ij T + Δ F ij eI d ( F ij T + F ij eI ) d t Δ t = d d t ( u ( r i j ) r i j r i j r i j ) Δ t = 2 u ( r i j ) ( r i j ) 2 d r i j d t Δ t r i j r i j u ( r i j ) r i j d d t ( r i j r i j ) Δ t = 2 u ( r i j ) ( r i j ) 2 Δ r i j r i j r i j u ( r i j ) r i j 1 r i j 2 ( d r i j d t r i j d r i j d t r i j ) Δ t = 2 u ( r i j ) ( r i j ) 2 Δ r i j r i j r i j u ( r i j ) r i j 1 r i j ( Δ r i j Δ r i j r i j r i j ) = 2 u ( r i j ) ( r i j ) 2 ( Δ r i j e + Δ r i j T ) r i j r i j u ( r i j ) r i j 1 r i j [ ( Δ r i j e + Δ r i j T ) ( Δ r i j e + Δ r i j T ) r i j r i j ] ,
where Δ r ij e and Δ r ij e are the change of the scalar and vector distance r ij and r ij caused by the external conditions, respectively; Δ r ij T and Δ r ij T are the change of the scalar and vector distance r ij and r ij without external loads, respectively. From (8) we know that the increment of the interaction force between the ith particle and the jth particle caused by the external conditions is
Δ F ij eI = 2 u ( r i j ) ( r i j ) 2 Δ r i j e r i j r i j u ( r i j ) r i j 1 r i j ( Δ r i j e Δ r i j e r i j r i j ) .
Thus, at the kth time step, the total interaction force of the ith particle, caused by the external conditions, can be expressed by
F ik eI = j F ijk eI = F ik - 1 eI + j Δ F ijk - 1 eI = F ik - 1 eI j [ 2 u ( r i j ) ( r i j ) 2 Δ r i j k 1 e r i j k 1 r i j k 1 + u ( r i j ) r i j 1 r i j k 1 ( Δ r i j k 1 e Δ r i j e r i j k 1 r i j k 1 ) ] ,
Δ r ijk - 1 e and Δ r ijk - 1 e are the change of the scalar distance r ij and the vector distance r ij at the k-1th time step caused by the external conditions. Δ r ijk - 1 e and Δ r ijk - 1 e can be calculated according to the positions of particle i and j. At the beginning of simulation, fluid particles are allowed to move without applying the external force until a thermodynamic equilibrium state is reached. At this stage, it is an equilibrium problem. Usually, this process needs thousands of time steps or more. Then the external forces are applied and the non-equilibrium simulation starts.
According to the geometrical relation in Figure 1 we have the following scalar distance expressions:
( r ij + Δ r ij ) 2 = r ij 2 + ( Δ v ij Δ t ) 2 + 2 Δ t Δ v ij . r ij ,
and
( r ij + Δ r ij T ) 2 = r ij 2 + ( Δ v ij T Δ t ) 2 + 2 Δ t Δ v ij T . r ij .
Subtracting (12) from (11) and considering Δ r ij e = Δ r ij Δ r ij T , after some mathematical manipulations, we obtain:
Δ r ij e ( 2 r ij + Δ r ij + Δ r ij T ) = ( r ij + Δ r ij ) 2 ( r ij + Δ r ij T ) 2 = ( Δ v ij Δ t ) 2 ( Δ v ij T Δ t ) 2 + 2 Δ t Δ v ij . r ij 2 Δ t Δ v ij T . r ij = ( Δ v ij 2 ( Δv ij T ) 2 ) Δ t 2 + 2 Δ t Δ v ij e . r ij = Δ v ij e . ( 2 Δ v ij Δ v ij e ) Δ t 2 + 2 Δ t Δ v ij e . r ij .
On the left hand side of (13) both Δ r ij e and r ij are scalar lengths. Usually Δ t is small and both Δ r ij and Δ r ij T are much smaller than r ij . We should notice here the smaller Δ r ij and Δ r ij T are ensured by the size of Δ t , so we do not need to assume the flow velocity is small. Considering that all the variables in (10) correspond to the k-1th time step, we add a subscript k-1 to represent the value at this time step and then equation (13) can be simplified to
Δ r ijk - 1 e = ( Δ v ijk - 1 e . ( 2 Δ v ijk - 1 Δ v ijk - 1 e ) Δ t 2 + 2 Δ t Δ v ijk - 1 e . r ijk - 1 ) / ( 2 r ijk - 1 ) .
The procedure of the algorithm is: at the kth time step, for any particle i, we calculate Δ r ijk - 1 e by formulae (14) and Δ r ijk - 1 e = Δ v ijk - 1 e Δ t , then the force exerted on this particle at this moment F ijk e and F ik e by (10). Finally, the flow velocity increment Δ v ik e and the flow velocity v ik e can be calculated by (5) and (6), respectively. When using (14) to calculate Δ r ijk - 1 e the term Δ v ijk - 1 e is the velocity increment corresponding to the k-1th time step and it is already known.
The considered Poiseuille flow is a steady state problem. The system temperature should be fixed throughout the flow process. To avoid increases of the system temperature during the simulation due to external forces or other unphysical effects, a constant temperature constraint is added. If the external forces are included in the system, it becomes a non-equilibrium molecular dynamics problem. Reference [25] is a good review regarding nonequilibrium MDS in flow problems. In this case the movement equation of particle i is
m i v ˙ i = F i + m i g i α ( v i v i e ) .
This is so called Gaussian thermostat method for nonequilibrium flow problems [25,26]. For equilibrium flow problems, both g i and v i e are zero vectors in (15). Because the flow velocity should not affect the temperature of the system, in the last term of (15), we deduct it from the total velocity. Here v ˙ i is the derivative vector of the velocity of the ith particle with respect to the time variable. m i g i is the external force vector F i eII . F i is the total interaction force exerted on the ith particle. It is composed of two parts, F i = F i T + F i eI . F i T is the interaction force without external loads. F i eI is the additional interaction force caused by the external force. We do not need to calculate F i separately. This means that during the simulation, we do not need to use the sum of F i T and F i eI to get F i . We can calculate F i according to F i = j i ( F ij T + F ij eI ) = j i u ( r i j ) and the position of particles directly. α is a thermostat multiplier to keep the kinetic temperature constant in the system. It can be obtained from the constant temperature condition [25,26]. This is one of the usually used methods in thermostat molecular dynamics. When calculate the system temperature, for each particle, the flow velocity should be subtracted from the total velocity.

Simulations of 2-D Poiseuille flow

In this Poiseuille flow we assume that the external force exerted on fluid particle i is m i g i , that is F ik eII = m i g i . Here gi is the equivalent constant gravity acceleration vector. There are N particles within a length L along the flow direction and a height h of the flow channel. The initial velocity is given by the original temperature. In these simulations the fluid material is liquid argon and the Lennard-Jones potential is used. The cutoff radius is taken as rc = 2.5σ. The interaction between the liquid and the solid wall is considered by the interaction strength parameter of the wall particles to the liquid particles. As mentioned before, in our simulations, the solid wall particles are fixed. Certainly this fixed boundary particle assumption will influence the accuracy of the flow velocity a little but it is not affect the examination of the effectiveness of our new velocity extraction method. During the simulation we fix the positions of these four rows of particles on each side.
The values of the parameters used in the calculations are: the size of particles σf = σw = 3.4x10-10 m; the mass of the particles mf = mw = 6.69x10-26 kg; the temperature T = 84K; the interaction strength parameter between argon molecules is εf = 1.657x106 kg·m-1·s-1. In the above, the subscript f represents fluid and the subscript w represents solid wall. In our simulations, the size of each time step is chosen as Δt = 10-16 s. In many papers Δt = 10-14 s is recommended, but to ensure our linearized algorithm is stable Δt = 10-16 s is necessary. In principle, we can use different interaction strength parameters εw to represent different solid wall materials [24,27]. We used εw = 3εf in our simulations. The interaction strength parameter for the interaction between a wall particle and a fluid particle εwf can be evaluated [24,27] by εwf = (εfεw)½. In all the simulations the density number is 0.82, corresponding to the liquid density of argon at T=84 K.
We calculated two different cases with different channel heights. For the first case, L=13.61nm; h=2.31 nm; N=243. For the second case, L=13.61 nm; h= 4.36 nm; N=459. The velocity profiles calculated in our simulations for these two cases are plotted in Figure 2, Figure 3, Figure 4 and Figure 5.
Because in this new method the total flow velocity ve is the linear accumulation of velocity increment of each time step we call it the LA method. The traditional method is based on the idea of time average, we call this the TA method. Figure 2 shows the simulation of case 1. In this simulation the equivalent acceleration in the flow direction is 1013 m·s-2 and the simulated maximum velocity which is calculated in the x direction by the method of this paper (the LA method) is 39.9 m·s-1. The MDS result where the velocity is calculated by the traditional method (TA method) is also shown in the Figure. Because the velocity is big both the TA and LA methods can be used to calculate the flow velocity. Figure 3 is another example of case 1, where the pressure difference is smaller. The equivalent acceleration in the flow direction is 109 m·s-2 and the maximum velocity in the x direction is 3.92x10-3 m·s-1.Because the flow velocity is much smaller than the thermal velocity, the traditional method (TA method) is not available, but, as shown in the Figure, the LA method proposed in this paper is still effective. Figure 4 shows the simulation of case 2. In this simulation the equivalent acceleration in the flow direction is 1012 m·s-2 and the maximum velocity in the x direction is 13.92 m·s‑1. Figure 5 is another example of case 2, where the equivalent acceleration in the flow direction is 108 m·s-2 and the maximum velocity in the x direction is 1.37x10-3 m·s-1. In Figure 3 and Figure 5 we can only get results for the LA method because in these simulation situations or for such small flow velocities the TA method is not available. In fact, when we simulated the same problems with the TA method under the same conditions as in Figure 3 and Figure 5, the maximal flow velocity is always great than 2 m·s-1. Certainly this is not the correct solution of the problems. Up to now, we have not found any algorithm which can be used to solve such a low velocity nanoscale flow problem. That is why the LA method is important for low flow velocity problems. However, we can compare the results of TA and LA methods for high flow velocity cases, as shown in Figure 2 and Figure 4. Because, in principle, the LA method can be used for both high and low flow velocity problems the conclusions obtained from the high flow velocity comparisons can be extended to low flow velocity problems.
These numerical examples demonstrate that the velocity extraction method proposed in this paper is a suitable approach for both low and high speed flow problems. When the flow velocity is big the result of LA method is the same as the result of TA method. When the flow velocity is small TA method is not available we have to use LA method to extract ve.

Convergence and Efficiency analysis

In this new algorithm, we only used the linear terms of the development of the force function with respect to Δ r ij e to approximate the force increment caused by the external conditions. What is the influence of the truncate errors of this approximation? Is it convergent and what are the conditions of convergence? These are the questions which should be answered for a reliable algorithm. We do a series of simulations with different size of time step from Δt = 10-14 s to Δt = 10-18 s. We find that when Δt = 10-14 s the simulation of our new velocity extraction method is not good. This means for the linear approximation the size of time step should be smaller. If the time step is big the small Δ r ij assumption will be broken. When Δt = 10-18 s the result is almost the same with the result of Δt = 10-17 s. So we have reason to believe that when Δt = 10-17 s the simulation is already convergent. We show the simulation results with time step size from Δt = 10-15 s to Δt = 10-17 s in Figure 6. From this figure we see that when Δt = 10-15 s the simulation result is far from convergence and when Δt = 10-16 s we obtain a quite satisfactory result. If we let Δt = 10-17 s the result is only improved a little, so for practical applications, the recommended time step size is Δt = 10-16 s.
(LA Method). Low flow velocity simulations of case 2. The vertical coordinate is the flow velocity and the horizontal coordinate is the height of the channel.
Figure 7 is the analysis of the error of the maximal flow velocity when the size of the time step is changed (LA Method). The simulation conditions are the same as the simulations in Figure 4. The vertical coordinate is the error of the maximal flow velocity | ( v maxLA v maxTA ) / v maxTA | 100% and the horizontal coordinate is size of the time step Δ t (s), where vmax LA is the maximal velocity of this paper and vmax TA is the maximal velocity of TA method.
Although the purpose of this new method is to extract the low flow velocity caused by the external conditions, which the conventional MDS method can not do, the efficiency of an algorithm is always an important part to be considered. In our new method we use Δt = 10-16 s. For the similar problems, in conventional MDS methods, Δt = 10-14 s is usually used. Because in this new method we do not need extra iterations to calculate the flow velocity, even 100 times smaller time step size is used, the total number of time steps is still less than that of the conventional MDS method. The simulation software was executed on a Pentium 2.8GHz PC. We compared the CPU time and the number of time steps as follows. In Figure 2, we used 89800 time steps (51 minutes CPU time) for LA method and 3.2x105 time steps (108 minutes CPU time) for the TA method; in Figure 4, we used 288400 time steps (565 minutes CPU time) for the LA method and 7.8 x105 time steps (1022 minutes CPU time) for the TA method. The efficiency of the new method is evident.

Conclusions

Using this method we can separate very small ve from the total velocity. It is very useful because in many cases we are more interested in the real movement of the medium. The new algorithm is also efficient compared to the traditional way. Because we do not need so many extra iterations as we usually do the new algorithm can reach a stable solution faster. This is important because usually the MDS use so much computation time that the number of particles is limited.
In this paper only the steady Poiseuille flow is simulated. Because, in this method, we do not use the time average algorithm the velocity computation is totally coincidence with the real transient state so the dynamic effective can be exhibited exactly. Hopefully this method will be suitable to simulate transient phenomena. It will be our later work. This method is also valid in other MDS problems where ve is important. In fact, in many MDS problems we want to know the small change of velocity caused by the exerted external conditions. The method will provide us a good way to solve this kind of problems.

Acknowledgements

This work was supported by grants from Yanshan University, China.

References

  1. Polubotko, A.M.; Hatta, A.; Suzuki, Y.; Ciofalo, M.; Collins, M.W.; Hennessy, T.R. Modelling nanoscale fluid dynamics and transport in physiological flows. Med. Eng. Phys. 1996, 18, 437–451. [Google Scholar] [CrossRef]
  2. Giordano, N.; Cheng, J-T. Microfluid mechanics: progress and opportunities. J. Phys. Condens. Mat. 2001, 13, 271–295. [Google Scholar] [CrossRef]
  3. Travis, K.P.; Todd, B.D.; Denis, J.E. Departure from Navier-Stokes hydrodynamics in confined liquids. Phys. Rev. 1997, 55, 4288–4295. [Google Scholar]
  4. Travis, K. P.; Gubbins, K. E. Poiseuille flow of Lennard-Jones fluids in narrow slit. J. Chem. Phys. 2000, 112, 1984–1994. [Google Scholar] [CrossRef]
  5. Thompson, P.A.; Robbins, M.O. Shear flow near solids: Epitaxial order and flow boundary conditions. Phys. Rev. A 1990, 41, 6830–6839. [Google Scholar] [CrossRef]
  6. Heinbuch, U.; Fischer, J. Liquid flow in pores: slip, no-slip or multiplayer sticking. Phys. Rev. A 1989, 40, 1144–1146. [Google Scholar] [CrossRef]
  7. Zhang, L.; Balasundaram, R.; Gehrke, S.H.; Jiang, S. Nonequilibrium molecular dynamics simulations of confined fluids in contact with the bulk. J. Chem. Phys. 2001, 114, 6869–6877. [Google Scholar] [CrossRef]
  8. Koplik, J.; Banavar, J.R.; Willemsen, J.F. Molecular Dynamics of Poiseuile flow and moving contact line. Phys. Rev. Lett. 1988, 60, 1282–1285. [Google Scholar] [CrossRef]
  9. Koplik, J.; Banavar, J.R. Corner flow in the sliding problem. Phys. Fluids 1995, 7, 3118–3125. [Google Scholar] [CrossRef]
  10. Todd, B.D.; Evans, D.J.; Daivis, P.J. Pressure tensor for inhomogeneous fluids. Phys. Rev. E 1995, 52, 1627–1638. [Google Scholar] [CrossRef]
  11. Moseler, M.; Landman, U. Formation, Stability and Breakup of Nanojets. Science 2000, 289, 1165–1169. [Google Scholar] [CrossRef]
  12. Eggers, J. Dynamics of liquid nanojets. Phys. Rev. Lett. 2002, 89, 084502. [Google Scholar] [CrossRef]
  13. Yang, J.X.; Koplik, J.; Banavar, J. R. Molecular Dynamics of Drop Spreading on a Solid Surface. Phys. Rev. Lett. 1991, 67, 3539–3542. [Google Scholar] [CrossRef]
  14. Klien, J.; Perahia, D.; Warburg, S. Forces between polymer-bearing surfaces undergoing shear. Nature 1991, 352, 143–145. [Google Scholar] [CrossRef]
  15. Fang, T.H.; Weng, C.I.; Chang, J.G. Molecular dynamics simulation of nano-lithography process using atomic force microscopy. Surf. Sci. 2002, 501, 138–147. [Google Scholar] [CrossRef]
  16. Fukui, K.; Sumpter, B.G.; Runge, K.; Kung, C.Y.; Barnes, M.D.; Noid, D.W. Collision dynamics and surface wetting of nano-scale polymer particles on substrates. Chem. Phys. 1999, 244, 339–349. [Google Scholar] [CrossRef]
  17. Landman, U.; Luedtke, W.D.; Burnham, N.A.; Colton, R.J. Microscopic Mechanisms and Dynamics of Adhesion, Microindentation and Fracture. Science 1990, 248, 454–461. [Google Scholar]
  18. Kohlhoff, S.; Gumbsch, P.; Fischmeister, H.F. Crack propagation in b.c.c. crystals studied with a combined finite-element and atomistic model. Philos. Mag. A 1991, 64, 851–878. [Google Scholar] [CrossRef]
  19. Inoue, H.; Akahoshi, Y.; Harada, S. Molecular dynamics simulation on fracture mechanisms of nanoscale polycrystal under static and cyclic loading. Int. J. Fatigue 1997, 19, 266–266. [Google Scholar]
  20. Sham, T.L.; Tichy, J. Scheme for hybrid molecular dynamics/finite element analysis of thin film lubrication. Wear 1997, 207, 100–106. [Google Scholar] [CrossRef]
  21. Wang, Y.U.; Jin, Y.M.; Cuitino, A.M.; Khachaturyan, A.G. Nanoscale phase field microelasticity theory of dislocations: model and 3D simulations. Acta Mater. 2001, 49, 1847–1857. [Google Scholar] [CrossRef]
  22. Matthey, T.; Hansen, J.P. Molecular dynamics simulation of sliding friction in a dense granular material. Model. Simul. Mater. Sci. Eng. 1998, 6, 701–707. [Google Scholar] [CrossRef]
  23. Somers, S. A.; Davis, H. T. Microscopic dynamics of fluids confined between smooth and atomically structured solid surfaces. J. Chem. Phys. 1992, 96, 5389–5407. [Google Scholar] [CrossRef]
  24. Rapaport, D.C. The art of molecular dynamics simulation; Cambridge University Press: Cambridge, 1995; p. 16. [Google Scholar]
  25. Todd, B.D. Computer simulation of simple and complex atomistic fluids by nonequilibrium molecular dynamics techniques. Comput. Phys. Commun. 2001, 142, 14–21. [Google Scholar] [CrossRef]
  26. Evans, D.J.; Morris, G.P. Statistical mechanics of nonequilibrium liquids; Academic Press: New York, London, 1990. [Google Scholar]
  27. Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular forces; Oxford University Press: New York, 1981; p. 137. [Google Scholar]
Figure 1. Distance between particles i and j at moment t and t + Δt.
Figure 1. Distance between particles i and j at moment t and t + Δt.
Ijms 07 00405 g001
Figure 2. Poiseuille flow simulations of case 1 (high flow velocity simulations). The vertical coordinate is the flow velocity and the horizontal coordinate is the height of the channel. Ijms 07 00405 i001 is the result of this paper (LA ). Ijms 07 00405 i002 is the result of the traditional (TA) method.
Figure 2. Poiseuille flow simulations of case 1 (high flow velocity simulations). The vertical coordinate is the flow velocity and the horizontal coordinate is the height of the channel. Ijms 07 00405 i001 is the result of this paper (LA ). Ijms 07 00405 i002 is the result of the traditional (TA) method.
Ijms 07 00405 g002
Figure 3. Poiseuille flow simulations of case 1 (low flow velocity simulations). Ijms 07 00405 i001 is the result of this paper (LA ).
Figure 3. Poiseuille flow simulations of case 1 (low flow velocity simulations). Ijms 07 00405 i001 is the result of this paper (LA ).
Ijms 07 00405 g003
Figure 4. Poiseuille flow simulations of case 2 (high flow velocity simulations). Ijms 07 00405 i001 is the result of this paper (LA ). Ijms 07 00405 i002 is the result of the traditional (TA) method.
Figure 4. Poiseuille flow simulations of case 2 (high flow velocity simulations). Ijms 07 00405 i001 is the result of this paper (LA ). Ijms 07 00405 i002 is the result of the traditional (TA) method.
Ijms 07 00405 g004
Figure 5. Poiseuille flow simulations of case 2 (low flow velocity simulations). Ijms 07 00405 i001 is the result of this paper (LA ).
Figure 5. Poiseuille flow simulations of case 2 (low flow velocity simulations). Ijms 07 00405 i001 is the result of this paper (LA ).
Ijms 07 00405 g005
Figure 6. Comparison of simulation results with different Δt (LA Method). Low flow velocity simulations of case 2. The vertical coordinate is the flow velocity and the horizontal coordinate is the height of the channel.
Figure 6. Comparison of simulation results with different Δt (LA Method). Low flow velocity simulations of case 2. The vertical coordinate is the flow velocity and the horizontal coordinate is the height of the channel.
Ijms 07 00405 g006
Figure 7. Illustration of the relationship between the error of the maximal flow velocity and the size of the time step (LA Method).
Figure 7. Illustration of the relationship between the error of the maximal flow velocity and the size of the time step (LA Method).
Ijms 07 00405 g007
Back to TopTop