Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 4, pp. 1141-1153, Warsaw 2017
DOI: 10.15632/jtam-pl.55.4.1141

ANALYSIS OF CONTROLLED MOLECULAR DYNAMIC FLOW IN

A CHANNEL WITH NON-EQUAL INLET AND OUTLET

CROSS-SECTIONAL AREAS

Hamed R. Najafi, S.M. Hossein Karimian

Department of Aerospace Engineering, Amirkabir University of Technology, Tehran, Iran

e-mail: hamedramezani@aut.ac.ir; hkarim@aut.ac.ir

Thermo-fluid properties are required for numerical modeling of nano/micro devices. These
properties are mostly obtained from the results of molecular dynamics (MD) simulations.
Therefore, efforts have been made to develop methods for numerical evaluation of fluid
properties such as pressure and velocity. One of the main challenges faced by numerical
simulations is to simulate steadymolecular flow in channels with non-equal inlet and outlet
boundaries. Currently, periodic boundary conditions at the inlet and outlet boundaries are
an inevitable condition in many steady flow molecular dynamics simulations. As a result,
a nano-channel with different cross sectional areas at the inlet and outlet could not be
simulated easily. Here, a method is presented to generate and control steadymolecular flow
in a nano-channel with different cross sectional areas at the inlet and outlet. The presented
method has been applied to a converging-diverging channel, and its performance has been
studied through qualitative and quantitative representation of flow properties.

Keywords: molecular dynamics, nano-channel, steady flow, non-equal inlet and outlet,
pressure

1. Introduction

Flow behavior at micro and nano scales has been the subject of interest in the recent years. As
a typical reduced-size fluid flow system, flow in a nanochannel can demonstrate different aspects
of nanofluid systems. Therefore, several investigations have been focused on the simulation of
nanochannel flows (Mi and Chwang, 2003).
Molecular Dynamics (MD) is a deterministic method to calculate the position of molecules

and their dynamic properties. Different potential functions have been introduced to represent
molecular forces in MD simulations (see for example Leach, 1999; Sadus, 2002). Extracting
macroscopic properties such as velocity, temperature and pressure from microscopic data has
also been a challenging issue (Allen and Tildesley, 1987; Karniadakis et al., 2002; Karimian et
al., 2011).
In earlier studies, the characterization of domain boundaries was not of the primary impor-

tance and thus they were treated as simple as possible. In most of the cases, therefore, periodic
boundary conditions were used in MD simulations (Stillinger and Rahman, 1974; Travis and
Evans, 1977; Koplik et al., 1989; Travis andGubbins, 2000; Nagayama andCheng, 2004). Howe-
ver, due to the rapid advancement of computational facilities in the recent years, new problems
with more complex boundaries have become a subject of MD simulations. As a result, imple-
mentation of realistic boundary conditions in MD simulation is now the center of attention,
especially in the field of fluid mechanics (Sun and Ebner, 1992; Huang et al., 2004; Hanasaki
and Nakatani, 2006; Huang et al., 2006).
One of the MD applications today that can practically be very beneficial is to present a

method for simulation of a specific known or predefined flow in a nanochannel with non-equal



1142 H.R. Najafi, S.M.H. Karimian

inlet andoutlet cross-sectional areas.Different approaches canbe found in the literature to create
flow fields in a nanochannel. The most simple and straightforward method is to implement an
external force or acceleration to move molecules in a specified direction (Koplik et al., 1988;
Fan et al., 2002; Ziarani and Mohamad, 2005; Sofos et al., 2009; Kamali and Kharazmi, 2011).
Other than this, some has used motion of a piston or plate upstream of the inlet to create flow
in nanotubes (Huang et al., 2004, 2006; Hanasaki and Nakatani, 2006). Implementation of flow
gradients such as the temperature gradient can also be used to create nanoflow (Han, 2008; Liu
and Li, 2010; Darbandi et al., 2011; Bao et al., 2015). In addition to these methods, in another
method,wallmotion has been used to createmolecular flow in nanochannels (Zhang et al., 2009;
Kim et al., 2010).

Implementation of a periodic boundary condition at the inlet and outlet is inevitable in all
the above-mentionedmethods except for the second one in which a piston upstream of the inlet
is used to move the molecules. However, the disadvantage of this approach is that the solution
is limited to the period at which the pistonmoves in its stroke. Also, since additional molecules
are to be taken into account in stroke space upstream of the inlet, extra computational efforts
are required in this method.

