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

Abstract

A new algorithm to extract the velocity caused by the external forces in molecular dynamic simulation of nanoscale flow problems is proposed. The flow velocity, an important component in these type of problems, is usually obtained from the average value in the time space because the accumulation of the thermal velocity will approach 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 the idea of the linear accumulation of the flow velocity, in this study a new algorithm is derived to extract the flow velocity. This algorithm can be used to calculate nanoscale flow problem no matter whether the value of the flow velocity is big or small. Using this new algorithm, the 2-D liquid flow of argon in nanochannels was simulated. The numerical result demonstrates the effectiveness of the new algorithm.

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 $\text{v}$ represents the total velocity vector of the fluid, it can be expressed by

v = v^T + v^e .

(1)

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 v^T. When external forces are exerted, a real flow velocity of the fluid will exist. We represent this part of velocity by v^e. Theoretically, v^T should be zero, but numerically it is not zero. In fact, in MDS v^T 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 v^e according to equation (1).

In MDS methods [3,4,5], the flow velocity v^e can be calculated by

$${\text{v}}^{\mathrm{e}}=\underset{{\mathrm{N}}_{\mathrm{t}}\to \infty }{\lim }\frac{1}{{\mathrm{N}}_{\mathrm{t}}}{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{t}}}{\text{v}}_{\mathrm{k}}}=\underset{{\mathrm{N}}_{\mathrm{t}}\to \infty }{\lim }\frac{1}{{\mathrm{N}}_{\mathrm{t}}}{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{t}}}({\text{v}}_{\mathrm{k}}^{\mathrm{T}}+{\text{v}}_{\mathrm{k}}^{\mathrm{e}}),}$$

(2)

where ${\mathrm{N}}_{\mathrm{t}}$ is the number of time steps. The subscript k refers to the k^th time step. Equation (2) is based on the concept that if ${\mathrm{N}}_{\mathrm{t}}$ approaches infinity

$${\text{v}}^{\mathrm{T}}=\underset{{\mathrm{N}}_{\mathrm{t}}\to \infty }{\lim }\frac{1}{{\mathrm{N}}_{\mathrm{t}}}{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{t}}}{\text{v}}_{\mathrm{k}}^{\mathrm{T}}}=\mathbf{0}.$$

(3)

This is the traditional method used in most MDS algorithms. During the calculation we consider a small volume with a specific point inside and v_k 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 v^e is not less than several meters per second, but when v^e is much smaller than the scalar velocity v^T, the above method does not work because, numerically, v^T = 0 is not true even when ${\mathrm{N}}_{\mathrm{t}}$ is sufficiently big. Therefore we have to devise a new way to calculate the velocity v^e.

If , at the k-1^th time step, we know the velocity of the i^th particle v_ik-1, then, at the k^th time step, the velocity of the i^th particle v_ik can be expressed by

$${\text{v}}_{\text{ik}}={\text{v}}_{\text{ik}-1}^{\mathrm{T}}+{\text{v}}_{\text{ik}-1}^{\mathrm{e}}+\Delta {\text{v}}_{\text{ik}}^{\mathrm{T}}+\Delta {\text{v}}_{\text{ik}}^{\mathrm{e}}.$$

(4)

