Molecular Dynamics Simulations of Pressure-Driven Flows and Comparison with Acceleration-Driven Flows

Abstract

We use molecular dynamics to simulate fluid flows between two parallel plates with constant wall temperature. Unlike the usual approach in molecular dynamics, instead of applying an external force on the molecules, the periodic boundary conditions are modified to create a pressure difference between the inlet and the outlet sections of the computational domain. The simulation results include velocity, pressure, density, and temperature profiles obtained by the new method. These results are compared with approximate solutions for nonisothermal Poiseuille flows. The method is also applied to simulate a flow in a rib-roughened channel.

1. Introduction

Over the past decades, nano/microfluidic systems have seen a rapid development and contribute to a revolution in many areas such as biotechnology and medicines. At the same time, calculation methods are also needed to analyze the reliability of the systems. The main challenge of the nano/microfluidics computations is that the studied system has a considerable surface/volume ratio so that the surface phenomena become highly important. These size effects cannot be solved by the traditional computational method for fluid dynamics.

Depending on the nature of the problem, there are different suitable computational methods, such as those based on the Navier-Stokes and energy equations, the Burnett or super-Burnett models [1–3], and the direct simulation Monte Carlo methods [4]. For nanoflows, the molecular dynamics (MD) method is one of the most accurate methods since realistic interactions between particles or between particles and boundaries may be accounted for.

One of the fundamental and practical problems in fluid mechanics is the steady-state flow between two parallel plates, the plane Poiseuille flow. Under the incompressibility assumption, the problem has a straightforward analytical solution. It is also studied extensively by other methods (Karniadakis et al. [5]). Most of the molecular dynamic simulations reported in the literature ([6–9], e.g.) concern flows driven by acceleration body force, rather than pressure driven flow. These limitations are due to the usage of the periodic boundary conditions which cannot generate pressure gradient inside the simulation domain. Several authors (Lupkowski and Van Swol [10], Li et al. [11], and Sun and Ebner [12]) have developed new techniques to mimic flows induced by pressure difference. Lupkowski and Van Swol [10] placed two rigid walls at the inlet and outlet and applied external forces on them. Li et al. [11] used a fictitious membrane which allows atoms to pass from one direction and forces atoms from the other direction to be elastically reflected with a given probability. Sun and Ebner [12] simulated 2D fluid flows in a long box while controlling the temperature at the two ends.

The present paper aims at presenting a new molecular dynamics approach based on pressure differences. Section 2 gives the problem statement and presents some approximate analytical solutions reported in the archival literature for Poiseuille flows. The pressure expression and the pressure difference from the atomistic point of view are also discussed. Section 3 is focused on the implementation of periodic boundary conditions linked to the pressure difference concept and its implementation in an MD code. The results of simulations are then compared with approximate analytical solutions of the Navier-Stokes and energy equations in Section 4. The application of the present method to fluids other than ideal gases and to rib-roughened channels is also discussed. Finally, some conclusions and perspectives are given in Section 5.

2. Problem Description and Theoretical Background

2.1. Problem Description

The objective of the present work is to study the fluid flows through a nanochannel. In the cartesian coordinate system $O x y z$ , we assume that the flow direction is along the x-axis. The length, width, and height of the channel are denoted by L, B, and H in directions x, y, and z, respectively. The channel width B being very large with respect to the other dimensions, it will be treated as infinite. To produce a fluid flow from the inlet to the outlet, we keep the inlet pressure p_i at higher level than the outlet p_o. The fluid can exchange energy with the channel walls which are maintained at constant temperature, T_w.

The Knudsen number, Kn, is a quantity to identify the fluid flow regime: from continuum (small Kn) to molecular flow regime (high Kn). The range of Kn chosen in the present study, $0.02 \leq Kn \leq 0.1$ , belongs to the slip flow regime according to [5], that is, the fluid continuum mechanics is still useful provided that slip boundary conditions are applied at the walls. The validity of the Navier-Stokes equations, as shown by many authors (see [7], e.g.), can therefore be extended to Poiseuille and Couette flows for the present range of Kn. The Knudsen number is defined as

Kn = \frac{λ}{H},

(1)

where λ is the mean free path, and H is a characteristic scale of the channel. In this work, λ is estimated by

λ = \frac{1}{\sqrt{2} n π σ^{2}} .

(2)

The quantities n and σ in (2) are the density number and the diameter of the fluid molecules, respectively. For Lennard Jones (LJ) fluids, the molecular diameter value σ is approximately taken equal to the σ parameter appearing in the LJ potential formula introduced in Section 4.1.

In what follows, the terms “pressure-driven flow” and “acceleration-driven flow” are employed to define two particular cases: (i) the body force ργ applied on a unit fluid volume is zero, and the pressure p decreases along the flow direction, $ρ γ = 0$ , $d p / d x \neq 0$ ; (ii) the inverse situation where $ρ γ \neq 0$ , $d p / d x = 0$ . We shall also consider some available analytical solutions of the Navier-Stokes equations reported in the literature. For both flows, the fluid velocity v in a cross-section which accounts for slip effects (Navier boundary conditions) at the wall admits the parabolic form

\begin{matrix} v (\hat{z}) = v (0) [1 - a_{2} {\hat{z}}^{2}], \hat{z} = \frac{z}{H}, \end{matrix}