While all of the above-mentioned methods achieve the goal of creating molecular flow in
nanochannels, they fail to model steady flow in a channel with different inlet and outlet cross-
-sectional areas since they use a periodic boundary condition at the inlet and outlet. Also, none
of the above methods offers a proper approach to produce specific known or predefined flow
in nanochannels. Therefore, the main objective of present study is to investigate flow behavior
in a channel of non-equal inlet and outlet cross-sectional areas with controllable flow. In the
following sections, we first introduce the calculation method used here to evaluate macroscopic
properties in nanoflow. Then we will discuss details about the nanoflow controlled by a tempe-
rature gradient. Finally, we introduce amethod to solve molecular dynamic flow in a non-equal
inlet and outlet cross-sectional-area channel controlled by the temperature gradient.

Finally, a short review on boundary conditions is required since they are implemented on the
nanotube wall. The most common option is to model the walls by two rows of molecules. The
wall molecules are allowed to oscillate about their initial positions where they are fixed. This is
set up either by assigning heavyweights to thewallmolecules (Koplik et al., 1988), rescaling the
velocity of wall molecules to their initial values (Huang et al., 2004), or using fictitious springs
(Fan et al., 2002; Huang et al., 2006; Sofos et al., 2009; Branam and Micci, 2009; Kamali and
Kharazmi, 2011). Meanwhile, some studies can be found where walls are modeled differently.
An example of such investigations is the study of Ziarani andMohamad (2005), where the walls
were modeled as reflective. The investigation of wall boundary conditions is not the subject of
the present study, and the interested scientists are encouraged to refer to the literature formore
details. Aswill be described later, two rows ofmolecules that are allowed to oscillate about their
initial fixed positions represent channel walls.

2. Pressure and temperature formula in molecular dynamics

In this paper, the Non-EquilibriumMolecular Dynamics (NEMD) technique is used to simulate
flow of Argon. In molecular dynamics, positions of molecules are determined using Newton’s
second law. Intermolecular forces are calculated using Lennard-Jones 12-6 potential equation
(Rapaport, 2004) with a small correction by truncating at the cut-off distance and shifting to
eliminate the discontinuity (Leach, 1999)

V (rij)=







4ε
[( σ

rij

)12
−

( σ

rij

)6
−

(σ

rc

)12
+
(σ

rc

)6]

for rij ¬ rc

0 for rij >rc

(2.1)



Analysis of controlled molecular dynamic flow... 1143

where rij is the intermolecular distance between atoms i and j, ε is the energy scale.We choose a
cut-off distance of rc =2.5σ herein, which has commonly been used in other studies to simulate
Argon flow (Koplik et al., 1989; Priezjev, 2007; Sofos et al., 2009; Bao et al., 2015). σ is the finite
distance atwhich the intermolecular potential is zero, and it is equal tomolecular diameter in the
hard sphere model. In this paper, wall and fluid molecules are formed from the same material.
Therefore, all interaction between wall molecules, fluidmolecules or wall-fluidmolecules follows
the above equation with Argon constants.
In amacroscopic analysis, pressure is calculated solely based on the virial equation of state.

In a microscopic system however, pressure is calculated from the virial equation of state on all
atoms plus the summation of intermolecular forces multiplied by the corresponding distances
between them. For pair-additive force fields, the fluid pressureP can be estimated through the
virial equation of state given below (Allen and Tildesly, 1987)

P =
NkBT

V
+
1

3V

N
∑

i=1

N
∑

j>i

rijfij (2.2)

whereN is the number of molecules in the computational domain, kB is Boltzmann’s constant,
and V and T are volume and temperature of the computational domain, respectively. rij is the
distance between two interacting molecules i and j, and fij is the force of molecule j acting on
molecule i. The first term in Eq. (2.2) is called the kinetic term and contains temperature of the
computational domain. In a microscopic view, temperature is defined as

T =
m

3NkB

N
∑

n=1

3
∑

i=1

(Vni−V i)
2 (2.3)

where n denotes the molecule number in the domain, i=1,2,3 denotes the x, y and z compo-
nents of the atomic velocityVn, V i is the i-th component of themean flow velocityV, andm is
the atom mass. Mean velocity can be obtained from each of the several averaging methods like
SAM(Tysanner andGarcia, 2004), CAM(Tysanner andGarcia, 2005) or SMC(Karimian et al.,
2011) which reduces statistical errors in its calculation. In the present study, the SMC method
is used to calculate mean velocity. The velocity difference within parentheses that represents
instantaneous velocity due to thermal fluctuations is called virial velocity.
The pressure calculation method defined by Eqs. (2.2) and (2.3) has been implemented

in an in-house MD code developed in FORTRAN to extract pressure values from the results
of MD simulations. This code has been used in authors’ previous studies, where its accuracy
in calculation of macroscopic properties was demonstrated (Karimian et al., 2011; Karimian
and Namvar, 2012; Namvar andKarimian, 2012; Karimian and Izadi, 2013; Hasheminasab and
Karimian, 2015; Najafi andKarimian, 2016).
It should be noted that in regions near the boundary, all of the molecules including those