Let $F_{\text{ik}}$ represent the total force exerted on the i^th particle at the end of the k-1^th time step; $F_{\text{ik}}^{\mathrm{T}}$ and $F_{\text{ik}}^{\mathrm{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_{\text{ik}}^{\mathrm{e}}=0$, but $F_{\text{ik}}^{\mathrm{T}}$ is still nonzero. We have $F_{\text{ik}}=F_{\text{ik}}^{\mathrm{T}}+F_{\text{ik}}^{\mathrm{e}}$. In equation (4) the increment of velocity $\Delta {\text{v}}_{\text{ik}}^{\mathrm{T}}+\Delta {\text{v}}_{\text{ik}}^{\mathrm{e}}$ depends on the force $F_{\text{ik}}$. Here $\Delta {\text{v}}_{\mathrm{i}\mathrm{k}}^{\mathrm{T}}$ is the velocity increment caused by $F_{\text{ik}}^{\mathrm{T}}$ at the k^th time step; $\Delta {\text{v}}_{\text{ik}}^{\mathrm{e}}$ is the velocity increment caused by $F_{\text{ik}}^{\mathrm{e}}$ at the k^th time step. So the flow velocity of the i^th particle at the k^th time step ${\text{v}}_{\text{ik}}^{\mathrm{e}}$ can be calculated by

$${\text{v}}_{\text{ik}}^{\mathrm{e}}={\text{v}}_{\text{ik}-1}^{\mathrm{e}}+\Delta {\text{v}}_{\text{ik}}^{\mathrm{e}}={\text{v}}_{\text{ik}-1}^{\mathrm{e}}+\frac{F_{\text{ik}}^{\mathrm{e}}}{{\mathrm{m}}_{\mathrm{i}}}\Delta \mathrm{t}$$

(5)

and we have the flow velocity

$${\text{v}}_{\text{ik}}^{\mathrm{e}}={\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{t}}}\Delta {\text{v}}_{\text{ik}}^{\mathrm{e}}}={\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{t}}}\frac{F_{\text{ik}}^{\mathrm{e}}}{{\mathrm{m}}_{\mathrm{i}}}\Delta \text{t}},$$

(6)

where the initial value of the flow velocity ${\text{v}}_{\text{i}0}^{\mathrm{e}}=0$ and ${\mathrm{m}}_{\mathrm{i}}$ is the mass of the i^th particle.

Figure 1 illustrates the positions of particles i and j at moment t and t + $\Delta \text{t}$, corresponding to the k-1^th and the k^th time steps, respectively. The vector distance between these two particles at moment t is ${\text{r}}_{\mathrm{i}\mathrm{j}}$. The vector velocity difference between these two particles is $\Delta {\text{v}}_{\text{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 + $\Delta \text{t}$, because of the velocity difference $\Delta {\text{v}}_{\text{ij}}$, the particle j moves from point j to point j₂. The distance between these two particles is ${\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}}$, where $\Delta {\text{r}}_{\text{ij}}=\Delta {\text{v}}_{\text{ij}}\Delta \text{t}$. The distance change is $\Delta {\text{r}}_{\text{ij}}=\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{e}}$, where $\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}}=\Delta {\text{v}}_{\text{ij}}^{\mathrm{T}}\Delta \mathrm{t}$; $\Delta {\text{r}}_{\text{ij}}^{\mathrm{e}}=\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}\Delta \mathrm{t}$; $\Delta {\text{v}}_{\text{ij}}^{\mathrm{T}}$ and $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}$ are the velocity increments caused by the heat movement and the external forces, respectively.