(3)

where a₂ is a positive dimensionless constant related to pressure gradient (or body force), viscosity, slip length, and so forth. It is noted that (3) is generally obtained for incompressible fluids, but it is also valid for compressible fluid flows in a long channel [13]. In the latter case, v and a₂ are functions of the streamwise coordinate x.

An approximation solution of the energy equation can be obtained by neglecting the convective term [14]. The temperature T becomes then a quartic function of the coordinate z,

T (z) = T_{w} [b_{0} - b_{4} {\hat{z}}^{4}] .

(4)

Equation (4) is based on the incompressibility assumption. For compressible fluid flows in a long channel, Cai et al. [15] used a perturbation technique to derive the following temperature distribution:

T (x, z) = T_{w} [c_{0} - c_{2} {\hat{z}}^{2} - c_{4} {\hat{z}}^{4}] .

(5)

The dimensionless coefficients $c_{i}, b_{i}$ in (4) and (5) are related to pressure gradient (or body force), viscosity, conductivity, specific heat, and so forth. In (5), the coefficients c_i depend on the position x. We also note that for acceleration-driven flow, Todd and Evans [6] suggested a correction to the temperature profile given by (4) as a sextic function of z,

T (z) = T_{w} [d_{0} - d_{2} {\hat{z}}^{2} - d_{4} {\hat{z}}^{4} - d_{6} {\hat{z}}^{6}] .

(6)

In Section 4, MD simulations are used to reexamine the validity of the velocity and temperature profiles given by (3)–(6) and to determine the coefficients a_i, b_i, c_i, and d_i by curve fitting.

2.2. Pressure Tensor and Pressure Difference

Before proceeding the molecular dynamics simulation, let us look at the definition of pressure tensor from the atomistic point of view. Using statistical mechanic theory, Irving and Kirkwood [16] derived the following pressure tensor decomposition (IK):

P = P_{K} + P_{V} .

(7)

The kinetic term $P_{K}$ depends on the molecular square velocity, while the potential term $P_{V}$ depends on pairwise interactions between molecules. In a system of N molecules, $P_{K}$ and $P_{V}$ read

\begin{gathered} P_{K} (r) = \sum_{i = 1}^{N} 〈 m (v_{i} - \overset{̅}{v}) \otimes (v_{i} - \overset{̅}{v}) δ (r_{i} - r) 〉, \\ P_{V} (r) = - \frac{1}{2} \sum_{i}^{N} \sum_{j \neq 1}^{N} 〈 r^{(i j)} \otimes f^{(i j)} O_{i j} δ (r_{i} - r) 〉 . \end{gathered}

(8)

The term inside the angular bracket $〈 \cdot 〉$ denotes the ensemble average, and ⊗ denotes the dyadic product. The terms $r_{i}$ , $v_{i}$ , and $\bar{v}$ are respectively the position vector, the velocity of the particle i, and the mean velocity. The distance vector and interaction force vector between two molecules $(i)$ and $(j)$ are denoted, respectively, as $r^{(i j)}$ and $f^{(i j)}$ . For LJ fluids, the two vectors $r^{(i j)}$ and $f^{(i j)}$ are colinear, and the force magnitude $f^{(i j)}$ is derived from the interaction potential in function of the distance between the two molecules

f^{(i j)} = \frac{\partial V (r^{(i j)})}{\partial r^{(i j)}}, r^{(i j)} = ‖ r_{i} - r_{j} ‖ .

(9)

In (8), the notation $δ (r_{i} - r)$ is used for the delta Dirac function at $r_{i}$ , and the expression of the operator O_i,j reads

\begin{matrix} O_{i j} = 1 - \frac{1}{2!} r^{(i j)} \frac{\partial}{\partial r} + \dots + \frac{1}{n!} {[- r^{(i j)} \frac{\partial}{\partial r}]}^{n - 1} + \dots . \end{matrix}

(10)

Equations (10) and (8) show that the derivation of IK pressure tensor involves an infinite sum of high-order derivatives of the delta function and ensemble average, not suitable for MD computations. A more convenient form of the pressure tensor and an associated calculation method, the method of plane (MOP), was proposed by Todd et al. [17] and by Evans and Morriss [18]. When the fluid density is uniform, the operator O_i,j is reduced to unity, and (8) becomes (see [17, 19–22])

\begin{gathered} P_{K} (r) = \frac{1}{V} \sum_{i} m (v_{i} - \overset{̅}{v}) \otimes (v_{i} - \overset{̅}{v}), \\ P_{V} (r) = - \frac{1}{V} \sum_{i < j} r^{(i j)} \otimes f^{(i j)}, \end{gathered}

(11)

where V is the volume of the fluid element located at $r$ . In the literature, (11) is sometimes referred to as IK1 pressure [17, 18]. In (8) and (11), the pressure can be decomposed into two scalar quantities p_K and p_V as follows:

p = p_{K} + p_{V}, p_{K} = \frac{1}{3} tr P_{K}, p_{V} = \frac{1}{3} tr P_{V} .

(12)

It is clear that for low-density fluids, high temperature, and small interaction energy ∊, the pressure is reduced to the kinetic part p_K, as for ideal gases. For the other fluids, the potential part p_V cannot be neglected. The pressure difference Δp between two points at distance Δl reads