outside domain boundary shouldbe taken into accountwithin their cut-off region.This is crucial
for correct calculation of pressure (Najafi and Karimian, 2016).

3. Controlling flow with the temperature gradient in a channel with non-equal

inlet and outlet cross-sectional area

In (Najafi andKarimian, 2016), it was shown that uniform flow of molecules can be created by
enforcing two macroscopic properties such as temperature and mean flow velocity at a control
point; e.g. at a bin. Najafi and Karimian (2016) also showed that flow of molecules with a
gradient of flow properties, such as temperature, can be created if one flow property is defined
on an additional control point aswell. In fact, theywere able to set up a desired periodic flow. In



1144 H.R. Najafi, S.M.H. Karimian

thepresent study,wewould like to apply the approach ofNajafi andKarimian to set upadesired
flow field in a channel with non-equal inlet and outlet cross-sectional areas. To investigate this
idea, the following test case is defined.

Geometry of thenanochannel consideredhere is anon-symmetric converging-diverging chan-
nel with length of 126.525Å, inlet cross-sectional area of 36.15Å×36.15Å, outlet cross-sectional
area of 29.05Å×36.15Å, and ratio of throat to inlet cross-sectional areas of 0.6, as shown in
Fig. 1. Seven bins with equal lengths of 18.075Å are located along the channel. Note that the
bin length cannot be very small. In this case, a sufficient number of molecules will not exist
in each bin and, therefore, fluctuations will appear in the calculated properties. On the other
hand, although smooth results are obtained in larger binswith amoderate number ofmolecules,
they cannot represent local properties accurately. More details about the effect of bin size on
sampling and averaging can be found in (Karimian and Izadi, 2013).Wall molecules in bin 1 are
arranged in a straight line. From the beginning of bin 2 to the end of bin 7, the wall molecules
are located on a half-sinusoidal profile. The throat of this channel is located at the boundary of
bins 5 and 6 that is 90.375Å from the inlet. Thewall profile of bins 6 and 7 is in symmetrywith
that of bins 5 and 4. This non-symmetric channel is in fact the first part of a fully symmetric
converging-diverging channel with a length of 180.75Å. As mentioned above, the desired flow
field can be developed in this channel by enforcing temperature andmean velocity at bin 2, and
temperature at bin 7.

Fig. 1. Non-symmetric channel with 7 bins along the x direction and different inlet and outlet
cross-sectional areas

In order to verify our results, we need to know the solution ofmolecular flow in this channel
for the desired values of mean velocity and temperatures at these two bins. For this purpo-
se, we can assume that the above-defined non-symmetric channel is a part of the symmetric
converging-diverging channel for which its solution is available. In this case, the solution of the
non-symmetric channel would be equal to the solution of this symmetric channel in its first
126.525Å provided that values of mean velocity and temperatures extracted from the solution
of the symmetric converging-diverging channel are enforced at bins 2 and 7.

Therefore, we first solvemolecular flow in a 180.75Å length symmetric converging-diverging
channel, which has the exact wall profile of our test case in its first 126.525Å. A schematic of
this channelwith its initial configuration ofmolecules is shown inFig. 2.This symmetric channel
is solved with the periodic boundary condition, and has exactly the same 7 bins of our test case
plus 3 extra bins in its 54.225Å extended part. Obviously, this symmetric converging-diverging
channel has an inlet/outlet cross-sectional area of 36.15Å×36.15Å, and the ratio of throat to
inlet cross-sectional areas of 0.6.

This channel contains mono-atomic molecules of Argon gas with an atomic diameter of
σ=3.4Å and energy parameter of ε=1.67·10−21J.Argonmolecules are initially arranged in a
FCC lattice form. Solution starts from the initial temperature of 120K and zero mean velocity.
The periodic boundary condition is applied in the x and z directions. In the y direction, the
channel is restricted at each side by two rows of fluidmolecules as channel walls.Wall molecules
are allowed to oscillate about their initial fixed positions. Temperature of the wall molecules



Analysis of controlled molecular dynamic flow... 1145

Fig. 2. Schematic of molecule positions in the symmetric converging-diverging channel with 1062 fluid
molecules and 1092 wall molecules

is maintained at 120K by scaling of the molecule velocities. This channel contains 1062 flow
molecules (N) and 1092wallmolecules (Nw). As seen inFig. 2, the flowmolecules are initialized
over the channel in a rectangular cube. Each wall has 546 molecules, arranged in a FCC lattice
form inahalf-sinusoidal profile.Thus,molecular densities of thewall andfluidare 2237.79kg/m3