From Figure 1 we see that $\Delta {\text{r}}_{\text{ij}}=\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{e}}=\Delta {\text{v}}_{\text{ij}}^{\mathrm{T}}\Delta \mathrm{t}+\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}\Delta \mathrm{t}$. It means the distance between particles i and j is influenced by both $\Delta {\text{v}}_{\text{ij}}^{\mathrm{T}}$ and $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}$. Obviously, the distances between particles i and j are different for $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}=0$ and $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}\ne 0$. Also, the interaction forces are different for $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}=0$ and $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}\ne 0$ because they are functions of ${\text{r}}_{\text{ij}}+\Delta {\text{r}}_{i\mathrm{j}}$. At moment t + $\Delta \text{t}$, the interaction force is $f({\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}})=f({\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{e}})$. The form of $f({\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}})$ depends on the potential function. Thus, we think the interaction force as consisting of two parts. The first part corresponds to $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}=0$, that is $f({\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}})$, we represent it by $F_{\text{ij}}^{\mathrm{T}}$. The second part is the additional force increment because $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}\ne 0$, that is $f({\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{e}})-$ $f({\text{r}}_{\text{ij}}+\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}})$, we represent it by $F_{\text{ij}}^{\text{eI}}$. If there is no external force, we have $\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}=0$ and then $F_{\text{ij}}^{\text{eI}}=0$. So $F_{\text{ij}}^{\text{eI}}$ is caused by external force. Now $F_{\text{ik}}^{\mathrm{e}}$ consists of two parts, $F_{\text{ik}}^{\mathrm{e}}=F_{\text{ik}}^{\text{eI}}+F_{\text{ik}}^{\text{eII}}$, where $F_{\text{ik}}^{\text{eI}}$ represents the interaction force due to the external load caused additional distance change between different particles; $F_{\mathrm{i}\mathrm{k}}^{\mathrm{e}\mathrm{I}\mathrm{I}}$ represents the external load exerted directly on particle i, usually ${\mathrm{m}}_{\mathrm{i}}{g}_{\mathrm{i}}$. $F_{\mathrm{i}\mathrm{k}}^{\mathrm{e}\mathrm{I}\mathrm{I}}$ depends on the form of the external load. If the potential function between the i^th particle and the j^th particle is $\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})$ we have the following interaction force function

$$F_{\text{ij}}^{\mathrm{T}}+F_{\text{ij}}^{\text{eI}}=-\nabla \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})=-\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}}.$$

(7)

The interaction force increment between t and t+Δ t is

$$\begin{array}{ll}\Delta F_{\text{ij}}^{\mathrm{T}}+\Delta F_{\text{ij}}^{\text{eI}}\hfill & \frac{\text{d}(F_{\text{ij}}^{\mathrm{T}}+F_{\text{ij}}^{\text{eI}})}{\text{d}\mathrm{t}}\Delta \mathrm{t}\hfill \\ & =-\frac{\text{d}}{\text{d}\mathrm{t}}(\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}})\Delta \mathrm{t}\hfill \\ & =-\frac{{\partial }^2\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{{(\partial {\text{r}}_{\mathrm{i}\mathrm{j}})}^2}\frac{\text{d}{\text{r}}_{\mathrm{i}\mathrm{j}}}{\text{d}\mathrm{t}}\Delta \mathrm{t}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}}-\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{\text{d}}{\text{d}\mathrm{t}}(\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}})\Delta \mathrm{t}\hfill \\ & =-\frac{{\partial }^2\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{{(\partial {\text{r}}_{\mathrm{i}\mathrm{j}})}^2}\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}}-\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{1}{{\text{r}}_{\mathrm{i}\mathrm{j}}^2}(\frac{\text{d}{\text{r}}_{\mathrm{i}\mathrm{j}}}{\text{d}\mathrm{t}}{\text{r}}_{\mathrm{i}\mathrm{j}}-\frac{\text{d}{\mathrm{r}}_{\mathrm{i}\mathrm{j}}}{\text{d}\mathrm{t}}{\text{r}}_{\mathrm{i}\mathrm{j}})\Delta \mathrm{t}\hfill \\ & =-\frac{{\partial }^2\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{{(\partial {\text{r}}_{\mathrm{i}\mathrm{j}})}^2}\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}}-\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{1}{{\text{r}}_{\mathrm{i}\mathrm{j}}}(\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}-\Delta {\mathrm{r}}_{\mathrm{i}\mathrm{j}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\mathrm{r}}_{\mathrm{i}\mathrm{j}}})\hfill \\ & =-\frac{{\partial }^2\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{{(\partial {\text{r}}_{\mathrm{i}\mathrm{j}})}^2}(\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}+\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{T}})\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}}-\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{1}{{\text{r}}_{\mathrm{i}\mathrm{j}}}[(\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}+\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{T}})-(\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}+\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{T}})\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\mathrm{r}}_{\mathrm{i}\mathrm{j}}}],\hfill \end{array}$$

