Synchronized molecular dynamics simulation via macroscopic heat and momentum transfer: an application to polymer lubrication

0.28 MB·

by  Shugo Yasuda

Synchronized molecular dynamics simulation via macroscopic heat and momentum
transfer: an application to polymer lubrication
Shugo Yasuda1 and Ryoichi Yamamoto2
1Graduate School of Simulation Studies, University of Hyogo, Kobe 650-0047, Japan,
2Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan
A synchronized molecular dynamics simulation via macroscopic heat and momentum transfer
is proposed to model the non-isothermal flow behaviors of complex fluids. In this method, the
molecular dynamics simulations are assigned to small fluid elements to calculate the local stresses
and temperatures and are synchronized at certain time intervals to satisfy the macroscopic heat-
4 and momentum-transport equations. This method is applied to the lubrication of a polymeric
1 liquid composed of short chains of ten beads between parallel plates. The rheological properties
0 and conformation of the polymer chains coupled with local viscous heating are investigated with a
2 non-dimensionalparameter,theNahme-Griffithnumber,whichisdefinedastheratiooftheviscous
heating tothe thermal conduction at thecharacteristic temperature required to sufficiently change
l
u the viscosity. The present simulation demonstrates that strong shear thinning and a transitional
J behavioroftheconformationofthepolymerchainsareexhibitedwitharapidtemperaturerisewhen
6 theNahme-Griffithnumberexceedsunity. Theresultsalsoclarifythatthereentranttransitionofthe
1 linear stress-optical relation occurs for large shear stresses dueto the coupling of the conformation
of polymer chains and dueto theheat generation undershear flows.

]
t
f
o I. INTRODUCTION for the efficient numericalcomputation ofproblems with
s multiscale characteristics. HMM has been applied to
.
t To predict the transportphenomena of complex fluids various problems, such as the simple polymeric flow[9],
a
m causedbythe coupledheatandmomentumtransferpro- coarseningofabinarypolymerblend[10]andthechannel
cesses is challenging from both a scientific and engineer- flow ofa simple Lennard-Jonesliquid[11]. The equation-
-
d ingpointofview. Moleculardynamics(MD)simulations free multiscale computation proposed by Kevrekids et
n are often used to predict material properties (e.g., the al. is also based on a similar idea and has been ap-
o rheological, thermal, and electrical properties), in which plied to various problems.[12, 13] De et al. proposed the
c the simulation is performed for a very small piece of the scale-bridgingmethod,whichcancorrectlyreproducethe
[
material under a certain ideal environment.[1, 2] How- memory effect of a polymeric liquid, and demonstrated

3 ever, in actual engineering and biological systems, the the non-linear viscoelastic behavior of a polymeric liq-
v macroscopicfeaturesofcomplexfluidsarehighlyaffected uid in slab and cylindrical geometries.[14, 15] The mul-
4 by the spatial heterogeneity caused by the macroscopic tiscale simulationfor polymeric flows with the advection
4
transport phenomenon involved in the boundary condi- of memory in two and three dimensions was developed
2
1 tions. A typical example is the generation of heat in by Murashima and Taniguchi.[16–18] We have also de-
. lubrication systems.[3] To predict such complicated be- veloped a multiscale simulation of MD and computa-
1 haviorincomplexfluids,theentiresystem,includingthe tionalfluiddynamics(CFD).Themultiscalemethodwas
0
boundary conditions, must be considered on the basis of firstdeveloped for simple fluids[19] andsubsequently ex-
4
1 an appropriate molecular model. In principle, full MD tended to polymeric liquids with the memory effect.[20–
: simulations of the entire system can meet these require- 23]
v
ments. However,it is difficult in practice to perform full However, a multiscale simulation for the coupled heat
i
X MD simulations on the macroscopic scale, which is com- and momentum transfer of complex fluids has yet to be
r mon in actual engineering systems and is far beyond the proposed. Forcomplexfluids,itisusuallydifficulttode-
a
molecular size. Multiscale modeling is a promising can- scribe the heat generation coupled with the momentum
didate to address this type of problem. transport using only macroscopic quantities. The spa-
The multiscale simulation for the flow behaviors of tial variation in the temperature also becomes notable
complexfluidswasfirstadvancedintheCONNFFESSIT atthe macroscopicscaledue to localviscousheating un-
approach for polymeric liquids by Laso and O¨ttinger[4– der shear flow. Thus, multiscale modeling is important
6], where the local stress in the fluid solver is calcu- for the coupled heat and momentum transfer of complex
lated using the microscopic simulation instead of using fluids. In the present paper, we propose a multiscale
any constitutive relations. The GENERIC approach is simulation,termedthesynchronizedmoleculardynamics
also presentedfor the non-isothermalpolymeric flows by simulation (SMD), for the coupled heat and momentum
making important corrections and clarifications to the transfer in complex fluids by extending the multiscale
CONNFFESSITscheme.[7]Thestrategyexploitedinthe simulation for momentum transport and apply it to the
CONNFFESSITapproachisalsointroducedintothehet- polymer lubrication. Using this method, we investigate
erogeneous multiscale modeling (HMM), which was pro- the rheological properties and conformation of polymer
posed by E and Enquist,[8] as a general methodology chainsforthethermohydrodynamiclubricationofapoly-
2
meric liquid composed of short chains in a gap between The macroscopic behavior of the polymeric liquid is
parallel plates, in which the width of the gap is suffi- described by the following transport equations:
cientlylargecomparedtothe characteristiclengthofthe
∂v ∂σ
flow behaviors, e.g., the length of the viscous boundary ρ x = xy, (3a)
0
layer, such that the macroscopic quantities, e.g., the ve- ∂t ∂y
locity, stress and temperature, become spatially hetero- ∂e ∂2T
ρ =σ γ˙ −λ , (3b)
geneous. 0∂t xy ∂y2
The full MD simulations for confined short-chain
molecules in slab geometry were previously studied in where vα is the velocity, σαβ is the stress tensor, e is
Refs. [24] and [25]. The viscous heating of simple liq- the internal energy per unit mass, and γ˙ is the shear
uids in the same geometry was also studied using a full rate, i.e., γ˙ = ∂vx/∂y. Hereafter, the subscripts α,
MD simulation in Refs. [26] and [27]. However, these β, and γ represent the index in Cartesian coordinates,
results areexclusivelyfor the molecularlythin films, i.e., i.e., {α,β,γ}∈{x,y,z}. Here, we assume that the macro-
the slab width is approximately ten times the molecular scopic quantities are uniform in the x and z directions,
size,whereasinthepresentstudy,awidthofoverathou- ∂/∂x=∂/∂z=0, and the density of the polymeric liquid
sand times the molecular size is considered in the SMD is constant. Fourier’s law for a heat flux with a con-
simulation. stant and uniform thermal conductivity λ is also consid-
ered in Eq. (3b). Note that the thermal conductivity
In the following, we first describe the problem consid-
of polymeric liquids is anisotropic under shear flows in
ered in the present paper. The simulation method is ex-
general[29–32], and some experimental studies have re-
plained after the presentationof the problem. The SMD
ported that the linear stress-thermal relation between
simulation of polymer lubrication is performed, and the
the stress tensor and thermal conductivity holds.[33–
results are discussed; these results are mainly the rheo-
36] However, in the present study, we only consider the
logical properties and the coupling of the conformations
isotropic thermal conductivity as the first step because
of the polymer chains with heat generation under shear
the effect of shear thinning of the viscosity is thought
flows. We also present a critical analysis of the concep-
to be more crucial to viscous heating under strong shear
tual and technical issues of the SMD method. Finally, a
flows than that of the anisotropy of the thermal conduc-
short summary and a perspective on the future of SMD
tivity. Involving the anisotropic thermal conductivity in
are given.
theSMDsimulationisanimportantfuturework. Wealso
assume that the velocity and temperature of the poly-
meric liquid are the same as those of the plates at the
II. PROBLEM
boundaries, i.e., the non-slip and non-temperature-jump
boundary conditions.
We consider a polymeric liquid contained in a gap of
The effect of viscous heating is estimated using the
widthH betweenparallelplateswithaconstanttemper-
ratio of the first and second terms in Eq. (3b) to be
ature T0 (see Fig. 1(a)). The polymeric liquid is com- σ Γ˙H2/λ∆T . Here, Γ˙ is the gross shear rate of the sys-
posed of short Kremer-Grest chains[28] of ten beads, in 0 0
tem, which is defined by the ratio of the velocity of the
which all of the bead particles interact with a truncated upperplatev tothewidthofthegapH,Γ˙ =v /H,and
Lennard-Jones potential defined by w w
∆T isacharacteristictemperatureriseforthepolymeric
0
liquid. Inthepresentproblem,weconsideracharacteris-
4ǫ σ 12− σ 6 +ǫ, (r ≤21/6σ), tictemperaturenecessarytosubstantiallychangethevis-
U (r)= r r (1)
LJ
(0, h(cid:0) (cid:1) (cid:0) (cid:1) i (r ≥21/6σ), cosityofthepolymericliquid,i.e.,∆T0 =|η0/(∂η0/∂T0)|,
where η is the characteristic viscosity of the polymeric
0
andconsecutivebeadsoneachchainareconnectedbyan liquid at a temperature of T0. Thus, the Nahme-Griffith
anharmonic spring potential, number Na, defined as
UF(r)=−12kcR02ln 1− Rr 2 , (2) Na= λ|∂logσ(0ηΓ˙0H)/2∂T0|−1 (4)
" (cid:18) 0(cid:19) #
representsthe effect of viscousheating onthe changesin
with k =30ǫ/σ2 and R =1.5σ. The polymeric liquid is the rheological properties.[3, 37] Usually, in lubrication
c 0
in aquiescentstate witha uniform density ρ anda uni- systems and in high-speed processing operations with
0
form temperature T before a time t = 0. Hereafter, polymericliquids,theNahme-Griffithnumberisnotneg-
0
the y-axisisperpendiculartothe parallelplates,andthe ligibly small because of the large viscosity and the small
boundaries between the upper and lower plates and the thermal conductivity of the polymeric liquid.[37] For ex-
polymericliquidarelocatedaty =H and0,respectively. ample, when a lubrication oil in a gap with a width of 1
The upper plate starts to move in the x direction with a µm is subjected to shear deformation with a strain rate
constant shear stress σ at a time t=0, while the lower of 1×106 s−1, the Nahme-Griffith number is estimated
0
plate is at rest. to be Na & 0.1. Thus, the rheological properties of the
3
Aconstantshearstressσ0isappliedtotheupperplate.
ShortK-Gchains
are involved in
each MD cell at
aseverlyjammed
condition.
MDcellsareasignedtoeachsmallfluidelementtocalculatethe MDsimulationsaresynchronizedattimeintervalsof
localquantities. ∆tviamacroscopicheatandmomentumtransfer.
(a) Geometry of problem (b) Time evolution
FIG. 1. The schematic of thegeometry of theproblem (a) and thetime evolution (b).
lubricant in such micro devices must be significantly af- calculate δK global mesh
∆t
fectednotonlybythe largevelocitygradientbutalsoby
thetemperatureincreasecausedbylocalviscousheating.
Topredicttherheologicalpropertiesofthepolymericliq- R Tdτ δK
uid in these systems, one must consider the temperature ∆τ
variation in Eq. (3b) coupled with Eq. (3a). local MD
Tn−1 T✲nT′n Tn+1
rescaling
III. SIMULATION METHOD
FIG. 2. A schematic for the calculation of the temperature
in the SMD method. The upper side represents the progress
Inthepresentsimulation,thegapbetweentheparallel
on the global mesh system, and thelower side representsthe
plates is divided into M mesh intervals with a uniform progress at each MD cell.
width of ∆y = H/M, and the local velocities are calcu-
lated at each mesh node through the typical finite vol-
ume scheme shown in Eq. (3a). The local shear stresses
σ (y) are calculated from the local shear rates in the the SMD method. The instantaneoustemperaturesT at
xy
MD cells associated with each mesh interval using the eachMDcellandtheirintegralsoverthedurationofeach
NEMD simulation with the SLLOD algorithm. The MD MDrun,i.e., n∆t T(τ)dτ,arecalculatedateachMD
(n−1)∆t
simulations are performedin a time interval ∆t, and the
cell. The heat fluxes between neighboring MD cells δK
time integrals of the instantaneous shear stresses P in R
xy are calculated on the global mesh system (depicted on
each MD cell are used to update the local velocities at
the upper side in Fig. 2) as
the next time step in accordance with the macroscopic
momentum transport Eq. (3a), λ ∂2 n∆t
δK=− T(τ)dτ, (6)
vn(y)=vn−1(y)+ ∂ n∆t P (τ;γ˙n−1(y))dτ. (5) ρ0∂y2 Z(n−1)∆t
x x ∂y Z(n−1)∆t xy and the instantaneous kinetic energies K at each MD
cell are corrected by rescaling the molecular velocities
Here,the superscriptn representsthe time stepnumber, ′
P (τ;γ˙n−1(y)) is the instantaneous shear stress in the according to the corrected temperature T (depicted on
NxEyMD simulation with a shear rate of γ˙n−1(y), and τ the lower side in Fig. 2) via
is the temporal progress of the NEMD simulation. Note 2
that the time-step size of the MD simulation is different T′n =Tn+ δK. (7)
3
from ∆t. The local viscous heating caused by the shear
flow, i.e., the first term of Eq. (3b), is calculated in Thus, the MD simulations assigned to each fluid ele-
the NEMD simulations without the use of a thermostat ment are synchronized at time intervals of ∆t to satisfy
algorithmineachMDcell;however,ateachtimeinterval the macroscopic heat and momentum transport equa-
∆t, the instantaneous kinetic energies of the molecules tions (see Fig. 1(b)). Note that it is difficult to rewrite
per unit mass K in eachMD cell are correctedaccording Eq. (3b) as the time evolution of temperature in gen-
to the heat fluxes between neighboring MD cells. Figure eralbecausethe internalenergydepends onnotonly the
2 is a schematic of the calculation of the temperature in macroscopic variables but also on the conformations of
4
the polymer chains. In the present SMD simulation, the 0.1 2
σ0=0.08 σ0=0.08
temperature increase caused by local viscous heating is
0.08
1.5
calculated autonomously in the MD simulation and sat- σ0=0.06
0.06
isfies the macroscopic energy balance of Eq. (3b).
Hereafter, we measure the length, time, temperature σxy 0.04 σ0=0.05 T 1
anddensityinunitsofσ,τ = mσ2/ǫ,ǫ/k ,andm/σ3, 0.5 σ0=0.06
respectively. Here,k isth0eBoltzmannconBstant,andm 0.02 σ0=0.01 σ0=0.05
is the mass of the LBJ particlep. In the following simula- 00 y/H 1 00 y/σH0=0.01 1
tions, the density and thermal conductivity of the poly- (a) (b)
meric liquid are fixed to be ρ = 1 and λ=150, respec-
0 0.01
0.05
tively,andthetemperatureoftheplatesandthewidthof σ0=0.08
σ0=0.05
the gap between the plates are fixed to be T = 0.2 and 0.04
0 σ0=0.06
H =2500, respectively, whereas the shear stress applied 0.001
0.03 σ0=0.06
to the upper plate σ0 varies. At this number density ρ0 γ˙ Qxy
and this temperature T0, the conformation of the bead 0.0001 σ0=0.05 0.02 σ0=0.08
particles becomes severely jammed and results in com- 0.01
plicated rheologicalproperties.[22, 38] σ0=0.01
1e-05 0
We have carriedout the numericaltests of the present 0 y/H 1 0 y/H 1
(c) (d)
method for variouscalculationconditions by varyingthe
numberofmeshintervalsM,time interval∆t,andnum-
ber of polymer chains in each MD cell Np. The results FIG. 3. The spatial variation in the shear stress σxy(y) (a),
andacriticalanalysisofthe presentmethodaregivenin temperatureT(y)(b),shearrateγ˙(y)(c),andxy component
Sec. V.Inthepresentpaper,unlessotherwisestated,the ofthebondorientationtensorinEq. (8)(d)forvariousvalues
following calculation parameters are used: The number of σ0.
of mesh intervals M and the time interval ∆t are M=32
and ∆t=1, respectively. A total of 100 polymer chains,
behaviors of the spatial variations in the temperature
i.e., N =100,of ten beads,i.e., 1000beadparticles,are
p
and in the shear rate are different. The spatial varia-
contained in each MD cell. Thus, the mesh width and
tion in the shear rate is only enhanced in the vicinity of
side length of the MD cell are ∆y=78.125 and l =10,
MD
the plate and is rather moderate, except in the vicinity
respectively. The time-step size in the MD simulation
ofthe plate, whereasthe spatialvariationinthe temper-
∆τ is set to be ∆τ = 0.001. Thus, MD simulations are
ature is a parabolic curve throughout the region. Note
performed for 1000 time steps in the time interval ∆t,
thatinFig. 3(c),weomittheresultforσ =0.01because
i.e., ∆t=1000∆τ. 0
the amplitude of the shear rate is very small and occa-
sionallyexhibits negativevalues as aresultofthe fluctu-
ations; the spatial averageof the shear rate for σ =0.01
IV. RESULTS 0
is 7.2×10−6. In Fig. 3(d), we show the xy component
of the bond orientation tensor Q , which is defined as
xy
WeperformedSMDsimulationswithvariousvaluesfor
the shear stress applied to the plate σ0, i.e., σ0 =0.002, 1 1 Nb−1 b b
0.01, 0.03, 0.05, 0.055, 0.06, 0.07, 0.08, and 0.09, and Q = jα jβ , (8)
αβ
investigatedthebehaviorsofapolymericliquidatsteady Np Nb−1 bminbmin
chain j=1
X X
state. Inthefollowing,wepresentquantitiesaveragedfor
a long period of time at steady state, in which the shear where N is the number of polymer chains in each MD
p
stress is spatially uniform and the time derivative of the cell, N is the number of bead particles in a polymer
b
local internal energy, i.e., the left-hand side of Eq. (3b), chain, b for 1≤j ≤N −1 is the bond vector between
j b
is negligible. consecutive beads in the same chain, and b is the dis-
min
Figure3showsthespatialvariationsinlocalquantities tance at which the sum U (r)+U (r) has a minimum
LJ F
forvariousvaluesofthe shearstressappliedtotheplate, andiscalculatedtobeb ≃0.97. Thespatialvariation
min
i.e., for σ =0.01, 0.05, 0.06, and 0.08. The shear stress in Q is also small for small applied shear stresses, i.e.,
0 xy
σ is found to be spatially uniform, and this fact also σ ≤0.05, but this spatial variationbecomes notable for
xy 0
indicates that the condition necessary for the polymeric σ ≥ 0.06. In the vicinity of the plate, Q increases
0 xy
liquid to be at steady state is satisfied in the present monotonically with increasing applied shear stress σ ,
0
simulation. The spatial variation in the temperature T whereas in the middle of the plates, Q varies non-
xy
is small when the applied shear stress σ is smaller than monotonically with the applied shear stress. This com-
0
0.05, i.e., σ ≤ 0.05, whereas for σ ≥ 0.06, the spatial plicated behavior arises from the temperature variation;
0 0
variation in T becomes notable. The spatial variation the conformationofthe polymerchainsin the vicinity of
in the shear rate γ˙ also increases rapidly when the ap- the plate is mainly influenced by the shear rate because
plied shear stress σ is larger than 0.05. However, the thevariationinthetemperaturefromthereferencevalue
0
5
0.1 10 Ref. 38 as ∆T = 0.15. In Fig. 4(a), the apparent vis-
σ¯xy 0
1000 T¯ cosityiscomparedtotheviscosityofthepolymericliquid
with the uniform and constant temperature T obtained
∝Γ˙−0.7 0
η σ¯xy0.01 1 T from the NEMD simulation; this case corresponds to a
100 ∝Γ˙−1 negligible Nahme-Griffith number, i.e., Na = 0. Strong
shear thinning occurs when the Nahme-Griffith number
τR−1 Na=1 Na=1 exceeds unity, i.e., Na>1. The slope in the power law,
10 1e-06 1e-05 0.0001 0.001 0.01 0.0011e-06 1e-05 0.0001 0.001 0.010.1 i.e., the index ν inthe approximaterelationη(Γ˙)∝Γ˙−ν,
Γ˙ Γ˙ is almost unity (but never exceeds unity), i.e., ν ≃1, for
(a) (b) Γ˙ = 1×10−4 ∼ 1×10−3. A discrepancy between the
FIG. 4. The apparent viscosity η, defined as η = σ0/Γ˙, (a) apparent viscosity and that for Na=0 is thought to be
causedbythefactthatthetemperatureslightlyincreases
and thespatial averages of theshear stressand temperature,
σ¯xy andT¯,(b)asfunctionsofthegrossshearrateΓ˙. Infigure inthepresentsimulationevenatthesmallestgrossshear
(a), the asterisks “∗” show the results of the viscosity for a rate(the averagetemperatureT¯ increasesby2.5percent
uniform temperature T = 0.2, i.e., Na = 0, obtained by the of the wall temperature T ), and the temporal-spatial
0
NEMD simulations. The downward arrows in both figures fluctuations of the temperature and shear rate are also
representtheshearrateatwhichtheNahme-Griffithnumber induced in the present simulation, whereas the temper-
Na definedin Eq. (4) equals unity. ature and the shear rate are kept constant and uniform
in the NEMD simulation. Figure 4(b) shows the spa-
0.38 0.05 tialaveragesofthe localshearstressandtemperature as
Q¯xx functionsofthegrossshearrate. Wenotethatthespatial
Q¯yy averageoftheshearstressσ¯ coincideswiththeapplied
0.36 Q¯xy 0.04 shearstressσ0 ontheplatebxeycausethelocalshearstress
is almost spatially uniform at steady state, as shown in
0.03 Fig. 3(b). The average shear stress σ¯ is observed to
¯Qyy0.34 xy increase monotonically with the grossxsyhear rate. The
,xx ¯Q plateau region of the curve corresponds to the strongly
¯Q 0.02
shear-thinningregimewiththeindexν ≃1inFig. 3(a).
The average temperature increases rapidly with increas-
0.32
0.01 ing gross shear rate when the Nahme-Griffith number
Na=1 exceeds unity. The rate of increase of the average tem-
perature is lowerthan thatof the averageshear stressat
0.3 0
1e-06 1e-05 0.0001 0.001 0.01 small gross shear rates, i.e., Γ˙ . 1×10−4, but reverses
Γ˙ at large gross shear rates, i.e., Γ˙ &1×10−4.
Figure 5 shows the conformation changes in the poly-
FIG.5. Thespatialaverageofthebond-orientationtensorof
mer chains as a function of the gross shear rate. A non-
the polymer chains Q¯ as a function of the gross shear rate
αβ monotonic dependence of the conformation of the poly-
Γ˙. In addition, see thecaption in Fig. 4
mer chains onthe grossshear rate is observed. A transi-
tion occurs at the entrance of the strong shear-thinning
regime, with ν ≃1 in Fig. 4(a), i.e., Γ˙ ≃ 1×10−4. The
T0issmallthere,butatthemiddleoftheplates,thecon- polymer chains are stretched in the x direction, and the
formation of the polymer chains is greatly influenced by xy componentofthebond-orientationtensorincreasesas
thetemperaturebecausethevariationinthetemperature Γ˙ forsmallgrossshearrates,i.e.,Γ˙ .1×10−4. However,
from the reference value T0 is large there. The relations for large gross shear rates, i.e., Γ˙ & 1×10−4, the align-
among the temperature, shear rate, and conformationof ment of the polymer chains is disturbed, and the confor-
the polymer chains are also discussed in detail later. mation of the polymer chains becomes more isotropic as
In the following, we investigate the rheological prop- thegrossshearrateincreases. This transitionalbehavior
erties of the system. Hereafter, the spatial average of a is caused by the temperature variation; when the tem-
quantity a(y) is represented as a¯, i.e., a¯ = 1 Ha(y)dy. perature increases sufficiently for the large gross shear
H 0
Figure 4 shows the apparent viscosity of the system and rates, the coherent structure becomes disturbed by the
the spatial averages of the shear stress and teRmperature thermal motion of the bead particles.
as a function of the gross shear rate. The downward Figure 6 shows the stress-optical relation for the
arrows in Fig. 4 indicate the reference shear rate at presentproblem.[39–41]Forthepresentmodelpolymeric
whichthe Nahme-Griffithnumber Na defined in Eq. (4) liquid, the NEMD simulations for the isothermal shear
equals unity, where the reference temperature is calcu- flows, i.e., Na=0, showed that a universal curve in the
lated from the viscosities of the model polymeric liquid stress-optical relation holds in both the linear (Q¯ .
xy
that were obtained by the NEMD simulation at the uni- 0.05) and nonlinear (Q¯ .0.05) regimes.[38] In the fig-
xy
form and constant temperatures of T = 0.2 and 0.4 in ure, the present results are compared with those for the
6
10 0.02 0.208
C1
C2 0.207
present 0.015 C3 0.206
T =0.2 C4
T =0.4 C5 0.205
1 C6
T =1.0 vx 0.01 T 0.204
0.203
0.005 0.202
¯/T 0.1 ng˙Γ re-entrant 0.201
σ¯xy ncreasi 00 y/H 1 0.20 y/H 1
i
(a)velocity (b)temperature
0.01
FIG.7. Thecomparisonsofthevelocity(a)andtemperature
(b) profiles under the various calculation conditions for the
applied shear stress σ0 =0.01.
0.001
0.0001 0.001 0.01 0.1 1
2
Q¯xy 78 CC12 1.8
C3 1.6
FIG. 6. The stress-optical relation σ¯xy/T¯ vs. Q¯xy. The 56 CC45 11..24
squares (cid:3) indicate the present results, and the symbols +, vx4 T 1
×, and ∗ indicate the results obtained in Ref. 38 at uniform 3 0.8
temperatures. 2 0.6
1 0.4
0 0.2
0 y/H 1 0 y/H 1
case Na=0. The leftmost square symbol in the present
(a)velocity (b)temperature
study represents the result for the smallest gross shear
rate, i.e., Γ˙ ≃ 1×10−6. It is observed that σ¯ /T¯ in- FIG.8. Thecomparisonsofthevelocity(a)andtemperature
xy
creaseswith Q¯ , while the grossshearrate increasesup (b) profiles under the various calculation conditions for the
xy
to Γ˙ . 1×10−4 because the shear stress increases more applied shear stress σ0 =0.08.
rapidly with the gross shear rate than does the temper-
ature (see Fig. 4(b)); moreover, Q¯ also increases with
xy
Gavis and Laurence are also presented.[42]
thegrossshearrateinthisregime(seeFig. 5). However,
forlargegrossshearrates,i.e.,Γ˙ &1×10−4,thetemper-
ature increases more rapidly than does the shear stress,
TABLE I.Calculation Conditions
andthexy componentofthebond-orientationtensorde-
creaseswiththegrossshearrate. Interestingly,thelinear M ∆t Np
C1 32 1.0 100
stress-opticalrelationisrecoveredforshearstresseslarger
C2 64 1.0 100
than that for the transitional behavior of the conforma-
C3 16 1.0 100
tion tensor, although the temperature, shear stress, and C4 32 0.5 100
conformation of the polymer chains exhibit very compli- C5 32 5.0 100
catednonlinearbehavior. Thisreentranttransitionofthe C6 32 1.0 1000
linear stress-optical relation for large shear stresses can
never be reproduced in the NEMD simulations using a
Figures 7 and 8 show the comparisons of the veloc-
thermostatbecauseboth the shearstress σ¯ andthe xy
xy ity and temperature profiles obtained under the differ-
componentofthe bondorientationtensorQ¯ monoton-
xy ent calculation conditions for the applied share stresses
icallyincreasewiththeshearrateataconstanttempera-
σ = 0.01 and σ = 0.08, respectively. The comparisons
ture T¯, but Q¯ saturates to a limiting value Q¯ ∼0.1; 0 0
xy xy between C1, C2, and C3 show the effect of changing the
therefore, the nonlinear stress-optical relation forms, as
mesh interval ∆x, and the comparisonsbetween C1, C4,
shown in Fig. 6, for the case Na=0.
and C5 show the effect of changing the time interval ∆t.
Note that the time-step size of the MD simulation ∆τ
is fixed to be ∆τ = 0.001. Thus, the MD simulations
V. CRITICAL ANALYSIS OF THE METHOD
are performed for 500 and 5000 time steps for C4 and
C5, respectively, in each time interval ∆t. The compari-
In this section, numerical tests for various calculation son between C1 and C6 shows the effect of changing the
conditions are implemented by varying the number of numberofparticlesineachMDcell,i.e.,theeffectofthe
meshintervalsM (i.e.,∆x=H/M),timeinterval∆t,and noise intensity arising from the MD simulations.
numberofpolymerchainsineachMDcellN ,whichare It is observed in Fig. 7 and in 8 that for σ = 0.01,
p 0
givenin Table I. The comparisonsbetweenthe results of the velocity profiles for C2 and C5 significantly deviate
the present method and the analytic solutions given by from those under other calculation conditions, whereas
7
for σ0 = 0.08, the deviations under different calculation 1 C1
conditions arenotnotable. The deviationsofC2 andC5 C6
σ0=0.05 analytic
for σ = 0.01 are considered to be a result of the nu-
0 0.1
merical stability condition for the momentum transport
equation Eq. (3a), which is written as ∆t <∆x2/2η for /T0 σ0=0.03
a Newtonian fluid with a constant viscosity η. The ex- T)0 0.01
−
act stability condition for the SMD method is unknown T σ0=0.01
(
because in the SMD simulations, the fluctuations are in-
0.001
volved in local stresses, and the local viscosities are also
autonomously modulated according to the local flow ve-
locities. In the present numerical tests, only the calcu- 0.0001
lation conditions C2 and C5 do not satisfy the condition 0 y/H 1
∆t < ∆x2/2η, where η is the proper viscosity obtained
FIG.9. Thecomparisonsoftemperatureprofilesbetweenthe
underothercalculationconditions. Asimilarsituationis
SMDsimulations andtheanalyticsolutions forvariousshear
alsofoundinthepreviousstudyinRef. [21],inwhichnu-
stresses.
merical tests for the creep motion of the model polymer
melt with a uniform temperature are performed. In the
previous study, we found that the multiscale simulation temperature is written as[42, 43]
in which the time-step size ∆t is larger than the viscous
diffusion time ∆x2/η, i.e., ∆t > ∆x2/η, reproduces the T(yˆ)=T0+∆T0×
velocity profile with a pseudo viscosity that is smaller
N˜a N˜a
than the proper viscosity. These facts indicate that the log 1+ sech2 arcsinh (2yˆ−1) ,
time interval ∆t must be smaller than the viscous diffu-  8 !  s 8  
sion time ∆x2/η, within which the viscous force propa-  
  (9)
gatesforthemeshinterval∆x,toobtainthepropersolu-  
tionsintheSMDsimulations. Thisalsoexplainswhythe where yˆ = y/H, and N˜a is obtained by replacing σ
deviations of the solutions between different calculation 0
with η Γ˙ in Eq. (4). Here, the characteristicviscosity η
conditions for σ =0.08 are small because for σ =0.08, 0 0
0 0 is estimated using the NEMD simulation with a uniform
thelocalviscositiesdecreaseasaresultofboththe shear
temperatureT atashearrateγ˙ =1×10−6asη =1930.
thinning and the temperature increase so that the rela- 0 0
InFig. 9,itcanbeobservedthatthetemperaturepro-
tion η ≪ ∆x2/∆t is satisfied everywhere, except at the
file for σ =0.01 under the calculation condition C1 de-
vicinity of the plate, where the temperature is close to 0
viates considerably from the analytic solution; however,
the reference temperature. Incidentally, a similar condi-
underthe calculationconditionC6,wherethe numberof
tion is also required for the energy transport equation,
polymer chains in the MD cell is increased tenfold, the
but it is usually satisfied as long as that for the momen-
temperature profile becomes much closer to the analytic
tum transport equation is satisfied because the thermal
solution. Note that for σ =0.01, the effect of the shear
conductivity is usually small compared to the viscosity 0
thinning of the viscosity is very small, as observed in
(i.e., the Prandtl number is small) for polymeric liquids.
Fig. 4. Thus, the analytic solution is thought to repre-
sent an accurate temperature profile for which the noise
effects in the SMD simulation are completely neglected.
The temperature profiles for σ = 0.01 do not coin- Thisfactindicatesthatthe temperature profileobtained
0
cide with each other, although the absolute differences by the SMD simulation for small applied shear stresses
are small, i.e., at most one percent of the wall temper- is greatly affected by the noise even after the long-time
ature T , while for σ = 0.08, the deviations between average is taken. Thus, to obtain an accurate temper-
0 0
differentcalculationconditions areverysmall. This may ature profile for a small applied shear stress, one needs
be caused by the fluctuations of local shear rates and a large number of molecules in each MD cell, although
stresses because the signal-to-noise ratio is smaller for the cell size l should be smaller than the mesh inter-
MD
σ = 0.01 than for σ = 0.08. However, by increasing val ∆x for an efficient computation. The temperature
0 0
the particle numbers in each MD cell, the fluctuations profile obtained using the SMD simulation for σ =0.05
0
due to the noise arising from the MD simulations can alsodeviatesfromtheanalyticsolution;thisiscausedby
be reduced. Figure 9 shows the comparisons between shear thinning, which is not considered for the analytic
the SMD simulations and the analytic solutions for dif- solution.
ferent applied shear stresses. Here, the analytic solution Inadditiontotheoveralltechnicalissues,therearealso
is obtained for a Newtonian liquid with an exponential concernsabouttheconceptualissuesoftheSMDmethod.
dependence of the viscosity on temperature[42, 43], i.e., The synchronous scheme via the macroscopic transport
η(T) = η exp[−(T − T )/∆T ], where η is a charac- equations imposes ignoring the molecular transports of
0 0 0 0
teristic viscosity at a uniform temperature T and ∆T constituents across the mesh interval. Thus, the SMD
0 0
is defined above Eq. (4). The analytic solution for the methodisnotapplicableforthedilutepolydispersefluids
8
butisratherdesignedfordensefluidssuchasthepolymer convection can be ignored at those scales, i.e., the
melt. In the concept of locality of the SMD method, the local Reynolds numbers at those scales must be
viscous diffusion and the vortex structures are resolved very small.
in the global mesh system. Thus, the mesh interval and
the size of the MD cell must be sufficiently small such
that the fluid inertia and convection can be ignored at The first issue concerns the numerical stability con-
those scales,i.e., the localReynolds number atthe mesh dition for the macroscopic transport equations Eq. (3),
interval ∆x must be very small. This is also related to althoughtheexactstabilityconditionisunknownforthe
the technical aspects of the SMD method because this SMD method. The second issue is the consequence of
condition warrantsexploiting the SLLOD algorithmand the comparisons between the SMD simulations and the
thehomogeneousrescalingofthe kineticenergiesineach analytic solutions. For the present problem, an SMD
MD cell. simulation using 10,000 particles in each MD cell can
successfully reproduce an accurate temperature profile
that is described by the analytic solution. The third is-
VI. SUMMARY sue is concerned with the locality concept of the SMD
method,andthisalsoinvolvesthetechnicalaspectsofthe
SMDmethodbecausetheuseofthesmalllocalReynolds
We have proposeda synchronizedmolecular dynamics
number warrants exploiting the SLLOD algorithm and
simulationviamacroscopicheatandmomentumtransfer
thehomogeneousrescalingofthekineticenergiesateach
andappliedthismethodtotheanalysisofthelubrication
time interval ∆t, which are described in Sec. III.
of a polymeric liquid, coupled with viscous heating. The
Although these issues must be considered, the SMD
rheologicalpropertiesandtheconformationsofthepoly-
simulation has two distinctive advantages over the full
mer chains are investigated using a non-dimensional pa-
MD simulations: First, the SMD simulation can reduce
rameter,i.e.,theNahme-Griffithnumber. TheSMDsim-
the computational effort (i.e., the number of molecules)
ulationdemonstratesthatstrongshearthinning,whichis
compared to that of the full MD simulation by a factor
almost inversely proportional to the shear rate, and the
of (l /∆x)d, where d is the dimension number of the
transitional behavior for the conformation of the poly- MD
macroscopic transport equation. In the present simula-
merchainsoccurwitharapidtemperatureincreasewhen
tion, the factor is 0.128 because d is one. However, the
the Nahme-Griffith number exceeds unity. The results
extensionto the two-and three-dimensionalcases is also
show that the linear stress-optical relation holds despite
possible by incorporating the algorithms developed in,
thecomplicatedbehaviorsofthetemperature,shearrate,
for example, Refs. 16–19. Second, almost perfect par-
and conformation of the polymer chains.
allelization efficiency (i.e., using N CPUs in the parallel
Wehavealsocarriedoutnumericaltestsundervarious
computation speeds up the calculation compared to us-
calculationconditionsbyvaryingthenumberofmeshin-
ing a single CPU by N times) is achieved in a parallel
tervalsM (i.e.,∆x=H/M),timeinterval∆t,andnumber
computation with as many CPUs as there are MD cells
of polymer chains in eachMD cell N and found the fol-
p
byassigningeachCPUto anMD cell.[21]Thisefficiency
lowing critical issues in the implementation of the SMD
isobtainedbecauseeachMDsimulationisperformedin-
simulations:
dependently in the time interval ∆t. The advantage in
1. The time interval∆t mustbe smallerthanthe vis- parallel computation holds not only for the short-range
cousdiffusiontime∆x2/η,withinwhichtheviscous interaction molecular models but also for any compli-
force propagates for the mesh interval ∆x; other- cated molecular models for which the parallelization of
wise, the SMD simulations reproduce the solutions MD simulations is difficult. These advantages enable us
withapseudoviscositythatissmallerthanthetrue toanalyzethecomplicatedflowbehaviorsofcomplexliq-
viscosity. uids at the macroscopic scales found in actual engineer-
ingandbiologicalsystemsonthebasisoftheappropriate
2. The temperature profile at a small applied shear
molecular model.
stress is strongly affected by the noise arising from
the localMDcells. Toobtainanaccuratesolution,
taking the long-time averages as well as increasing
ACKNOWLEDGMENTS
thenumberofparticlesineachMDcellisrequired.
3. Concerning the concept of locality in the SMD This study was financially supported by the Hyogo
method, the flow behaviors involving the viscous Science and Technology Association and by the Grant
dissipation and the vortex structures are resolved for Basic Science ResearchProjects from The Sumitomo
intheglobalmeshsystem. Thus,themeshinterval Foundation. The computations havebeen performed us-
∆x and the size of the MD cell l must be ing the facilities of the Supercomputer Center at the In-
MD
sufficiently small such that the fluid inertia and stitute for Solid State Physics, the University of Tokyo.
9
[1] M.P.AllenandD.J.Tildesley, Computer Simulation of polymer melt between rapidly oscillating plates: an ap-
Liquids, (Oxford University Press, NewYork, 1989). plication of multiscale modeling”, Europhys. Lett. 86,
[2] D. J. Evans and G. Morris, Statistical mechanics of 18002 (2009).
nonequilibriumliquids,(Cambridgeuniversitypress,New [21] S. Yasuda and R. Yamamoto, “Multiscale modeling and
York,2008). simulationforpolymermeltflowsbetweenprallelplates”,
[3] R.B.Bird,R.C.Armstrong,andO.Hassager,Dynamics Phys. Rev.E 81, 036308 (2010).
of polymeric liquids Vol. 1 (John Wiley and Sons, New [22] S.YasudaandR.Yamamoto,“Dynamicrheologyofasu-
York,1987). percooled polymer melt in nonuniform oscillating flows
[4] M.LasoandH.C.O¨ttinger,“Calculation ofviscoelastic between rapidly oscillating plates”, Phys. Rev. E 84,
flow using molecular models: the CONNFFESSIT ap- 031501 (2011).
proach”, J. Non-Newtonian Fluid Mech. 47, 1 (1993). [23] T. Murashima, S. Yasuda, T. Taniguchi, and R.
[5] K.Feigl,M.Laso,andH.C.O¨ttinger,“CONNFFESSIT Yamamoto, “Multiscale Modeling for polymeric flow:
approach for solving a two-dimensional viscoelastic fluid particle-fluidbridgingscalemethods”,J.Phys.Soc.Jpn.
problem”, Macromolecules 28, 3261 (1995). 82, 012001 (2013).
[6] M. Laso, M. Picasso, H. C. O¨ttinger, “2-D time- [24] A.Jabbarzadeh,J.D.Atkinson,andR.I.Tanner,“Wall
slipinthemoleculardynamicssimulationofthinfilmsof
dependent viscoelastic flow calculations using CON-
NFFESSIT”, AIChEJ. 43, 877 (1997). hexadecane”, J. Chem. Phys. 110, 2612 (1999).
[7] M.Dressler,B.J.Edwards,O¨ttinger,“Macroscopicther- [25] R.Khare, J. dePablo, and A.Yethiraj, “Molecular sim-
ulation and continuum mechanics investigation of vis-
modynamics of flowing polymeric liquids”, Rheol. Acta
38, 117 (1999). coelasticpropertiesoffluidsconfinedtomolecularly thin
films”, J. Chem. Phys.114, 7593 (2001).
[8] W. E and B. Engquist, “The heterogeneous multi-scale
methods”, Comm. Math. Sci. 1, 87 (2003). [26] R.Khare, J. dePablo, and A.Yethiraj, “Molecular sim-
ulation and continuum mechanics study of simple fluids
[9] W. Ren and W. E, “Heterogeneous multiscale method
in non-isothermal planar couette flows”, J. Chem. Phys.
forthemodelingofcomplexfluidsandmicro-fluidics”,J.
Compt. Phys. 204, 1 (2005). 107, 2589 (1997).
[27] B. H. Kim, A. Beskok, and T. Cagin, “Viscous heat-
[10] M. M¨uller and K. C. Daoulas, “Speeding Up Intrinsi-
ing in nanoscale shear driven liquid flows”, Microfluid
cally Slow Collective Processes in Particle Simulations
Nanofluid 9, 31 (2010).
by Concurrent Coupling to a Continuum Description”,
[28] K.KremerandG.S.Grest,“Dynamicsofentabgledlin-
Phys.Rev.Lett. 107, 227801 (2011).
ear polymer melts: A molecular-dynamics simulation,”
[11] M. K. Borg, D. A. Lockerby,J. M. Reese, “A multiscale
J. Chem. Phys. 92, 5057 (1990).
method for micro/nano flows of high aspect ratio”, J.
[29] B.H.A.A. van den Brule and S.B.G. O’Brien,
Compt. Phys. bf233, 400 (2013).
“Anisotropic conduction of heat in a flowing poly-
[12] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P.
meric material”, RheolActa 29, 580 (1990).
G. Kevrekidis, O. Runborg, and C. Theodoropoulos,
[30] H.C.O¨ttingerandF.Petrillo,“Kinetictheoryandtrans-
“Equation-free, coarse-grained multiscale computation:
portphenomenaforadumbbellmodelundernonisother-
enablingmicroscopicsimulationstoperformsystem-level
analysis”, Comm. Math. Sci. 1, 715 (2003). mal conditions”, J. Rheol. 40, 857 (1996).
[31] R.B. Bird andC. F.Curtiss, “Nonisothermal polymeric
[13] I. G. Kevrekidis and G. Samaey, “Equation-free multi-
fluids”, Rheol. Acta 35, 103 (1996).
scale computation: algorithms and applications”, Annu.
Rev.Phys.Chem. 60, 321 (2009). [32] R. B. Bird, C. F. Curtiss, and K. J. Beers, “Polymer
contributiontothethermalconductivityandviscosityin
[14] S. De, J. Fish, M. S. Shephard, P. Keblinski, and S.
a dilute solution”, Rheol. Acta 36, 269 (1997).
K. Kumar, “Multiscale modeling of polymer rheology”,
Phys.Rev.E 74, 030801(R) (2006). [33] D.C. Venerus,J.D.Schieber, H.Iddir,J. Guzma´n,and
A.Broerman,“AnisotropicThermalDiffusivityMeasure-
[15] S. De, “Computational study of the propagation of the
ments in Deforming Polymers and the Stress-Thermal
longitudinalvelocityin apolymermelt containedwithin
Rule”, Int. J. Thermophysics 22, 1215 (2001).
a cylinderusing a scale-bridging method”, Phys. Rev.E
88, 052311 (2013). [34] J. D. Schieber, D. C. Venerus, K. Bush, V. Balasubra-
manian, and S. Smoukov, “Measurement of anisotropic
[16] T.MurashimaandT.Taniguchi,“MultiscaleLagrangian
energy transport in flowing polymers by using a holo-
Fluid Dynamics Simulation for Polymeric Fluid” J.
Polym. Sci. B 48, 886 (2010). graphic technique”,PNAS 101, 13142 (2004).
[35] J. D. Schieber, D. C. Venerus, and S. Gupta, “Molecu-
[17] T.MurashimaandT.Taniguchi,“Multiscale Simulation
lar origins of anisotropy in the thermal conductivity of
ofHistoryDependentFlowinPolymerMelt”,Europhys.
Lett.96, 18002 (2011). deformed polymer melts: stress versus orientation con-
tributions”. Soft Matter 8, 11781 (2012).
[18] T.MurashimaandT.Taniguchi,“Multiscale Simulation
[36] S. Gupta, J. D. Schieber, and D. C. Venerus,
of History Dependent Flow in Polymer Melt”, J. Phys.
Soc. Jpn. 81, SA013 (2012). “Anisotropicthermalconductioninpolymermeltsinuni-
axial elongation flows”, J. Rheol. 57, 427 (2013).
[19] S. Yasuda and R. Yamamoto, “A model for hybrid sim-
[37] C.J.Pipe, T.S.Majmudar, andG.H.McKinley,“High
ulation of molecular dynamics and computational fluid
dynamics”, Phys. Fluids 20, 113101 (2008). shear rate viscometry”, Rheol Acta 47, 621 (2008).
[38] R. Yamamoto and A. Onuki, “Dynamics and rheology
[20] S.Yasudaand R.Yamamoto, “Rheological propertiesof
of a supercooled polymer melt in shear flow,” J. Chem.
10
Phys.117, 2359 (2002). 2007).
[39] H. Janeschitz-Kriegl, Polymer Melt Rheology and Flow [42] J. Gavis and R. L. Laurence, “viscous heating in plane
Birefringence (Springer, Berlin, 1983). and circular flow between moving surfaces”, I&EC Fun-
[40] M. Doi and S. F. Edwards, The theory of polymer dy- damentals 7, 232 (1968).
namics(Clarendon, Oxford, 1986). [43] L.E. Becker and G.H. McKinley, “The stability of vis-
[41] G.Strobl,ThePhysicsofpolymers(Springer,Heidelberg, coelastic creeping plane shear flows with viscous heat-
ing”, J. Non-Newton.Fluid Mech. 92, 109 (2000).