and 1108.18kg/m3, respectively. Low amplitude oscillations of the wall molecules, managed by
a spring stiffness coefficient of 5000N/m, decrease the possibility of molecule escape from the
domain boundaries.

The equation of motion is integrated over time with time steps of 10−16 s using the Verlet
scheme (Verlet, 1967). Note that the time step value should be chosen appropriately. Very
small time steps increase computational cost and solution time, and large time steps may cause
solution divergence. The solution domain is divided into 10 equal bins. The solution is preceded
for 15000000 time steps to reach equilibrium.

During the solution, temperature of 300K and mean velocity of 20m/s are continuously
enforced in the second bin. At each time step, mean velocity in this bin is modified by adding
the difference between desired and calculated mean velocities to the velocity of each molecule
inside the bin. In a similar fashion, temperature in the bin is updated by rescaling the virial
velocity of molecules at each time step. The SMC averaging method (Karimian et al., 2011) is
applied to calculate pressure and velocity at each time step using Eq. (2.2). Pressure,mean flow
velocity, molecular density and temperature variations along the channel are shown in Fig. 3.
Convergence rates of velocity and temperature in the bins are shown in Fig. 4.

Fig. 3. Variations of flow properties along the symmetric converging-diverging channel generated by
enforcement of temperature andmean velocity in bin 2, (a) velocity, (b) pressure and (c) molecular

density (d) temperature

It can be seen in Fig. 3a that the mean flow velocity increases through the channel from
the inlet to throat and, after that, decreases towards the end of the channel. The throat of this
channel is located between bins 5 and 6. Velocity enforcement at bin 2 results in molecules to
be driven towards the flow direction and leave the bin. Therefore, this bin would have the least
molecular density. Higher concentration ofmolecules in bin 3 leads to lower velocities and higher
pressures at this bin.



1146 H.R. Najafi, S.M.H. Karimian

Fig. 4. Convergence rates of: (a) velocity and (b) temperature, along the symmetric
converging-diverging channel generated by enforcement of temperature andmean velocity in bin 2

As shown in Fig. 3c, the highest value ofmolecular density occurs in bin 4 before the throat.
Temperature decreases from the control bin through the end of the channel, and the minimum
temperature occurs in bin 8.The temperature setting at the control bin influences upstreambins
at lower flowvelocities. So the temperature rise could be observed in bins 1, 9 and 10. According
to Eq. (2.2), the pressure increases with temperature only through the first term while the
effect of molecular density appears in both terms, giving the temperature the dominant role in
determining the pressure value and its variation. Therefore, the highest pressure has occurred
at bin 3 near the highest temperature bin, as can be seen in Fig. 3b. It should be noted that
the flow energy decreases through the channel by restricting and scaling temperature of wall
molecules to a fixed value.

As mentioned before, these results will be used to verify the solution of our test case with
non-equal inlet and outlet cross-sectional areas, defined in Fig. 1. This non-symmetric channel
includes 7 bins. The equation of motion is integrated over time with time steps of ∆t= 10−16

second using the Verlet scheme. The solution proceeds for 15000000 time steps to reach equili-
brium. Similar to the symmetric channel flow, temperature of 300Kandmean velocity of 20m/s
are continuously enforced in the second bin. In addition, temperature of 234K is enforced at
the outlet of the channel in bin 7. This temperature has been extracted from the solution of the
symmetric converging-diverging channel introduced in Fig. 3. There are 756 wall molecules and
803 fluidmolecules between bins 1 to 7which are equal to those in the same bins of the symme-
tric converging-diverging channel at the equilibrium. The number of molecules is not changed
during the solution. This means that for each molecule leaving the end of the channel, either
it should return back into the channel at the inlet or one molecule should be inserted into the
domain at the inlet. Insertion of molecules into the domain creates drastic intermolecular forces
and change the internal energy in the domain. Such instantaneous forces sometimes cause une-
xpecteddisorder in themotion ofmolecules, whichnormally results in divergence of the solution.
Somemethods have been presented to perform the insertion process with minimum changes in
the domain internal energy (Sun and Ebner, 1992; Li et al., 1998). In those methods however,
inserted molecules may overlap with other molecules and, therefore, due to high internal forces
between them, one or both of themmay be ejected out of the domain.

Due to difficulties associated with insertion of molecules at the inlet, it has been decided
to use a periodic boundary condition at the inlet and outlet of the channel. In this case, every
molecule that leaves the solution domain at the outlet will be returned back into the solution
domain at the inlet without any difficulty. However, the problem is that cross-sectional areas of
the inlet and outlet are not equal to each other. To resolve this problem, the molecules leaving



Analysis of controlled molecular dynamic flow... 1147