(8)

where ${\Delta \text{r}}_{\text{ij}}^{\mathrm{e}}$ and $\Delta {\text{r}}_{\text{ij}}^{\mathrm{e}}$ are the change of the scalar and vector distance ${\mathrm{r}}_{\text{ij}}$ and ${\mathrm{r}}_{\text{ij}}$ caused by the external conditions, respectively; ${\Delta \text{r}}_{\text{ij}}^{\mathrm{T}}$ and $\Delta {\text{r}}_{\text{ij}}^{\mathrm{T}}$ are the change of the scalar and vector distance ${\mathrm{r}}_{\text{ij}}$ and ${\mathrm{r}}_{\text{ij}}$ without external loads, respectively. From (8) we know that the increment of the interaction force between the i^th particle and the j^th particle caused by the external conditions is

$$\Delta F_{\text{ij}}^{\text{eI}}=-\frac{{\partial }^2\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{{(\partial {\text{r}}_{\mathrm{i}\mathrm{j}})}^2}\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\text{r}}_{\mathrm{i}\mathrm{j}}}-\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{1}{{\text{r}}_{\mathrm{i}\mathrm{j}}}(\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}-\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}}}{{\mathrm{r}}_{\mathrm{i}\mathrm{j}}}).$$

(9)

Thus, at the k^th time step, the total interaction force of the i^th particle, caused by the external conditions, can be expressed by

$$F_{\text{ik}}^{\text{eI}}={\displaystyle \sum _{\mathrm{j}}{F}_{\text{ijk}}^{\text{eI}}}=F_{\text{ik}-1}^{\text{eI}}+{\displaystyle \sum _{\mathrm{j}}\Delta F_{\text{ijk}-1}^{\text{eI}}}=F_{\text{ik}-1}^{\text{eI}}-{\displaystyle \sum _{\mathrm{j}}[}\frac{{\partial }^2\mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{{(\partial {\text{r}}_{\mathrm{i}\mathrm{j}})}^2}\Delta {\text{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}^{\mathrm{e}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}}{{\text{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}}+\frac{\partial \mathrm{u}({\text{r}}_{\mathrm{i}\mathrm{j}})}{\partial {\text{r}}_{\mathrm{i}\mathrm{j}}}\frac{1}{{\text{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}}(\Delta {\text{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}^{\mathrm{e}}-\Delta {\text{r}}_{\mathrm{i}\mathrm{j}}^{\mathrm{e}}\frac{{\text{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}}{{\mathrm{r}}_{\mathrm{i}\mathrm{j}\mathrm{k}-1}})],$$

(10)

${\Delta \text{r}}_{\text{ijk}-1}^{\mathrm{e}}$ and ${\Delta \text{r}}_{\text{ijk}-1}^{\mathrm{e}}$ are the change of the scalar distance ${\mathrm{r}}_{\text{ij}}$ and the vector distance ${\text{r}}_{\text{ij}}$ at the k-1^th time step caused by the external conditions. ${\Delta \text{r}}_{\text{ijk}-1}^{\mathrm{e}}$ and ${\Delta \mathbf{\text{r}}}_{\text{ijk}-1}^{\mathrm{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:

$${{(\text{r}}_{\text{ij}}+\Delta {\mathrm{r}}_{\text{ij}})}^2={\mathrm{r}}_{\text{ij}}^2+{(\Delta {\mathrm{v}}_{\text{ij}}\Delta \text{t})}^2+2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ij}}.{\text{r}}_{\text{ij}},$$

(11)

and

$${{(\text{r}}_{\text{ij}}+\Delta {\mathrm{r}}_{\text{ij}}^{\mathrm{T}})}^2={\mathrm{r}}_{\text{ij}}^2+{(\Delta {\mathrm{v}}_{\text{ij}}^{\mathrm{T}}\Delta \text{t})}^2+2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ij}}^{\mathrm{T}}.{\text{r}}_{\text{ij}}.$$

(12)

Subtracting (12) from (11) and considering ${\Delta \text{r}}_{\text{ij}}^{\mathrm{e}}=\Delta {\mathrm{r}}_{\text{ij}}-\Delta {\mathrm{r}}_{\text{ij}}^{\mathrm{T}}$, after some mathematical manipulations, we obtain:

$$\begin{array}{ll}{\Delta \text{r}}_{\text{ij}}^{\mathrm{e}}{(2\text{r}}_{\text{ij}}+{\Delta \text{r}}_{\text{ij}}+\Delta {\mathrm{r}}_{\text{ij}}^{\mathrm{T}})\hfill & ={{(\text{r}}_{\text{ij}}+\Delta {\mathrm{r}}_{\text{ij}})}^2-{{(\text{r}}_{\text{ij}}+\Delta {\mathrm{r}}_{\text{ij}}^{\mathrm{T}})}^2\hfill \\ & ={(\Delta {\mathrm{v}}_{\text{ij}}\Delta \text{t})}^2-{(\Delta {\mathrm{v}}_{\text{ij}}^{\mathrm{T}}\Delta \text{t})}^2+2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ij}}.{\text{r}}_{\text{ij}}-2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ij}}^{\mathrm{T}}.{\text{r}}_{\text{ij}}\hfill \\ & =(\Delta {\mathrm{v}}_{\text{ij}}^2-{({\text{Δv}}_{\text{ij}}^{\mathrm{T}})}^2)\Delta {\mathrm{t}}^2+2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}.{\text{r}}_{\text{ij}}\hfill \\ & =\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}.(2\Delta {\text{v}}_{\text{ij}}-\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}})\Delta {\mathrm{t}}^2+2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ij}}^{\mathrm{e}}.{\text{r}}_{\text{ij}}.\hfill \end{array}$$