Δ p = Δ p_{K} + Δ p_{V} = k_{B} Δ (n T) + Δ p_{V} .

(13)

The pressure components P_x,x, P_y,y, P_z,z, number density n, and stream velocity $\bar{v}$ are invariant in direction y. From a microscopic point of view, taking $Δ l = Δ y$ , we must have

\begin{gathered} Δ 〈 v_{x}^{2} 〉 = Δ 〈 v_{y}^{2} 〉 = Δ 〈 v_{z}^{2} 〉 = 0, \\ Δ 〈 v_{x} 〉 = Δ 〈 v_{y} 〉 = Δ 〈 v_{z} 〉 = 0 . \end{gathered}

(14)

As it will be discussed in the following, all these conditions are satisfied in our MD simulation algorithm if we apply periodicity for the boundary conditions in the y direction. However, n, p, and T can all vary in the flow direction x. We note that in the steady-state regime, all quantities like T, p, p_V, and n are stable with time, and thus their differences $Δ T$ , Δp, $Δ p_{V}$ , and Δn are also stable. In particular, (13) shows that the finite variation of temperature $Δ T$ along x is related to those of other quantities via the expression

Δ T = \frac{Δ p - Δ p_{V} - k_{B} T Δ n}{k_{B} n + k_{B} Δ n} .

(15)

If the stream velocity is much smaller than the thermal velocity, we can obtain a simplified expression for the temperature difference as

Δ T = \frac{m}{3 k_{B}} Δ 〈 v^{2} - {\overset{̅}{v}}^{2} 〉 ≃ \frac{m}{3 k_{B}} Δ 〈 v^{2} 〉 .

(16)

The objective of the MD method discussed in what follows is to maintain the differences in squared velocity and, thus in temperature, in direction x. By that way, we can indirectly generate the pressure difference Δp. In the case where density ρ is invariant in x-direction (incompressible flow assumption), the temperature difference is proportional, for an ideal gas, to the pressure difference since

Δ p = \frac{ρ}{3} Δ 〈 v^{2} 〉

(17)

3. Numerical Procedure

In molecular dynamics, periodic boundary conditions (PBC) applied to velocities are powerful techniques to reduce the study of a large system to a smaller one far from the edge. Considering a simulation domain as a cube, PBC requires that if a molecule passes through one face, it reappears on the opposite face with the same velocity. Obviously, there is no difference in pressure, density, or temperature between any two opposite faces of the simulation domain. Consequently $Δ p = 0$ , as it will be shown in the next section.

3.1. Modified Periodical Boundary Conditions

In our problem, all the virial pressure components P_x,x, P_y,y, and P_z,z do not vary along y, which can be satisfied by the traditional PBC applied on the faces $y = \pm B / 2$ . However, to create a pressure difference along the x-axis, we must develop a strategy that produces a constant difference between the squared velocity of the molecules crossing the face × = 0 and those crossing the face $x = L$ .

Usually, to keep unchanged the total number of molecules in the domain whenever a molecule goes out of the domain, we must insert another molecule inside. Thus, if we relax the PBC constraints, there are many ways to assign velocities to ingoing molecules, either independent of or dependent on the outgoing ones. In this paper, we generalize the PBC to account for the pressure difference by maintaining the difference in squared velocity at the two opposite faces × = 0 and $x = L$ .