the outlet of the non-symmetric channel will be sent into an adaption zone which expands
the flow cross-sectional area from the outlet of the non-symmetric channel to an area equal to
the inlet cross-sectional area of the channel, i.e. from area of 29.05Å×36.15Å, to an area of
36.15Å×36.15Å. Note that the periodic boundary condition is applied between the outlet of
the non-symmetric channel and the inlet of the adaption zone, and also between the outlet of the
adaption zone and the inlet of the non-symmetric channel. Different profiles can be designed for
thewall geometry in this zone; this will be discussed later in this Section. Butwe know that the
best option is to select a wall profile similar to the part of the symmetric converging-diverging
channel that has been cut off from it, i.e. a part of the channel with bins 8 to 10, as shown
in Fig. 5. This adaption zone, denoted by 1st, is filled with 259 Argon molecules to bring the
number ofmono-atomic Argonmolecules to 1062 in the solution domain. Results obtained from
the numerical solution of molecular flow in the non-symmetric channel are compared with the
results of the symmetric converging-diverging channel in Fig. 6.

Fig. 5. 1st adaptation zone with 3 bins along the x direction

Fig. 6. Variations of flow properties along the non-symmetric channel with temperature of 300K and
mean velocity of 20m/s at bin 2 and temperature of 234K at bin 7 in comparison to those of the

symmetric channel, (a) mean flow velocity, (b) pressure, and (c) temperature

As shown in Fig. 6, the agreement between results obtained from the solutions of symmetric
and non-symmetric channels is excellent, except for the pressure variation. Note that since
temperature of bin 7 is set up to 234K by rescaling of the Virial velocity of molecules in this
bin at each time step, the arrangement of molecules in this bin will be changed in comparison
to that of the symmetric converging-diverging channel. Due to this rearrangement of molecules
in bin 7, temperature variation is slightly underestimated.As noted before, temperature has the
main role in the determination of pressure. Therefore, rearrangement of molecules in bin 7 has
affected pressure through both terms of Eq. (2.2). As a result, a small difference between the
results of pressure depicted in Fig. 6b is inevitable. Therefore, we have been able to generate
specific flow in a channel with non-equal inlet and outlet cross-sectional areas by enforcement
of temperature and mean velocity at one bin, e.g. bin 2 here, and temperature at another bin,
e.g. bin 7 here.
In order to produce extra validation data for comparison purposes, molecular flow of Argon

in the symmetric converging-diverging channel of Fig. 2 has been solved for two additional



1148 H.R. Najafi, S.M.H. Karimian

conditions at bin 2: (a) temperature of 300K andmean velocity of 50m/s, and (b) temperature
of 300Kandmeanvelocity of 100m/s.Other settings of these twonumerical solutions are exactly
the same as those of the first solution with temperature of 300K and mean velocity of 20m/s.
Solutions of cases (a) and (b) indicate that temperature of bin 7 in these two cases reaches
263K and 279K, respectively. Now, once again molecular flow is solved in the non-symmetric
channel, having enforced temperature and mean velocity at bin 2, and extracted temperature
at bin 7. All other settings of the solution are as before. Comparison of pressure, mean flow
velocity, molecular density and temperature variations along the non-symmetric channel with
those of the symmetric converging-diverging channel are shown in Fig. 7.

Fig. 7. Comparison of velocity, pressure, molecular density and temperature along the non-symmetric
channel with those of the symmetric converging-diverging channel for two cases of (a) temperature of
300K andmean velocities of 50m/s at bin 2 and temperatures of 263K at bin 7, and (b) temperature

of 300K andmean velocities of 100m/s at bin 2 and temperatures of 279K at bin 7

Variation of the mean flow velocity along the channel shows a slightly different trend at
100m/s in comparison to 50m/s case as seen in Fig. 7a. Due to the accumulation of molecules
in the throat region, the mean velocity of molecules decreases in the throat upstream for both
cases. Downstream the throat however, molecules gain velocity as a result of less concentration
in this region. It could be seen inFig. 7b that pressure increaseswithmolecular density upstream
the throat. In this region, higher temperature has also caused the pressure to increase, according
to Eq. (2.2). The agreement between the results of the symmetric and non-symmetric channels
is more realized at higher mean flow velocities. Therefore, it is clear that the presented method
is more accurate at higher velocities of mean flow.
Again, comparison of the results in Fig. 7 proved that the presented approach has been able

to generate specific flow in the channel with non-equal inlet and outlet cross-sectional areas by
enforcement of temperature andmean velocity at one bin, and temperature at another bin if a
proper adaption zone is used. In order to generalize this approach, it is needed to investigate
the effect of adaption-zone geometry on the results. It is significant to demonstrate that the
presented approach ismainly independent of the shape of the adaption zone. The adaption zone
is characterized by its length, its wall slopes at the inlet and outlet, and cross-sectional areas at
the inlet and outlet. For any channel, the latter one is fixed, i.e. from inlet of 29.05Å×36.15Å
to outlet of 36.15Å×36.15Å. Therefore, we define test cases to investigate the effect of other
two parameters, i.e. length and slope of the adaption zone.
First, two different adaption zones shown in Fig. 8 are investigated. Test case (b) where