(13)

On the left hand side of (13) both ${\Delta \text{r}}_{\text{ij}}^{\mathrm{e}}$ and ${\mathrm{r}}_{\text{ij}}$ are scalar lengths. Usually $\Delta \text{t}$ is small and both $\Delta {\mathrm{r}}_{\text{ij}}$ and $\Delta {\mathrm{r}}_{\text{ij}}^{\mathrm{T}}$ are much smaller than ${\mathrm{r}}_{\text{ij}}$. We should notice here the smaller $\Delta {\mathrm{r}}_{\text{ij}}$ and $\Delta {\mathrm{r}}_{\text{ij}}^{\mathrm{T}}$ are ensured by the size of $\Delta \text{t}$, so we do not need to assume the flow velocity is small. Considering that all the variables in (10) correspond to the k-1^th time step, we add a subscript k-1 to represent the value at this time step and then equation (13) can be simplified to

$$\Delta {\mathrm{r}}_{\text{ijk}-1}^{\mathrm{e}}=(\Delta {\text{v}}_{\text{ijk}-1}^{\mathrm{e}}.(2\Delta {\text{v}}_{\text{ijk}-1}-\Delta {\text{v}}_{\text{ijk}-1}^{\mathrm{e}})\Delta {\mathrm{t}}^2+2\Delta \mathrm{t}\Delta {\text{v}}_{\text{ijk}-1}^{\mathrm{e}}.{\text{r}}_{\text{ijk}-1})/(2{\mathrm{r}}_{\text{ijk}-1}).$$

(14)