From (15) and (16), we know that a constant pressure difference is related to a constant difference in squared velocity. A constant parameter δ is thus used to create a constant difference in squared velocity between the inlet and outlet faces (see Figure 1). In order to apply a velocity square difference equal to a constant $3 δ L$ between the faces × = 0 and $x = L$ , we modify the periodic boundary conditions so that if a molecule goes through the face $x = L$ with a velocity $v$ , we insert a molecule at × = 0 with a velocity $v^{'}$ satisfying the following conditions:

v^{'} = α v, α = \sqrt{1 + \frac{3 δ L}{v^{2}}} .

(18)

Analogously, if a molecule crosses the face × = 0 with a velocity $v$ , it reappears at $x = L$ with a velocity $v^{'}$ so that

v^{'} = α v, α = \sqrt{1 - \frac{3 δ L}{v^{2}}} .

(19)

Since $3 δ L$ is positive, it may happen that $v^{2} - 3 δ L < 0$ . It is then impossible to find v′ satisfying (19). The boundary conditions applied to the simulation domain are summarized in Table 1. However, this effect can be assumed negligible for flows with relatively small speed (of the order $1 %$ of the thermal speed), as seen in most MEMS/NEMS devices. In the present simulations, the number of molecules that do not satisfy (19) were found less than $1 %$ of the outgoing molecules. In order to simulate low-speed flows, small values of δ were used.

Table 1:

Summary of the boundary conditions applied on the simulation domain.

	Boundary conditions

Faces x	Equations (18)-(19)
Faces y	PBC
Faces z	Thermal walls

Figure 1:

Modified periodic boundary conditions. The velocity v between the outgoing molecules is different from the ingoing molecules ((18) and (19)).

4. Results

4.1. Molecular Dynamic Simulations of Ideal Gas Flows

To model the interaction between the molecules, the 6–12 Lennard Jones potential is used

V (r) = 4 ϵ [{(\frac{σ}{r})}^{12} - {(\frac{σ}{r})}^{6}],

(20)

where ∊ is the depth of the potential well, σ is the finite distance at which the interparticle potential is zero, and r is the distance between the particles.

The potential well depth parameter is $∊ = 0.1 ∊_{0}$ with ∊₀ being the reference potential well depth. Two cutoff distances were used in this study: $r_{c} = 2^{1 / 6} σ$ [17] corresponds to repulsive interactions between molecules, while $r_{c} = 2.5 σ$ [21, 23–25] corresponds to attractive-repulsive interactions. In general, the present method is valid for any pair potential. However, the simulation conditions of Table 2, that is, $(∊ = 0.1 ∊_{0}, r_{c} = 2^{1 / 6} σ)$ , lead to a very small potential part of the pressure p_V (less than 1% of the kinetic part p_K) making the behavior of the fluid close to an ideal gas.

Table 2:

Input data for the first series of MD simulations.

N	L [σ]	B [σ]	H [σ]	Kn	δ $[∊_{0} / m σ]$	γ $[∊_{0} / m σ]$

2400	86.2	6.5	43.1	0.05	10 $^{- 3}$	4.3 × $1 0^{- 4}$
600	43.1	6.5	21.5	0.1	2 × $1 0^{- 3}$	8.6 × $1 0^{- 4}$

$(∊, r_{c}) = (0.1 ∊_{0}, 2^{1 / 6} σ)$ , $\bar{n} = 0.1 σ^{- 3}$ , and $T = 2 ∊_{0} / k_{B}$ .

The two walls are modeled as thermally diffusive walls at the same and constant temperature. Whenever a molecule arrives at a wall, it is reflected with the velocity corresponding to the wall temperature $T_{w} = 2 ∊_{0} / k_{B}$ and with a random direction. In this work, we do not use thermostats, and the fluid can exchange energy with the wall due to the wall model used here.

In our simulations, the global number density $\bar{n}$ is kept fixed, $\bar{n} = 0.1 σ^{- 3}$ in every configurations. Other geometric parameters like L, B, H, and the number of the molecules N are changed as shown in Table 2. The values of δ are chosen such that the difference of the square of the molecular velocity between the inlet and outlet is small enough so that the stream velocity is much smaller than the thermal speed. From the microscopic viewpoint, the local stream velocity is the local average velocity of the molecules, and from the macroscopic viewpoint, it is the flow speed $v$ in the Navier-Stokes equations and thus connected to the pressure gradient $\nabla p$ (or volume force ργ). On the other hand, the thermal speed is the root mean square of the molecular velocity and is directly connected to the local temperature T. The slope of the simulated total pressure distribution along the channel (see Figures 2 and 3) is used to determine the pressure gradient $- Δ p / Δ x$ and the equivalent external acceleration $\bar{ρ} γ$ used in the acceleration-driven flow simulations. The quantity $\bar{ρ}$ is the average mass density of the system, $\bar{ρ} = m \bar{n}$ . The computation with Leap-Frog Verlet integration scheme is carried out for 10⁸ time steps from equilibrium, each of which is equal to 0.005 unit time $τ = \sqrt{m σ^{2} / ∊_{0}}$ . The height and length of the channel are divided into 50 layers to determine accurately the distribution of local velocity $v$ , temperature T, and pressure tensor $P$ . The local velocity $v$ measured in one bin located at $r$ is given by

v (r) = \frac{1}{n_{bin}} \sum_{r_{i} \in bin} v_{i},

(21)

n_bin being the particle number in the bin. The local pressure $P$ (IK1-model) is

\begin{matrix} P (r) & = \frac{1}{V_{bin}} \\ \times [\sum_{r_{i} \in bin} m (v_{i} - \overset{̅}{v}) \otimes (v_{i} - \overset{̅}{v}) - \frac{1}{2} \sum_{r_{i} \in bin} r^{(i j)} \otimes f^{(i j)}] . \end{matrix}

(22)

The pressure components were also computed by using the method of plane [17, 26]. With α, β being x, y, or z, the pressure component P_β,α along direction β and acting on the area element A_α normal to the axis α is defined by

\begin{matrix} P_{β α} (r) & = \frac{1}{τ A_{α}} \sum_{r_{i} \in A_{α}} \frac{m (v_{i β} - {\overset{̅}{v}}_{i β}) (v_{i α} - {\overset{̅}{v}}_{i α})}{v_{i α}} \\ + \frac{1}{τ A_{α}} \int d t \sum_{r_{i j} \cap A_{α}, i < j} f_{β}^{(i j)} . \end{matrix}

(23)

In (23), the first sum is for all molecules crossing the area element A_α over the period τ, and the second sum is for all pairs whose distance vectors cut the area element A_α. The temperature T of one bin is then calculated by

\begin{matrix} T (r) = \frac{m}{3 k_{B}} {〈 {‖ v - \bar{v} ‖}^{2} 〉}_{bin} . \end{matrix}

(24)

In what follows, the x-distributions of pressure, temperature, and density correspond to an average on y and z of the local quantities $p (r)$ and $n (r)$ . They are denoted as $\bar{p} (x) = {〈 p (r) 〉}_{y, z}$ and $\bar{n} (x) = {〈 n (r) 〉}_{y, z}$ , respectively. Analogously, the profile along z of velocity, temperature, and density corresponds to an average on x and y of local quantities with the notations: $\bar{v} (z) = {〈 v (r) 〉}_{x, y}$ , $\bar{T} (z) = {〈 T (r) 〉}_{x, y}$ , and $\bar{n} (z) = {〈 n (r) 〉}_{x, y}$ , respectively.

Figure 2:

Pressure distribution ( $\bar{p} (x)$ ) along the channel for pressure-driven flow case and acceleration-driven flow case (Kn = 0.05, $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ , $δ = 1 0^{- 3} ∊_{0} / m σ$ , and $γ = 4.3 \times 1 0^{- 4} ∊_{0} / m σ$ ). The full line represents the linear fit and is used to determine the pressure gradient $Δ p / Δ x$ . The pressure unit is $[∊_{0} / σ^{3}]$ .

Figure 3:

Distribution of pressure $\bar{p} (x)$ and pressure components ${\bar{P}}_{x x} (x)$ , ${\bar{P}}_{y y} (x)$ , and ${\bar{P}}_{z z} (x)$ along the channel for pressure-driven flow case (Kn = 0.05, $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ , and $δ = 1 0^{- 3} ∊_{0} / m σ$ ). The pressure unit is [ $∊_{0} / σ^{3}$ ].

From the simulations, the axial pressure distribution is plotted in Figure 2 for the pressure-driven and the acceleration-driven flows with $Kn = 0.05$ , $r_{c} = 2^{1 / 6} σ$ . As expected, for the acceleration-driven flow case, the pressure is constant. In contrast, the pressure distribution for pressure-driven flow decreases linearly. We note here that the slope of the pressure curve, equal to $Δ \bar{p} / Δ x$ , is used to compute the value γ in the acceleration flow simulation according to the formula $γ = - Δ \bar{p} / (\bar{ρ} Δ x)$ . Under the simulation conditions $\bar{n} = 0.1 σ^{- 3}$ , $T_{w} = 2 ∊_{0} / k_{B}$ , and $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ , the kinetic part p_K of the pressure dominates.

In order to know how the pressure difference is distributed in the x-, y-, z-directions, the pressure components ${\bar{P}}_{x x}$ , ${\bar{P}}_{y y}$ , and ${\bar{P}}_{z z}$ along x are plotted in Figure 3. All the three pressure components decrease with x, and the slopes are almost the same. The values of components ${\bar{P}}_{x x}$ , ${\bar{P}}_{y y}$ are very close, whereas the value of component ${\bar{P}}_{z z}$ is smaller because the fluid is confined in the z-direction. Although the difference between the pressure components is small (less than 4%), it implies that the velocity distribution deviates from the equilibrium distribution. Because the pressure $\bar{p}$ is the average value of ${\bar{P}}_{x x}$ , ${\bar{P}}_{y y}$ , ${\bar{P}}_{z z}$ , the $\bar{p}$ -curve lies between the others and has the same slope. It should be noted that (22) and (23) lead to relatively close results when considering the average axial pressure component ${\bar{P}}_{x x}$ (see Figure 4). Figure 5 shows that the fluid density distribution $\bar{n} (z)$ across the channel width varies slightly around the average density $\bar{n}$ , except near the walls where it can be as high as $1.4 \bar{n}$ . However, the phase diagram for LJ fluid [25] shows that the fluid near the wall is still in gas state. In our simulations, $\bar{n}$ and ∊ are chosen to be very small ( $\bar{n} = 0.1 σ^{- 3}$ , $∊ = 0.1 ∊_{0}$ ), and the fluid temperature near the wall is relatively large (the temperature of the wall is $T_{w} = 2 ∊_{0} / k_{B}$ ). The fluctuation of the density profile shows that the incompressibility assumption is no longer valid, especially near the channel walls. However, the fluctuations being localized, it is interesting to check if the analytical solutions presented in Section 2 still agree with the MD solutions.

Figure 4:

Distribution of pressure ${\bar{P}}_{x x}$ along the channel for pressure-driven flow case (Kn = 0.05, $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ , and $δ = 1 0^{- 3} ∊_{0} / m σ$ ) computed by (22) and (23). The pressure unit is [ $∊_{0} / σ^{3}$ ].

Figure 5:

Dimensionless density profile ( $\bar{n} (z)$ ) in half of the channel cross-section for pressure-driven simulations (Kn = 0.05, 0.1, $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ , $δ = 1 0^{- 3} ∊_{0} / m σ$ , and $2 \times 1 0^{- 3} ∊_{0} / m σ$ ).

Figure 6 shows that acceleration and pressure-driven flows exhibit parabola-like velocity profiles in the center region of the channel, in agreement with (3). Near the wall where the Knudsen layer dominates, the velocity distribution tends to deviate from the solution given by (3), and a velocity slip at the wall is predicted. Based on the kinetic theories, different models have been derived to explain the slippage [5], and generally, the slip effect becomes important when the Knudsen number increases. The little difference in velocity profile between the two types of flow shows that, despite their different microscopic natures, volume forces can be seen as equivalent to pressure gradients at the macroscopic scale.

Figure 6:

Dimensionless velocity profile in half of the channel cross-section for Kn = 0.05 and Kn = 0.1 ( $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ ). The full lines represent equation $v (\hat{z}) = v (0) (1 - a_{2} {\hat{z}}^{2})$ that best fit the numerical results ( $a_{2} = 0.18$ for Kn = 0.05 and $a_{2} = 0.29$ for Kn = 0.1).

About the temperature profile across the channel width, there is no visible difference between acceleration and pressure-driven formulations whatever Kn (Figure 7). The temperature profile cannot be fitted by any of the three approximate analytical solutions ((4)–(6)). The temperature is minimal at the center of the channel and increases rapidly near the wall, which does not agree with the approximate solutions obtained with the incompressible assumption ((4)-(5)). With the simulation conditions reported in Table 2, this reverse trend cannot thus be explained by using the incompressible flow equations. This anomaly was also observed by using DSMC method [27] and super-Burnett (SB) method [28]. Note that the flow speeds in our simulations are very small, as in [27, 28].

Figure 7:

Dimensionless temperature profile ( $\bar{T} (z)$ ) in half of the channel cross-section for two Kn values according to the flow model ( $r_{c} = 2^{1 / 6} σ$ , $∊ = 0.1 ∊_{0}$ ).

Next, we study the influence of γ and δ on the temperature profile. By increasing the acceleration parameter γ or the δ-parameter (squared velocity), the stream velocity is increased. A change in the form of the temperature profile is observed in Figures 8 and 9 as γ (or δ) increases. At high values of γ (or δ), the temperature profiles seem to be closer to the approximate solutions which predict a maximum at the center of the channel and a minimum at the walls. At small values of γ (or δ), the temperature profiles for the two flow types differ just a little. However, at high values of γ (or δ), the discrepancy becomes more important: the curvature at the center of channel for the pressure-driven flow case is higher than that for the acceleration-driven flow case. Using a curve-fitting procedure, we find that the temperature profile of the pressure-driven flow case agrees quite well with the quartic expression (5). The temperature for the acceleration-driven flow case does not fit well with (4). To explain this discrepancy, Todd and Evans [6] argued that the transport coefficients are not constant and that an additional cross-coupling between strain rate and heat flux exists. They proposed a correction to (5) with a $z^{6}$ -extraterm for the temperature profile (see (6)). The dashed line drawn in Figure 9 shows that the sextic polynomial fits very well the present numerical results.

Figure 8:

Variation of temperature profile when increasing the pressure gradient or parameter δ (Kn = 0.05, $r_{c} = 2^{1 / 6} σ$ , and $∊ = 0.1 ∊_{0}$ ). The dashed lines are based on equation $c_{0} - c_{2} {\hat{z}}^{2} - c_{4} {\hat{z}}^{4}$ that best fits the numerical results. $c_{0} = 3.07$ , $c_{2} = 4.16$ , and $c_{4} = 12.8$ for $γ = 0.1 ∊_{0} / m σ$ , $c_{0} = 1.74$ , $c_{2} = 2.5$ , and $c_{4} = 1.0$ for $γ = 0.05 ∊_{0} / m σ$ .

Figure 9:

Changes in the temperature profile for various values of the acceleration parameter γ (Kn = 0.05, $r_{c} = 2^{1 / 6} σ$ , and $∊ = 0.1 ∊_{0}$ ). The solid and the dashed lines correspond to equations $b_{0} - b_{4} {\hat{z}}^{4}$ and $d_{0} - d_{2} {\hat{z}}^{2} - d_{4} {\hat{z}}^{4} - d_{6} {\hat{z}}^{6}$ that best fit the numerical results. $b_{0} = 2.53$ , $b_{4} = 13.0$ , $d_{0} = 2.46$ , $d_{2} = 0.54$ , $d_{4} = 4.58$ , and $d_{6} = - 97.8$ for $γ = 0.1 ∊_{0} / m σ$ and $b_{0} = 1.48$ , $b_{4} = 4.57$ for $γ = 0.05 ∊_{0} / m σ$ .

4.2. Applications to More General Cases and Rib-Roughened Channels

The application of the method to cases for which the fluid particles interact strongly, that is, the potential pressure p_V is of the order of p_K, is considered in that section. The case of rib-roughened channels is also briefly discussed.

With the algorithm developed in the previous section, two cases were considered with the same parameters, except for the cutoff distance r_c and density number $\bar{n}$ . In the first case, $\bar{n} = 0.1 σ^{- 3}$ , and the cutoff distance is set to $r_{c} = 2.5 σ$ (attractive-repulsive interaction). In the second case, the fluid density is increased up to $\bar{n} = 0.8 σ^{- 3}$ while setting $r_{c} = 2^{1 / 6} σ$ as in the previous section. The well depths ∊ are the same in both cases and equal to the reference value, that is, $∊ = ∊_{0}$ . All the parameters for these two cases are summarized in Table 3.

Table 3:

Input data for MD simulations in Section 4.2.

$\bar{n}$ $[σ^{- 3}]$	N	L [σ]	B [σ]	H [σ]	δ $[∊_{0} / m σ]$	T _w $[∊_{0} / k_{B}]$

0.1	2400	86.2	6.5	43.1	$1 \times 1 0^{- 3}$	2
0.8	2400	86.2	3.2	43.1	$1 \times 1 0^{- 3}$	2

$(\bar{n}, ∊, r_{c}) = (0.1 σ^{- 3}, ∊_{0}, 2.5 σ)$ , $(0.8 σ^{- 3}, ∊_{0}, 2^{1 / 6} σ)$ .

In the first case ( $\bar{n} = 0.1 σ^{- 3}$ ), Figure 10 displays the kinetic part ${\bar{p}}_{K}$ and the virial pressure $\bar{p}$ according to (11) along the channel. Except a slight curvature in the ${\bar{p}}_{K}$ curve, both pressure curves decrease continuously and linearly. The relative difference between the slopes of both curves is less than 1%, which means that the gradient of the potential part $Δ {\bar{p}}_{V} / Δ x$ is also small in comparison with the kinetic part $Δ {\bar{p}}_{K} / Δ x$ even if ${\bar{p}}_{V}$ is not negligible ( $| {\bar{p}}_{V} / {\bar{p}}_{K} | ≃ 0.2$ ). In the second case ( $\bar{n} = 0.8 σ^{- 3}$ ), Figure 11 shows a totally reverse trend: both the virial pressure profile and its potential part decrease with x. The variation of the potential part contributes to $80 %$ of the virial pressure variation with $Δ \bar{p} / Δ x = - 0.0177$ and $Δ {\bar{p}}_{V} / Δ x = - 0.0159$ . We conclude that the method is convenient for generating virial pressure profiles which decrease in the axial direction.

Figure 10:

Axial variations of pressures $\bar{p} (x) = p_{K} + p_{V}$ and p_K for pressure-driven flow ( $\bar{n} = 0.1 σ^{- 3}$ , $r_{c} = 2.5 σ$ , and $∊ = ∊_{0}$ ). The pressure unit is [ $∊_{0} / σ^{3}$ ].

Figure 11:

Axial variations of pressures $\bar{p} (x) = p_{K} + p_{V}$ and p_V for pressure-driven flow ( $\bar{n} = 0.8 σ^{- 3}$ , $r_{c} = 2^{1 / 6} σ$ , and $∊ = ∊_{0}$ ). The pressure unit is [ $∊_{0} / σ^{3}$ ].

Next, we consider the velocity profiles across the channel section for these two cases. The velocity profiles (Figure 12) agree then with (3). The slip length for the case $\bar{n} = 0.8 σ^{- 3}$ is smaller than for the case $\bar{n} = 0.1 σ^{- 3}$ because the mean free path is decreased when the density number is increased. When computing the average pressure component ${\bar{P}}_{x x}$ along the channel, Figure 13 shows that there is no visible difference between (22) and (23).

Figure 12:

Dimensionless velocity profile in half of the channel cross-section for $\bar{n} = 0.1 σ^{- 3}$ and $\bar{n} = 0.8 σ^{- 3}$ . The dashed lines represent the equations $a_{0} - a_{2} z^{2}$ that best fit the numerical results.

Figure 13:

Axial variations of pressure component ${\bar{P}}_{x x}$ computed by (22) and (23). The pressure unit is [ $∊_{0} / σ^{3}$ ].

It should thus be emphasized that the present method is relevant to generate various fluid flows even if we control only the difference between the inlet- and outlet-squared molecular velocities. The potential part p_V, which depends on the interatomic interaction, can vary along the flow direction. As shown numerically, the variation of p_V contributes considerably to the pressure gradient.

The third case aimed at modeling the flow within a rib-roughened channel. The parameters are the same as in the second case ( $\bar{n} = 0.8 σ^{- 3}$ , $r_{c} = 2^{1 / 6} σ$ ), but the channel height is suddenly reduced in its middle part by inserting on both walls ribs of height h and length l. The rib-to-channel width ratio is equal to $h / H = 0.8$ and $l / L = 0.3$ (Figure 14). Figure 14 exhibits the onset of two vortices close to the rib corners at the upstream section. It is also shown that the fluid flow is highly nonuniform and characterized by wavy streamline patterns within the downstream region. From Figure 15, the axial pressure variations are rather different from those predicted for smooth channels. Although the variations of the axial pressure predicted by (22) and (23) are quite similar, there are considerable differences between them. The MOP pressure is $8 %$ smaller than the IK1 pressure. The changes in pressure at the downstream and upstream sections display also smoother pressure variations. This considerable difference is due to the strong inhomogeneity of the fluid within the channel. Generally, the MOP, based on the hydrodynamics equation, must be used in such situation.

Figure 14:

Streamlines in a rib-roughened channel. The color represents the fluid density field.

Figure 15:

Distribution of pressure component ${\bar{P}}_{x x}$ in a rib-roughened channel computed by (22) and (23). The pressure unit is [ $∊_{0} / σ^{3}$ ].

5. Conclusions

In this work, we have used a molecular dynamics method for simulating pressure-driven flows in channels. The main novelty lies in the modification of the periodic inlet and outlet velocity conditions: the difference in squared velocity between the ingoing and outgoing molecules of the simulation domain is maintained constant.

Plane Poiseuille flows were simulated, and the results were compared with approximate analytical solutions of the Navier-Stokes and energy equations reported in the literature. This new method has been proven to be realistic and could be applied to many problems with a nonconstant axial pressure gradient: flows around obstacles or in rough channels, for example. The method appears thus relevant for modeling compressibility effects as well as temperature variations along the flow direction, a domain unfulfilled by the usual molecular dynamics methods.

References

Struchtrup

, Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methods in Kinetic Theory, Springer, New York, NY, USA, 2005.

Liou

W. W.

and Fang

, Microfluid Mechanics, McGraw-Hill, New York, NY, USA, 2003.

Chapman

Cowling

T. G.

, and Burnett

, The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction, and Diffusion in Gases, Cambridge University Press, Cambridge, UK, 1990.

Bird

G. A.

, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon Press, 1994.

Karniadakis

Beskok

, and Aluru

, Microflows and Nanoflows: Fundamentals and Simulation, Springer, New York, NY, USA, 2005.

Todd

B. D.

and Evans

D. J.

, “Temperature profile for Poiseuille flow,” Physical Review E, vol. 55, no. 3, pp. 2800–2807, 1997.

Bhattacharya

D. K.

and Lie

G. C.

, “Molecular-dynamics simulations of nonequilibrium heat and momentum transport in very dilute gases,” Physical Review Letters, vol. 62, no. 8, pp. 897–900, 1989.

Barrat

J. L.

and Bocquet

, “Large slip effect at a nonwetting fluid-solid interface,” Physical Review Letters, vol. 82, no. 23, pp. 4671–4674, 1999.

Q. D.

Bercegeay

Lauriat

Léonard

, and Bonnet

, “A slip model for micro/nano gas flows induced by body forces,” Microfluidics and Nanofluidics, vol. 8, no. 3, pp. 417–422, 2010.

10.

Lupkowski

and Van Swol

, “Computer simulation of fluids interacting with fluctuating walls,” The Journal of Chemical Physics, vol. 93, no. 1, pp. 737–745, 1990.

11.

Liao

, and Yip

, “Coupling continuum to molecular-dynamics simulation: reflecting particle method and the field estimator,” Physical Review E, vol. 57, no. 6, pp. 7259–7267, 1998.

12.

Sun

and Ebner

, “Molecular-dynamics simulation of compressible fluid flow in two-dimensional channels,” Physical Review A, vol. 46, no. 8, pp. 4813–4818, 1992.

13.

Arkilic

E. B.

Schmidt

M. A.

, and Breuer

K. S.

, “Gaseous slip flow in long microchannels,” Journal of Microelectromechanical Systems, vol. 6, no. 2, pp. 167–178, 1997.

14.

Landau

L. D.

and Lifshitz

E. M.

, Fluid Mechanics, Pergamon Press, 1987.

15.

Cai

Sun

, and Boyd

I. D.

, “Gas flows in microchannels and microtubes,” Journal of Fluid Mechanics, vol. 589, pp. 305–314, 2007.

16.

Irving

J. H.

and Kirkwood

J. G.

, “The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics,” The Journal of Chemical Physics, vol. 18, no. 6, pp. 817–829, 1950.

17.

Todd

B. D.

Evans

D. J.

, and Daivis

P. J.

, “Pressure tensor for inhomogeneous fluids,” Physical Review E, vol. 52, no. 2, pp. 1627–1638, 1995.

18.

Evans

and Morriss

, Statistical Mechanics of Nonequilibrium Liquids, Cambridge University Press, Cambridge, UK, 2008.

19.

Hoover

W. G.

Hoover

C. G.

, and Lutsko

J. F.

, “Microscopic and macroscopic stress with gravitational and rotational forces,” Physical Review E, vol. 79, no. 3, Article ID 036709, 2009.

20.

Tsai

D. H.

, “The virial theorem and stress calculation in molecular dynamics,” The Journal of Chemical Physics, vol. 70, no. 3, pp. 1375–1382, 1979.

21.

Allen

M. P.

and Tildesley

D. J.

, Computer Simulation of Liquids, Oxford University Press, Oxford, UK, 1989.

22.

Hansen

J. P.

and McDonald

I. R.

, Theory of Simple Liquids, Academic Press, New York, NY, USA, 2006.

23.

Rapaport

D. C.

, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, UK, 2004.

24.

Bhattacharya

D. K.

and Lie

G. C.

, “Nonequilibrium gas flow in the transition regime: a molecular-dynamics study,” Physical Review A, vol. 43, no. 2, pp. 761–767, 1991.

25.

Frenkel

and Smit

, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, New York, NY, USA, 2002.

26.

Han

and Lee

J. S.

, “Method for calculating the heat and momentum fluxes of inhomogeneous fluids,” Physical Review E, vol. 70, no. 6, Article ID 061205, 2004.

27.

Zheng

Garcia

A. L.

, and Alder

B. J.

, “Comparison of kinetic theory and hydrodynamics for Poiseuille flow,” Journal of Statistical Physics, vol. 109, no. 3–4, pp. 495–505, 2002.

28.

, “Super-Burnett solutions for Poiseuille flow,” Physics of Fluids, vol. 15, no. 7, pp. 2077–2080, 2003.