mean flow velocity of 100m/s is enforced at bin 2 is solved with these adaption zones. The
adaptation zones proposed in Figs. 8a and 8b as 2nd and 3rd geometries include three and two
bins, respectively; i.e. equal to 54.225Å and 36.15Å in length, respectively. In both cases, the
slope of the wall profile is 0.196 in the first bin, and is zero in the other bins.
Thenumber of fluidmolecules in the 2ndand3rdadaptation zones ofFigs. 8a andFig. 8b are

equal to 259 and 153, respectively. The number of molecules in the geometries has been selected



Analysis of controlled molecular dynamic flow... 1149

Fig. 8. Modified adaptation zones: (a) 2nd geometry (b) 3rd geometry

based on the average density of molecules in the non-symmetric channel. Therefore, the added
number of molecules is proportional to the adaptation zone volume. The results obtained from
these two solutions are compared with those obtained from the symmetric converging-diverging
channel and from the non-symmetric channel with the 1st adaption zone in Fig. 9.

Fig. 9. Comparison of velocity, pressure and temperature variations of flow variables along the
non-symmetric channel generated by temperature of 300K andmean velocity of 100m/s at bin 2 and
temperatures of 279K at bin 7 with 1st, 2nd, and 3rd adaption zones, with those of the symmetric

channel

In comparison to the results of the symmetric channel and those obtained using the 1st
adaption zone, both adaption zones have predicted flow properties very well, however, between
the 2nd and 3rd adaptation zones, the 2nd adaption zone leads to a better estimation of flow
properties. We believe that this is due to the facts that 1) the slope of the 2nd adaption zone
from the inlet to outlet is smoother than that of the 3rd adaption zone, and 2) length of the 2nd
adaption zone is longer. But we would like to have an adaption zone with small length in order
to minimize computational cost. For this purpose, two additional adaption zones denoted by
fourth (4th) and fifth (5th), with 36.15Å in length are defined. Lengths of both adaption zones
are equal to the 3rd adaption zone. The 4th adaption zone shown in Fig. 10a has a smoother
wall slope. From the beginning of the first bin up to themiddle of the second bin, thewall slope
is constant and equal to 0.13, i.e. with an angle of 7.4◦. From themiddle of the second bin, the
slope of thewall is set equal to zero. The 5th adaption zone shown inFig. 10bwill be introduced
after analysis of the results obtained from the 4th adaption zone. The number of fluidmolecules
in both adaptation zones is equal to 259.

Having solved molecular flow in the non-symmetric channel using the 4th adaption zone
with temperature of 300K and mean velocity of 100m/s at bin 2 and temperatures of 279K
at bin 7, its results are compared with the results of the symmetric channel and those of the
non-symmetric channel with the 2nd and 3rd adaption zones in Fig. 11. Having compared with
the results of the symmetric channel, results of Fig. 11 show that using the 4th adaptation zone,
the flow properties are estimated better than when the 3rd adaptation zone is used.We believe
that this is due to the fact that the wall of the 4th adaption zone has a lower wall slope. This
once again proves that the adaption zone with smoother wall profiles are preferred. Note that



1150 H.R. Najafi, S.M.H. Karimian

Fig. 10. (a) 4th adaptation zone (b) 5th adaptation zone

Fig. 11. Comparison of velocity, pressure and temperature along the non-symmetric channel generated
by temperature of 300K andmean velocity of 100m/s at bin 2 and temperatures of 279K at bin 7 with

2nd, 3rd and 4th adaption zones, with those of the symmetric channel

based on Fig. 11, although the results obtained using the 4th adaption zone do not give a better
estimation of flow properties in comparison to the results obtained using the 2nd adaptation
zone, it costs computationally less since it is shorter in length.

Careful inspection of the 4th adaption zone reveals that the slope of this adaptation zone
at its inlet is very close to the slope of the non-symmetric channel at its outlet. Therefore, to
achievemore correct results, it is necessary tomatch the inlet slope of the adaptation zone with
the outlet slope of the non-symmetric channel.
To account for these requirements, it is proposed to generate thewall profile of the adaption

zonewith a 3rddegree (cubic) polynomialwhichwill have an inlet slope equal to the outlet slope
of the non-symmetric channel and an outlet slope equal to the inlet slope of the non-symmetric
channel. Two other constants of the cubic polynomial are determined by the definition of the
beginning and the end coordinates of the adaption zone. This profile that has length equal to
that of the 4th adaption zone (2 bins) is shown in Fig. 10b. The previous flow simulation is
repeated but using the 5th adaption zone. The results of this simulation including variation of
velocity, pressure and temperature along the channel are shown in Fig. 12.