The procedure of the algorithm is: at the k^th time step, for any particle i, we calculate $\Delta {\mathrm{r}}_{\text{ijk}-1}^{\mathrm{e}}$ by formulae (14) and $\Delta {\text{r}}_{\text{ijk}-1}^{\mathrm{e}}$ = $\Delta {\text{v}}_{\text{ijk}-1}^{\mathrm{e}}\Delta \mathrm{t}$, then the force exerted on this particle at this moment $F_{\text{ijk}}^{\mathrm{e}}$ and $F_{\text{ik}}^{\mathrm{e}}$ by (10). Finally, the flow velocity increment $\Delta {\text{v}}_{\text{ik}}^{\mathrm{e}}$ and the flow velocity ${\text{v}}_{\text{ik}}^{\mathrm{e}}$ can be calculated by (5) and (6), respectively. When using (14) to calculate $\Delta {\mathrm{r}}_{\text{ijk}-1}^{\mathrm{e}}$ the term $\Delta {\text{v}}_{\text{ijk}-1}^{\mathrm{e}}$ is the velocity increment corresponding to the k-1^th 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

$${\mathrm{m}}_{\mathrm{i}}{\dot{\text{v}}}_{\mathrm{i}}=F_{\mathrm{i}}+{\mathrm{m}}_{\mathrm{i}}{g}_{\mathrm{i}}-\mathsf{\alpha }({\text{v}}_{\mathrm{i}}-{\text{v}}_{\mathrm{i}}^{\mathrm{e}}).$$

(15)

This is so called Gaussian thermostat method for nonequilibrium flow problems [25,26]. For equilibrium flow problems, both $g_{\mathrm{i}}$ and ${\text{v}}_{\mathrm{i}}^{\mathrm{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 ${\dot{\text{v}}}_{\mathrm{i}}$ is the derivative vector of the velocity of the i^th particle with respect to the time variable. ${\text{m}}_{\mathrm{i}}{g}_{\mathrm{i}}$ is the external force vector$F_{\mathrm{i}}^{\text{eII}}$. $F_{\mathrm{i}}$ is the total interaction force exerted on the i^th particle. It is composed of two parts, $F_{\mathrm{i}}=F_{\mathrm{i}}^{\mathrm{T}}+F_{\mathrm{i}}^{\text{eI}}$. $F_{\mathrm{i}}^{\mathrm{T}}$ is the interaction force without external loads. $F_{\mathrm{i}}^{\text{eI}}$ is the additional interaction force caused by the external force. We do not need to calculate $F_{\mathrm{i}}$ separately. This means that during the simulation, we do not need to use the sum of $F_{\mathrm{i}}^{\mathrm{T}}$ and $F_{\mathrm{i}}^{\text{eI}}$ to get $F_{\mathrm{i}}$. We can calculate $F_{\mathrm{i}}$ according to $F_{\mathrm{i}}={\displaystyle \sum _{\mathrm{j}\ne \mathrm{i}}(}{F}_{\text{ij}}^{\mathrm{T}}+F_{\text{ij}}^{\text{eI}})=-{\displaystyle \sum _{\mathrm{j}\ne \mathrm{i}}\nabla }u(r_{ij})$ and the position of particles directly. $\mathsf{\alpha }$ 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 ${\mathrm{m}}_{\mathrm{i}}{g}_{\mathrm{i}}$, that is $F_{\text{ik}}^{\text{eII}}={\text{m}}_{\mathrm{i}}{g}_{\mathrm{i}}$. Here g_i 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 m_f = m_w = 6.69x10^-26 kg; the temperature T = 84K; the interaction strength parameter between argon molecules is ε_f = 1.657x10⁶ 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 v^e 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 10¹³ 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 10⁹ 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 10¹² 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 10⁸ 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 v^e.

Convergence and Efficiency analysis

In this new algorithm, we only used the linear terms of the development of the force function with respect to ${\Delta \text{r}}_{\text{ij}}^{\mathrm{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 $\Delta {\mathrm{r}}_{\text{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 $\left|({\text{v}}_{\mathrm{maxLA}}-{\text{v}}_{\mathrm{maxTA}})/{\text{v}}_{\mathrm{maxTA}}\right|$ 100% and the horizontal coordinate is size of the time step $\Delta \mathrm{t}$(s), where v_{max LA} is the maximal velocity of this paper and v_{max 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.2x10⁵ 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 x10⁵ 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 v^e 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 v^e 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.