1 Online
1 Online
1 Online
View Export
Online Citation
AFFILIATIONS
1
Department of Mechanical Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japan
2
Department of Mechanical Engineering, Osaka City University, 3-3-138 Sugimoto, Sumiyoshi-ku,
Osaka 558-8585, Japan
3
Water Frontier Research Center (WaTUS), Research Institute for Science and Technology, Tokyo University of Science,
1-3 Kagurazaka, Shinjuku-ku, Tokyo 162-8601, Japan
a)
Author to whom correspondence should be addressed: hiroki@nnfm.mech.eng.osaka-u.ac.jp
ABSTRACT
In this work, we developed a calculation method of local stress tensor applicable to non-equilibrium molecular dynamics (NEMD) systems,
which evaluates the macroscopic momentum advection and the kinetic term of the stress in the framework of the Method-of-Plane (MoP),
in a consistent way to guarantee the mass and momentum conservation. From the relation between the macroscopic velocity distribution
function and the microscopic molecular passage across a fixed control plane, we derived a method to calculate the basic properties of the
macroscopic momentum conservation law including the density, the velocity, the momentum flux, and the two terms of the stress tensor,
i.e., the interaction and the kinetic terms, defined on a surface with a finite area. Any component of the streaming velocity can be obtained
on a control surface, which enables the separation of the kinetic momentum flux into the advection and stress terms in the framework of
MoP, and this enables strict satisfaction of the mass and momentum conservation for an arbitrary closed control volume (CV) set in NEMD
systems. We validated the present method through the extraction of the density, velocity, and stress distributions in a quasi-one-dimensional
steady-state Couette flow system and in a quasi-2D steady-state NEMD system with a moving contact line. We showed that with the present
MoP, in contrast to the volume average method, the conservation law was satisfied even for a CV set around the moving contact line, which
was located in a strongly inhomogeneous region.
© 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license
(http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0062889
equation typically including the viscosity, and the local acceleration atomic stress,31 for instance, provided as the stress/atom command
of the fluid is given by the gradient of the stress tensor and the exter- in LAMMPS package,32 often used to simply visualize the stress field
nal field to satisfy the momentum conservation. On the other hand, (see Appendix C).
the molecular motion is governed by the intermolecular interaction, In NEMD systems with a non-zero local flow with the macro-
and the stress tensor is defined through the average of the molecu- scopic velocity u ≠ 0, local u must be defined on a surface S enclosing
lar motion and interaction. With this respect, the microscopic stress a control volume V in NEMD systems so that the stress is consistent
tensor is comprised of the kinetic term and the intermolecular inter- with the macroscopic momentum conservation. In other words, if
action (or configuration) term. The main scope of this article is one wants to examine the macroscopic mass and momentum bal-
the proper description of the kinetic term of the stress in NEMD ances with setting an arbitrary CV in a complex flow system, e.g., a
systems. CV set around a contact line, local macroscopic properties must be
Before the establishment of the MD method, Irving and defined on a finite-sized area on a CV.
Kirkwood6 (IK) put forward the connection between the macro- In this paper, we show a calculation method of the MoP-
scopic conservation laws (for details, see Appendix A) and micro- based local stress tensor applicable to NEMD systems with a proper
scopic molecular motion governed by the intermolecular interaction definition of the density and macroscopic velocity consistent with
through the statistical mechanical theory based on the distribution the mass and momentum conservation. We provide the formu-
function in the phase space and derived an expression of the local lation for systems consisting of single-component mono-atomic
pointwise stress comprised of kinetic and interaction parts including fluid molecules for simplicity, while the present framework is also
Taylor series expansion of differences in delta functions to express applicable to systems of multi-component or poly-atomic fluid
the microscopic particle feature. molecules, in which the intra-molecular interaction term (for details,
After the introduction of numerical MD simulations,7,8 the cal- see Appendix B) gives rise to a difficulty for the stress definition
culation of average stress in homogeneous bulk systems was enabled and resulting flux,33 not the kinetic term in Eq. (2) discussed in the
based on the virial theorem, which, indeed, corresponds to the aver- present study. Related to this point, the non-uniqueness problem
age of the IK form integrated in space and time.9,10 Regarding the of the local stress/pressure tensor due to the inhomogeneity of the
local stress implemented for equilibrium MD (EMD) simulations liquids typically at the interface34–39 is also an issue regarding the
with a discrete time step, a pragmatic scheme to calculate the aver- interaction term. For the derivation, we introduced the velocity dis-
aged stress defined on a flat plane in quasi-one-dimensional (1D) tribution function (VDF) to give the average of physical properties
planar systems11 and a spherical curved surface12 was proposed, defined on a fixed control plane. To check its validity, we performed
Sk with its normal vector pointing to the kth Cartesian direction is The RHS denotes the integral weighted with the VDF considering an
calculated by oblique pillar of a base area Sk and a height ∣v k ∣δt with its central axis
parallel to v, which is typically assumed upon the derivation of the
1 crossing Sk i i vki equilibrium pressure in the kinetic theory of gases.
τklkin ≡ − ⟨ ∑ m vl i ⟩, (2)
Sk δt i∈fluid,δt ∣vk ∣ With the limit δt → 0 and by rewriting the average of f (x, v, t)
and ξ(x, t) in the oblique pillar by f (Sk , v, t) and ξ(Sk , t), respec-
where mi and vli denote the mass and l-component of the velocity tively, the integral with respect to xk in the RHS of Eq. (4) writes
vector v i of fluid particle i, respectively. We also denote the bin face
by Sk , hereafter. The angular brackets denote the ensemble aver- ∣vk ∣δt
crossing S lim ∫ dxk Sk f (x, v, t)ξ(x, t) = lim Sk f (Sk , v, t)ξ(Sk , t)∣vk ∣δt,
age, and the summation ∑i∈fluid,δt k is taken for every fluid particle δt→0 0 δt→0
i passing through Sk within a time interval of δt, which is equal to (5)
the time increment for the numerical integration. Considering that and it follows for Eq. (4) that
we deal with single-component fluid molecules of an identical mass
crossing Sk
vki ∞
m, we substitute mi with m, hereafter. A sign function ∣vki ∣
equal lim ⟨ ∑ mξ i ⟩ = lim Sk ∭ dv f (Sk , v, t)ξ(Sk , t)∣vk ∣δt. (6)
δt→0 δt→0 −∞
to ±1 is multiplied with the momentum transfer mvli across Sk to i∈fluid,δt
local density ρ(x, t) by From Eqs. (8) and (10), the macroscopic velocity ul results in
∞ ∞ ∞
ρ(x, t) = ∫ dvx ∫ dvy ∫ dvz f (x, v, t) crossing Sk
mvli
−∞ −∞ −∞ ⟨ ∑ ⟩
∞ ρul (Sk , t) ∣vki ∣
i∈fluid,δt
≡∭ dv f (x, v, t), (3) ul (Sk , t) = = lim crossing Sk
. (11)
−∞ ρ(Sk , t) δt→0 m
⟨ ∑ ∣vki ∣
⟩
∞ ∞ ∞
where we rewrite ∫−∞ dvx ∫−∞ dvy ∫−∞ dvz as ∭−∞ dv. Then, a
∞ i∈fluid,δt
By expanding Eq. (12), we obtain equation (for details, see Appendix A). By subtracting ρul uk (Sk , t)
from the rightmost-HS and leftmost-HS of Eq. (14), we obtain
∞
τklkin (x, t) = −∭ dv f (x, v, t)vk (vl − ul (x, t))
−∞ 1 crossing Sk mvki vli
∞ τklkin (Sk , t) − ρul uk (Sk , t) = − lim ⟨ ∑ ⟩, (16)
+ uk (x, t)∭ dv f (x, v, t)(vl − ul (x, t)) δt→0 Sk δt i∈fluid,δt ∣vki ∣
−∞
∞
= −∭ dv f (x, v, t)vk (vl − ul (x, t)) + uk ρul (x, t) meaning that the microscopic total momentum transfer in the
−∞
RHS corresponds to the stress minus the advection term in the
− uk ρul (x, t) LHS. Technically, the summation in the RHS of Eq. (16) is calcu-
∞
= −∭ dv f (x, v, t)vk (vl − ul (x, t)). (13) lated during the MD simulation, and as the post-process, the stress
−∞ τklkin (Sk , t) is obtained by adding the advection term ρul uk to the total
microscopic momentum transfer as
Hence, by substituting ξ(Sk , t) and ξ i in Eq. (7) with − vk (v l −ul )
∣vk ∣
and
vki (vli −ul )
− , respectively, we obtain τklkin (Sk , t) = [τklkin (Sk , t) − ρul uk (Sk , t)] + ρul uk (Sk , t), (17)
∣vki ∣
TABLE I. Microscopic expressions for the calculation of the corresponding macroscopic properties defined as the average on
the bin face Sk in NEMD systems. The top four properties can be directly calculated from NEMD systems through the MoP
procedure, whereas the others below are derived from the four.
crossing Sk
mvki vli
τklkin (Sk , t) − ρul uk (Sk , t) 1
− lim Sk δt ⟨ ∑ ∣vki ∣
⟩ Eq. (16)
δt→0 i∈fluid,δt
FIG. 1. (a) Quasi-1D Couette-type flow system of a Lennard-Jones liquid confined between two solid walls. (b) Distributions of density ρ and velocity ux calculated by the
proposed Method-of-Plane (MoP) and the volume average (VA). Solid and dashed lines denote the results of MoP and VA, respectively, while the two lines almost overlap
in this scale. (c) Difference between the MoP and the VA regarding density ρMoP − ρVA and velocity uxMoP − uxVA with their error bars depicted with semi-transparent areas
around the average.
integrated using the velocity Verlet algorithm, with a time step velocity distributions based on a standard VA, where the time aver-
Δt of 5 fs. age in equally divided bin volumes parallel to the solid wall with a
The periodic boundary condition was set in the x and y direc- height of 0.150 nm was calculated.
tions, and 4000 LJ particles were confined between two parallel solid Figure 1(b) shows the distributions of density ρ and macro-
walls consisting of the fcc crystal located on the bottom and top sides scopic velocity in the x direction ux calculated by the proposed MoP
of the calculation cell, which directed (001) and (001) planes nor- and standard VA as a reference. Note that these distributions by
mal to the z direction. Both had eight layers so that the possible the MoP can be calculated both on x-normal bin faces and on z-
minimum distance between the fluid particle and the solid particle normal ones as shown later, while only the distributions obtained
in the outmost layer was longer than the cutoff distance. The rel- on x-normal bin faces are shown as the MoP results here. Over-
ative positions of the solid particles in the outmost layers of each all, the MoP well reproduced the results by the VA, and the two
base crystal were fixed, and the temperature of those in the second lines almost overlap in this scale. Regarding the density distribu-
outermost layers was controlled at a control temperature of 100 K tion, except near the walls where layered structures are observed,
by using the standard Langevin thermostat.46 The system was first a bulk liquid with almost constant density was formed. Note that
equilibrated for 10 ns using the top wall as a piston with a con- the bulk density was not completely constant because the tempera-
trol pressure of 4 MPa without shear so that a quasi-1D system ture was not constant due to the viscous heat dissipation induced by
with a LJ liquid confined between fcc solid walls was achieved. After the extreme shear imposed on this system. The shear velocity pro-
the equilibration, a further relaxation run to achieve a steady shear files are linear throughout almost all the liquid part except in layered
flow was carried out for 10 ns by moving the particles in the out- structures, which can be understood by the change in local viscos-
most layers of both walls with opposite velocities of ±100 m/s in ity there. The density and velocity differences between the MoP and
the x direction, using the top wall as a piston with a control pres- the VA are shown in Fig. 1(c). The density difference was within
sure of 4 MPa. Finally, steady shear flow simulation was carried 10 kg/m3 , which is less than 1% of the bulk density, and the veloc-
out, keeping their z position constant at the average position dur- ity difference was also within 0.5 m/s, showing that the proposed
ing the second relaxation run, where the system pressure resulted MoP can extract the density and velocity distribution consistent with
in 3.61 MPa. the VA.
int kin
We tested the MoP expression (Table I) in the steady state, Figure 2(a) shows the distributions of τzz (≡ τzz + τzz ), τxx
int kin
where the local density, velocity, advection term, and stress were − ρux ux [≡ τxx + (τxx − ρux ux )], and ρux ux , where the first two were
obtained as a time average of 200 ns on a grid with x-normal bin directly obtained with simple addition based on Eqs. (B1) and (14) as
that in equilibrium systems without a macroscopic flow. As clearly Fig. 1. More concretely, the densities ρ(Sx ) and ρ(Sz ) averaged on
observed, τ xx − ρux ux including the advection and τ zz are different x-normal and z-normal bin faces Sx and Sz , respectively, were calcu-
away from the solid walls, indicating that the flow effect should be lated by Eq. (8) with setting k = x and k = z, whereas the macroscopic
removed to properly evaluate the fluid stress. Figure 2(b) displays the mass fluxes ρux (Sx ) and ρux (Sz ) were obtained by Eq. (10) with l = x
distributions of the stress component τ xx , τ zz , τ zx , and τ xz , where τ zz on Sx and Sz , respectively. With these definitions, ux (Sx ) ≡ ρuρ(S
x (Sx )
x)
and τ zx were calculated on z-normal bins, whereas the others were
obtained on x-normal bins. As explained above, the stress value τ xx on Sx and ux (Sz ) ≡ ρuρ(S x (Sz )
z)
on Sz can be obtained as in Eq. (11).
was calculated by adding ρux ux to τ xx − ρux ux , whereas the advec- Note that in the present laminar flow system with uz = 0, the cal-
tion terms for the others can be neglected considering uz = 0. In the culation of ux is practically not needed for the stress separation of
kin kin kin
bulk region sufficiently away from the walls, τ xx = τ zz and τ zx = τ xz τzz , τzx , and τxz in Eq. (17). Figure 3 shows the distributions of
are satisfied as expected from the solution of a laminar Couette flow, the (a) density ρ, (b) mass flux ρux , and (c) velocity ux defined on
and the former indicates that the stress τ xx is adequately calculated x-normal and z-normal bin faces in the system in Fig. 1, where
by the proposed MoP with the resulting value −τ xx (= −τ zz ) equal to the values averaged on each bin face of x-normal and z-normal are
the external pressure value of 3.61 MPa. The wall-tangential diago- plotted with setting the z position at the center of each bin face,
nal stress τ xx fluctuates near the walls as typically observed also in respectively, i.e., they are staggered by Δz/2. In the bulk, ρ, ρux , and
equilibrium systems,17,29 because of the layered structure of the liq- resulting ux averaged on bin faces with different normal directions
uid as displayed in the density distribution in Fig. 1(b). On the other agreed well, indicating that the separation of the stress and advec-
hand, τ zz was constant except near the walls, where the solid–liquid tion terms in Eq. (17) is possible with the velocity values properly
(SL) interaction acts as the external force on the liquid. Regarding evaluated by the proposed method. The difference seen around the
the off-diagonal components τ zx (= τ xz ), they were constant except top and the bottom is due to the layered structures around the two
just around the walls, where friction from the solid is included in the walls, i.e., the values on Sz are the average on a surface parallel to
force balance even in the laminar flow. the layered structure, whereas those on Sx are the average across the
In addition to the normal velocity component uk on the MoP layers.
plane Sk , the calculation of ul (l ≠ k) tangentially to Sk is needed for
the separation of τklkin (Sk ) − ρul uk (Sk ) in Eq. (17) to properly define
the stress in general flows with ul ≠ 0 and uk ≠ 0. Including this tan- B. Quasi-2D shear flow with solid–liquid–vapor
gential velocity, we compared the distributions of the density ρ, the contact lines
FIG. 3. Comparison of the time-averaged distributions of the (a) density ρ, (b) mass flux ρux , and (c) velocity ux averaged on x-normal and z-normal bin faces Sx and Sz ,
respectively.
those in the quasi-1D system. The periodic boundary condition was After the relaxation run, the density, velocity, and stress dis-
set in the x and y directions, and 20 000 LJ particles were confined tributions were obtained by the present MoP expression in the
between two parallel solid walls with a distance of about 10.4 nm and steady state with a time average of 500 ns on x-normal bin faces
a dimension of x × y = 39.2 × 3.92 nm2 so that the LJ fluid may form with a length of Δz = 0.149 nm and z-normal ones with a length of
two quasi-2D menisci with contact lines on the walls upon the pre- Δx = 0.150 nm.
liminary equilibration at a control temperature of T = 85 K without The middle panel of Fig. 4 shows the distributions of density
shear. The static contact angles on both top and bottom walls were ρ(Sx ) calculated on the x-normal bin faces, the velocity vector with
about 57○ . After the equilibration, further relaxation runs to achieve components calculated on each bin face corresponding the compo-
a steady shear flow with asymmetric menisci were carried out for nent direction, and a stress component τ zx (Sz ) calculated on the
10 ns by moving the particles in the outmost layers of both the walls z-normal bin faces, where those for ρ(Sx ) and τ zx (Sz ) are displayed
with opposite velocities of ±10 m/s in the x direction. only for half of the system with respect to the center of mass of
Indeed, apparent flow features as shown in Fig. 4 can be qualita- DATA AVAILABILITY
tively visualized by other methods such as atomic stress,31 as shown
in Appendix C, but the present method provides the distributions The data that support the findings of this study are available
of physical properties defined on a surface establishing a direct link from the corresponding author upon reasonable request.
with the conservation laws for the arbitrary CV, as shown in Table II,
and is generally applicable to a wide range of nanoscale systems APPENDIX A: MACROSCOPIC CONSERVATION LAWS
with a liquid flow. One of our future research targets is dynamic
The macroscopic equation of continuity given with the density
wetting,3,4,49 for which we plan to examine the mechanical balance
ρ(x, t) and velocity vector u(x, t) both as functions of position x and
exerted on the fluid around a CV set around the moving contact
time t writes
line as in Fig. 4. Through the comparison with the static case,17 this ∂ρ
would enable the analysis of advancing and receding contact angle ∭ dV ∂t = −∬ dSρuk nk , (A1)
V S
from a mechanical point of view.
satisfied for an arbitrary volume V in an enclosing surface S, where
nk is the k direction component of the outward unit normal vector
IV. CONCLUDING REMARKS n with respect to the infinitesimal surface element dS. The Einstein
In this work, we developed a calculation method of local stress notation is used with a dummy index k for the vectors. Similarly, the
tensor applicable to non-equilibrium molecular dynamics (NEMD) macroscopic momentum equation, the Navier–Stokes equation for
systems, which evaluates the macroscopic momentum advection an arbitrary volume V enclosed by S writes
and the kinetic term of the stress in the framework of the Method-
∂ρul
of-Plane (MoP), in a consistent way to guarantee the mass and the ∭ dV ∂t = −∬ dSρul uk nk + ∬ dSτkl nk + ∭ dVρFl ,
V S S V
momentum conservation. From the relation between the macro- (A2)
scopic velocity distribution function and the microscopic molecular
passage across a fixed control plane, we derived a basic equation to where the fluid stress tensor component τ kl expresses the stress in
connect the macroscopic field variable and the microscopic molec- the l direction exerted on a surface element with an outward nor-
ular variable. Based on the connection, we derived a method to mal in the k direction and F l denotes the external force per mass.
calculate the basic properties of the macroscopic momentum con- Equation (A2) means that the total momentum in V in the LHS
servation law including the density, velocity, and momentum fluxes can be changed by the momentum flux passing the surface S as
4
Atomic stress for mono-atomic molecules interacting with a pair J. J. Thalakkttor and K. Mohseni, “Role of the rate of surface dilatation
potential is defined as a pointwise stress tensor per atom given by in determining microscopic dynamic contact angle,” Phys. Fluids 32, 012111
19 35
J.-G. Weng, S. Park, J. R. Lukes, and C.-L. Tien, “Molecular dynamics inves- A. Harashima, Molecular Theory of Surface Tension, Advances in Chemical
tigation of thickness effect on liquid films,” J. Chem. Phys. 113, 5917–5923 Physics Vol. 1 (Wiley, 1958), p. 203.
(2000). 36
D. Schofield and J. R. Henderson, “Statistical mechanics of inhomogeneous
20
J. Cormier, J. M. Rickman, and T. J. Delph, “Stress calculation in atomistic fluids,” Proc. R. Soc. London, Ser. A 379, 231–246 (1982).
simulations of perfect and imperfect solids,” J. Appl. Phys. 89, 99–104 (2001). 37
J. S. Rowlinson, “Thermodynamics of inhomogeneous systems,” Pure Appl.
21
D. M. Heyes, E. R. Smith, D. Dini, and T. A. Zaki, “The equivalence between Chem. 65, 873–882 (1993).
volume averaging and method of planes definitions of the pressure tensor at a 38
B. Hafskjold and T. Ikeshoji, “Microscopic pressure tensor for hard-sphere
plane,” J. Chem. Phys. 135, 024512 (2011). fluids,” Phys. Rev. E 66, 011203 (2002).
22 39
J. Z. Yang, X. Wu, and X. Li, “A generalized Irving–Kirkwood formula for the K. Shi, Y. Shen, E. E. Santiso, and K. E. Gubbins, “Microscopic pressure tensor
calculation of stress in molecular dynamics models,” J. Chem. Phys. 137, 134104 in cylindrical geometry: Pressure of water in a carbon nanotube,” J. Chem. Theory
(2012). Comput. 16, 5548–5561 (2020).
23 40
D. M. Heyes, D. Dini, and E. R. Smith, “Equilibrium fluctuations of liquid state H. Kusudo, T. Omori, and Y. Yamaguchi, “Extraction of the equilibrium pin-
static properties in a subvolume by molecular dynamics,” J. Chem. Phys. 145, ning force on a contact line exerted from a wettability boundary of a solid sur-
104504 (2016). face through the connection between mechanical and thermodynamic routes,”
24
E. R. Smith, D. M. Heyes, and D. Dini, “Towards the Irving–Kirkwood limit of J. Chem. Phys. 151, 154501 (2019).
the mechanical stress tensor,” J. Chem. Phys. 146, 224109 (2017). 41
Y. Imaizumi, T. Omori, H. Kusudo, C. Bistafa, and Y. Yamaguchi, “Wilhelmy
25
E. R. Smith and C. Braga, “Hydrodynamics across a fluctuating interface,” equation revisited: A lightweight method to measure liquid-vapor, solid-
J. Chem. Phys. 153, 134705 (2020). liquid, and solid-vapor interfacial tensions from a single molecular dynamics
26 simulation,” J. Chem. Phys. 153, 034701 (2020).
K. Shi, E. E. Santiso, and K. E. Gubbins, “Can we define a unique microscopic
42
pressure in inhomogeneous fluids?,” J. Chem. Phys. 154, 084502 (2021). H. C. Andersen, “Rattle: A ‘velocity’ version of the shake algorithm for molecu-
27 lar dynamics calculations,” J. Comput. Phys. 52, 24–34 (1983).
G. Bakker, Kapillarität und Oberflächenspannung (Wien-Harms, 1928), Vol. 6.
28 43
J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Dover, 1982). P. J. Daivis, K. P. Travis, and B. D. Todd, “A technique for the calcula-
29 tion of mass, energy, and momentum densities at planes in molecular dynamics
S. Nishida, D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima, and
H. Fujimura, “Molecular dynamics analysis of multiphase interfaces based on in simulations,” J. Chem. Phys. 104, 9651–9653 (1996).
44
situ extraction of the pressure distribution of a liquid droplet on a solid surface,” K. Ogawa, H. Oga, H. Kusudo, Y. Yamaguchi, T. Omori, S. Merabia, and L. Joly,
J. Chem. Phys. 140, 074707 (2014). “Large effect of lateral box size in molecular dynamics simulations of liquid-solid
30 friction,” Phys. Rev. E 100, 023101 (2019).
D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima, and H.
45
Fujimura, “Molecular dynamics analysis on wetting and interfacial properties of H. Oga, Y. Yamaguchi, T. Omori, S. Merabia, and L. Joly, “Green-Kubo mea-
water-alcohol mixture droplets on a solid surface,” J. Chem. Phys. 140, 034505 surement of liquid-solid friction in finite-size systems,” J. Chem. Phys. 151, 054502
(2014). (2019).