Fig. 12. Comparison between variations of flow variables along the non-symmetric channel generated by
temperature of 300K andmean velocity of 100m/s at bin 2 and temperatures of 279K at bin 7 with

4th and 5th adaptation zones, (a) velocity, (b) pressure and (c) temperature



Analysis of controlled molecular dynamic flow... 1151

As seen in Fig. 12, the results obtained using the 4th and 5th adaptation zones are very
similar to each other. Especially, accurate results have been achieved at the end of the nozzle.

To show that this method does not depend on the number of molecules, the previous flow
simulations are repeated with 200more and fewer molecules in the symmetric nanochannel and
the same trend is observed with the adaptation zones in variation of velocity, pressure and
temperature.

4. Conclusion

In this paper, it is shown that flow can be generated and controlled in a channel with different
inlet and outlet cross sections by setting temperature at two different regions and mean flow
velocity at one of these regions. In a channel with non-equal inlet and outlet, periodic boundary
conditions cannot be implemented. Therefore, to handle this difficulty, an adaption zone is used
to computationally expand the outlet of the channel to a cross sectional area equal to its inlet.
This adaption zone should have the following specifications to result in a more accurate flow
prediction.

• Length: If its length is too short then the accuracy of estimated flow properties along the
channel with variable cross-sectional area decreases. On the other hand, if the length of
the channel is too long, then the cost of computation will raise. In the non-symmetric
converging-diverging channel, for instance, the making use of the adaptation zone with
length and wall profile similar to the part that has been cut off from it, results in a more
accurate estimation of flow properties.

• Slope at junctions: To achieve more accurate results, it is necessary to set the inlet slope
of the adaption zone equal to the outlet slope of the channel, and also to set the outlet
slope of the adaption zone equal to the inlet slope of the channel.

• Wall profile: A 3rd degree (cubic) polynomial is used for the wall profile of the adaption
zone.With this profile, slopes at the inlet and outlet, and also length of the adaption zone
can be fixed to the desired values and, therefore, guarantee a smooth wall profile. Note
that the slope of the channel inlet (and thus the adaption zone outlet) is preferred to be
zero.

References

1. Allen M.P., Tildesley D.J., 1987,Computer Simulation of Liquids, Oxford, NewYork

2. Bao F., Huang Y., Zhang Y., Lin J., 2015, Investigation of pressure-driven gas flows in nano-
scale channels using molecular dynamics simulation,Microfluid Nanofluid, 18, 1075-1084

3. Branam R.D., Micci M.M., 2009, Comparison of wall models for the molecular dynamics simu-
lation of microflows,Nanoscale and Microscale Thermophysical Engineering, 13, 1-12

4. Darbandi M., Abbassi H.R., Khaledi Alidusti R., Sabouri M., 2011,Molecular Dynamics
Simulation of Nano Channel as Nanopumps, ICNMM, Edmonton, Alberta, Canada

5. Fan X.J., Nhan P.T., Ng T.Y„ Xu D., 2002, Molecular dynamics simulation of a liquid in a
complex nano channel flow,Physics of Fluids, 14, 3, 1146-1153

6. Han M., 2008, Thermally-driven nanoscale pump by molecular dynamics simulation, Journal of
Mechanical Science and Technology, 22, 157-165

7. Hanasaki I., Nakatani A., 2006, Fluidized piston model for molecular dynamics simulations of
hydrodynamics flow,Modelling and Simulation in Materials Science and Engineering, 14, S9-S20



1152 H.R. Najafi, S.M.H. Karimian

8. HasheminasabS.M.,Karimian S.M.H., 2015,New indirectmethod for calculationof flow forces
in molecular dynamics simulation, Journal of Molecular Liquids, 206, 183-189

9. HuangC., Choi P.Y.K., NandakumarK., Kostiuk L.W., 2006,Molecular dynamics simula-
tion of a pressure-driven liquid transport process in a cylindrical nanopore using two self-adjusting
plates, Journal of Chemical Physics, 124, 234701

10. HuangC., NandakumarK., Kwok D.Y., 2004,Non-equilibrium injection flow in a nanometer
capillary channel, ICMENS’04, 374-378

11. Kamali R., Kharazmi A., 2011, Molecular dynamics simulation of surface roughness effects on
nanoscale flows, International Journal of Thermal Science, 50, 3, 226-232

12. Karimian S.M.H., Izadi S., 2013, Bin size determination for the measurement of mean flow
velocity inmolecular dynamics simulation, International Journal for Numerical Methods in Fluids,
71, 7, 930-938

13. Karimian S.M.H., Izadi S., Barati Farimani A., 2011, A study on themeasurement of mean
velocity and its convergence in molecular dynamics simulations, International Journal for Nume-
rical Methods in Fluids, 67, 12, 2130-2140

14. Karimian S.M.H., Namvar S., 2012, Implementation of SMC averagingmethod in a channeled
molecular flow of liquids and gases, Journal of Physics: Conference Series, 362, 01, 2029

15. Karniadakis G.E., Beskok A., Aluru N., 2002,Micro Flows and Nano Flows, Springer, New
York, 641-648

16. Kim B.H., Beskok A., Cagin T., 2010, Viscous heating in nanoscale shear driven liquid flows,
Microfluid Nanofluid, 9, 31-40

17. Koplik J., Banavar J.,Willemsen J., 1988,Molecular dynamics of Poiseuille flow andmoving
contact lines,Physical Review Letters, 60, 1282-1285

18. Koplik J., Banavar J.R., Willemsen J.F., 1989, Molecular dynamics of fluid flow at solid
surfaces,Physics of Fluids A, 1, 781-794

19. Leach A.R., 1999,Molecular Modeling: Principles and Applications, Longman

20. Li J., Liao D., Yip S., 1998, Coupling continuum to molecular-dynamics simulation: reflecting
particle method and the field estimator,Physical Review E, 57, 6, 7259-7267

21. Liu C., Li Z., 2010, Molecular dynamics simulation of composite nanochannels as nanopumps
driven by symmetric temperature gradients,Physical Review Letters, 105, 174501

22. Mi X.B., Chwang A. T., 2003, Molecular dynamics simulations of nanochannel flows at low
Reynolds numbers,Molecules, 8, 193-206

23. Nagayama G., Cheng P., 2004, Effects of interface wettability onmicroscale flow bymolecular
dynamics simulation, International Journal of Heat and Mass Transfer, 47, 501-513

24. Najafi H.R.,Karimian S.M.H., 2016,Analysis of pressure behavior in a temperature controlled
molecular dynamic flow, Journal of Theoretical and Applied Mechanics, 54, 3, 881-892

25. Namvar S., Karimian S.M.H., 2012, Detailed investigation on the effect of wall spring stiffness
on velocity profile in molecular dynamics simulation, Journal of Physics: Conference Series, 362,
01, 2039

26. Priezjev N.V., 2007, Effect of surface roughness on rate-dependent slip in simple fluids, Journal
of Chemical Physics, 127, 144708

27. Rapaport D.C., 2004,The Art of Molecular Dynamics Simulation, Cambridge University Press

28. Sadus R.J., 2002, Molecular Simulation of Fluids: Theory, Algorithms and Object-Orientation,
Elsevier

29. Sofos F., Karakasidas T.E., LiakopoulosA., 2009,Non-equilibriummolecular dynamics in-
vestigation of parameters affecting planar nanochannel flows,Contemporary Engineering Sciences,
2, 6, 283-298



Analysis of controlled molecular dynamic flow... 1153

30. StillingerF.H.,RahmanA., 1974, Improvedsimulationof liquidswaterbymoleculardynamics,
Journal of Chemical Physics, 60, 1545-1557

31. Sun M., Ebner C., 1992, Molecular-dynamics simulation of compressible fluid flow in two-
-dimensional channels,Physical Review A, 46, 4813

32. Travis K.P., Evans D.J., 1997, Molecular spin in a fluid undergoing Poiseille flow, Physical
Review E, 55, 1566-1572

33. Travis K.P., Gubbins K.E., 2000, Poiseuille flow of Lennard-Jones fluids in narrow slit pores,
Journal of Chemical Physics, 112, 1984-1994

34. Tysanner M.W., Garcia A.L., 2004, Measurement bias of fluid velocity in molecular simula-
tions, Journal of Computational Physics, 196, 173-183

35. Tysanner M.W., Garcia A.L., 2005, Non-equilibrium behavior of equilibrium reservoirs in
molecular simulations, International Journal of Numerical Methods in Fluids, 2050, 1-12

36. Verlet L., 1967, Computer “experiments” on classical fluids. I. Thermodynamical properties of
Lennard-Jonesmolecules,Physical Review, 159, 98

37. Zhang Z.Q., Zhang H.W., Ye H.F., 2009, Pressure-driven flow in parallel-plate nanochannels,
Applied Physics Letters, 95, 154101

38. Ziarani A.S.,MohamadA.A., 2005,Amolecular dynamics study of perturbed Poiseulle flow in
a nanochannel,Microfluid Nanofluid, 2, 12-20

Manuscript received November 30, 2016; accepted for print April 15, 2017