Multiscale Simulations for Polymeric Flow

2.4 MB·

by  Takahiro Murashima

Multiscale Simulations for Polymeric Flow
Takahiro Murashima , Takashi Taniguchi , Ryoichi Yamamoto , and Shugo Yasuda
∗ † ‡ §
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan, and
CREST, Japan Science and Technology Agency, Kawaguchi 332-0012, Japan.
(Dated: January 7, 2011)
Multiscale simulation methods have been developed based on the local stress sampling strategy
andappliedtothreeflowproblemswithdifferentdifficultylevels: (a)generalflowproblemsofsimple
fluids, (b) parallel (one-dimensional) flow problems of polymeric liquids, and (c) general (two- or
three-dimensional) flow problems of polymeric liquids. In our multiscale methods, the local stress
of each fluid element is calculated directly by performing microscopic or mesoscopic simulations
1 according to the local flow quantities instead of using any constitutive relations. For simple fluids
1 (a), such as the Lenard-Jones liquid, a multiscale method combining MD and CFD simulations is
0 developed based on the local equilibrium assumption without memories of the flow history. The
2 results of the multiscale simulations are compared with the corresponding results of CFD with
or without thermal fluctuations. The detailed properties of fluctuations arising in the multiscale
n
simulations are also investigated. Forpolymericliquids in parallel flows (b),themultiscale method
a
J isextendedtotakeintoaccountthememoryeffectsthatariseinhydrodynamicstressduetotheslow
relaxationofpolymer-chainconformations. Thememoryofpolymerdynamicsoneachfluidelement
6
isthusresolvedbyperformingMDsimulationsinwhichcellsarefixedatthemeshnodesoftheCFD

simulations. The complicated viscoelastic flow behaviours of a polymeric liquid confined between
]
t oscillatingplatesaresimulatedusingthemultiscalemethod. Forgeneral(two-orthree-dimensional)
f
o flowproblemsofpolymericliquids(c),itisnecessarytotracethehistoryofmicroscopicinformation
s such as polymer-chain conformation, which carries the memories of past flow history, along the
t. streamlineofeachfluidelement. ALagrangian-based CFDisthusimplementedtocorrectly advect
a the polymer-chain conformation consistently with the flow. On each fluid element, coarse-grained
m
polymersimulationsarecarriedouttoconsiderthedynamicsofentangledpolymerchainsthatshow
- extremelyslowrelaxation comparedtomicroscopictimescales. Thismethodissuccessfullyapplied
d tosimulate a flow passing through a cylindrical obstacle.
n
o PACSnumbers: 31.15.xv46.15.-x
c Keywords: multiscale,simulation, modeling,polymericliquids,complex fluids,softmatters, viscoelastic flu-
[ ids,memoryeffect

1
v
I. INTRODUCTION are in local equilibrium at each instance. Therefore, one
1
1 must trace the “memory” or “history” of the deforma-
2 tions of each fluid element along its streamline. This
1 The prediction of flow behaviours of polymeric liquids coupling between microscale and macroscale dynamics
. usingcomputersimulationsisachallengingthemeinvar- hinders the simulation of polymeric liquids.
1
ious fields of science and engineering including physical
0 In this paper, we present a detailed explanation of a
science, materials science, and mechanical and chemical
1 multiscale method that bridges the hydrodynamic mo-
1 engineering. Polymeric liquids are known to exhibit pe- tionsoffluidsusingcomputationalfluiddynamics(CFD)
: culiar flow behaviours that are related to the microscale
v and the microscopic [or mesoscopic] dynamics of poly-
dynamics of their polymer chains; these dynamics affect
i mer configurations using molecular-dynamics (MD) [or
X theviscoelasticity,shearthinning/thickeningbehaviours, coarse-grained (CG)2–6] simulations. The concept of
r andflow-inducedphasetransitionsofthese liquids.1 The bridging microscale and macroscale dynamics is also im-
a
characteristic times of the microscale dynamics of poly-
portantforotherflowproblemsofsoftmatterswithcom-
mer chains tend to be very long and are often compara-
plex internal degrees of freedom (e.g., colloidal disper-
ble to the time scales of macroscalefluid motions. Thus,
sion, liquid crystal, and glass). The basic idea of our
forthesecompounds,onecannotseparatethemicroscale
multiscale method can be applied to those softmatter
and macroscale dynamics in the temporal domain, even
flows.
though a large gapexists between those length and time
In non-Newtonian fluid dynamics, many novel CFD
scales; one also cannot make the assumption usually
schemes have been proposed for the viscoelastic flows of
made for simple fluids under flow that the fluid elements
polymeric liquids.7–11 In CFD methods, a model consti-
tutive equation is used to determine the local stresses
at each instant from the history of previous velocity
fields. For dense polymeric liquids (melts), however, the
∗Electronicmail: [email protected]
detailed form of the constitutive relations is so compli-
†Electronicmail: [email protected]
‡Electronicmail: [email protected] catedthattheyaregenerallyunknown.12Thus,theusual
§Electronicmail: [email protected] CFD methods, which require constitutive equations, are
2
simple fluids complex fluids complex fluids
(general flow) (parallel flow) (general flow)
y
x
(a) (b) (c)
FIG. 1: Schematic illustrations of our multiscale methods bridging CFD and MD [or CG] simulations. The developments
were carried out in astep-by-stepmannerstarting from thebridgingmethod for simplefluids(a) andproceeding toamethod
applicable to general flow problems of complex fluids. The small cubic boxes placed in thefluid systems represent MD or CG
simulation cells within which the local stress of each fluid element is sampled. Case (a) represents a cavity flow of a simple
fluid for which the local stress σ of the fluid at a time t and a position r is given by a function of the local deformation rate
γ˙ at the same time and position, σ(r,t) = f(γ˙(r,t)). Case (b) represents a polymeric liquid subject to an oscillatory shear
flow where theflow is parallel to the x-direction and thevelocity gradient exists only in the y-direction. Here, the local stress
is a function of the history of the velocity gradient in the past t′ ≤ t at the position y, σ(y,t) = F[γ˙(y,t′); t′ ≤ t]. Case (c)
represents a polymeric liquid passing through a cylindrical obstacle. Here, the local stress is a function of the history of the
velocity gradient in thepast t′≤t along thestreamline R(t′) of thefluidelement, σ(r,t)=F[γ˙(R(t′),t′); t′≤t]. Thedetails
of our multiscale methods for cases (a), (b), and (c) are given in Secs. II, III,and IV,respectively.
usually not straightforwardly applicable to the compli- ids,(b)parallel(one-dimensional)flowproblemsofpoly-
cated flow problems of polymeric liquids. Instead, mi- mericliquids,and(c)general(two-orthree-dimensional)
croscopic simulations such as MD and CG simulations flowproblemsofpolymericliquids,asschematicallyillus-
are often used to investigate the rheological properties trated in Fig. 1.
of such materials. The microscopic simulations are usu- For the simple fluids shown in Fig. 1 (a), the local
ally performed only for a tiny piece of the material in stress σ of the fluid at a time t and a position r is given
the equilibrium or non-equilibrium state under uniform byafunctionofthe localdeformationrateγ˙ atthe same
external fields of shear velocity, temperature gradient, time and position,
and electric field. Although the microscopic simulations
are even applicable to macroscopic flow behaviours, the σ(r,t)=f(γ˙(r,t)). (1)
drawbackofthistypeofsimulationistheenormouscom-
We proposed a multiscale method that combines CFD
putationtime required. Thus, forproblems thatconcern
and MD based on the local equilibrium assumption, in
large-scale and long-time fluid motions far beyond the
which a local stress immediately attains a steady state
molecular scale, which are commonly encountered in en-
after a short transient time during which a strain rate is
gineering,fully microscopicsimulationsaredifficultfrom
subjected to the fluid element. In this method, a lattice-
a practical standpoint. To overcome the weaknesses of
mesh-basedCFDschemeisusedatthemacroscopiclevel,
the individual methods, we have developed a multiscale
and the non-equilibrium MD simulations are performed
simulation method that combines CFD and MD [or CG]
simulations.13–17 insmallMDcellsassociatedwitheachlatticenodeofthe
CFDtogeneratealocalstressaccordingtoalocalstrain
Inourmultiscalemethods,macroscopicflowbehaviour rate.
is calculated by a CFD solver; however, instead of us- ThemultiscalemethodcombiningCFDandMDisex-
ing any constitutive equations, a local stress is calcu- tendedinastraightforwardfashiontothepolymericflows
lated using the MD [or CG] simulation associated with in the one-dimensional geometry shown in Fig. 1 (b), in
a local fluid element according to the local flow variable which the macroscopic quantities, e.g., velocity, temper-
obtained by the CFD calculation. The multiscale meth- ature, and stress, are uniform in the flow direction par-
odsareappliedto threeflowproblemswithdifferentlev- allel to the plates. This situation allows us to neglect
els of difficulty: (a) general flow problems of simple flu- the advection of memory on a fluid element along the
3
streamline, so the local stress of a fluid element at t and
y is a functionalofthe history ofthe velocitygradientin
the past t ≤t at the same position y,
′
σ(y,t)=F[γ˙(y,t); t ≤t]. (2)
′ ′
Forgeneraltwo-orthree-dimensionalflowproblemsof
polymericliquids, showninFig.1 (c), onemustconsider
theadvectionofpolymer-chainconformations,whichcar-
ries memory effects on a fluid element along the stream-
line. Thelocalstressisthusafunctionalofthehistoryof
the velocity gradientin the past t ≤t along the stream-
′
line R(t) of the fluid element,
′
σ(r,t)=F[γ˙(R(t′),t′); t′ ≤t]. (3)
Tomeetthisrequirement,Lagrangianfluiddynamicsare
FIG.2: Thestressrelaxation functionG(t)fortheLJliquid
implementedfortheCFDcalculationtotracethe advec-
(dashed line) and model polymermelt (solid line).
tion of a fluid element that contains the memory of the
configuration of the polymer chains. The local stresses
oneachfluid particle are calculatedusing coarse-grained
to the one-dimensional flows of polymeric liquids con-
polymer simulations in which the dynamics of entangled
fined in parallel plates. The flow behaviours of a glassy
polymer chains are calculated.
polymer melt in the oscillating plates are calculated, the
The idea of using multiscale modelling to calculate
local rheological properties of the polymer melt in the
the local stress for the fluid solver using the micro-
rapidlyoscillatingplatesareinvestigated,andtheresults
scopic simulation instead of any constitutive relations
arecomparedwiththeanalyticalsolutionforaninfinites-
was first proposed for polymeric liquids by Laso and
O¨ttinger18–20, who presented the CONNFFESSIT ap- imally small strain to demonstrate the validity of the
multiscale method. In Sec. IV, the multiscale method
proach. The multiscale method was also proposed by
ofLagrangiandynamics andCG polymer-dynamicssim-
E and Enquist,21 who presented the heterogeneous mul-
ulations is presented for a two-dimensional flow problem
tiscale method (HMM) as a generalmethodology for the
of polymeric liquids. A transient flow passing a circular
efficient numerical computation of problems with multi-
object is demonstrated. Summaries of each section and
scale characteristics. HMM has been applied to the sim-
future perspectives are given in Sec. V.
ulationofcomplexfluids22,23 butcompletelyneglectsthe
advectionof memory. The equation-free multiscale com-
putationproposedbyKevrekidiset al. isbasedonasim-
ilar idea and has been applied to various problems.24,25 II. SIMPLE FLUID
The basic idea of our multiscale modelling is the same
as those earlier proposed. This type of bridging method We consider a simple liquid composed of particles in-
has also been developed recently by several researchers. teractingviatherepulsivepartoftheLennard-Jones(LJ)
Deetal. haveproposedthescale-bridgingmethod,which potential,
cancorrectlyreproducethememoryeffectofapolymeric
liquid, and demonstrated the non-linear viscoelastic be- 4ǫ (σ/r)12−(σ/r)6 +ǫ (r ≤21/6σ),
U (r)= ,
haviourofapolymericliquidbetweenoscillatingplates.26 LJ 0 (r >21/6σ).
(cid:26) (cid:2) (cid:3)
The methodology of the present multiscale simulations (4)
forone-dimensionalflowsofpolymericliquidsisthesame where σ and ǫ are the length and energy units of the LJ
as that used in the scale-bridging method. Kessler et al. potential,respectively. Inthis section,we measurespace
have recently developed a multiscale simulation for rar- and time in the units of σ and τ = mσ2/ǫ, respec-
0
efied gas flows based on a similar idea.27 tively. The temperature T is measured in the unit of
p
Inthefollowing,thebridgingmethodofMDandCFD ǫ/k . Here, m and k are the mass of the LJ particle
B B
simulations for simple fluids is presented in Sec. II. In and the Boltzmann constant, respectively. The temper-
this section, the results of the 2-D cavity flows are com- ature T and density ρ of the liquid are assumed to be
paredwith thoseofNewtonianfluids to demonstratethe uniform and fixed as T =1.0 and ρ=0.8.
validity of our multiscale simulation method. A special At this density and temperature, the LJ liquid does
focusisplacedonthefluctuationarisinginourmultiscale not have a long-time memory. That is, the shear stress
method. We carry out spectral analysis of the fluctua- dependsonlyontheinstantaneousstrainrate,regardless
tions and compare the results with those obtained us- of the history of the previous strain rates experienced
ing fluctuating hydrodynamics. In Sec. III, the multi- by the fluid element. In Fig. 2, we show the stress-
scale method of MD and CFD simulations is extended relaxation function G(t) of the present LJ liquid and of
4
a model polymeric liquid (the latter will be discussed in 1. CFD Scheme
the next section). The stress-relaxation function G(t) is
calculated as We use a lattice-mesh-based finite-volume method
with a staggered arrangement for vector and scalar
G(t)=hΠ (t+t )Π(t )i/k TV, (5)
xy 0 0 B quantities29 (Fig. 3). The control volume for a vector
where Π is the space integral of the microscopic stress quantityisaunitsquaresurroundedbydashedlines,and
xy
tensor in the volume V. It is seen that the stress relax- that for a scalar quantity is a unit square surroundedby
ationoftheLJliquidrapidlydecreasesandthatthetime solid lines. Equations. (6) and (7) are discretised by in-
correlation of stress almost disappears in t ≪ 1. The tegrating the quantities on each controlvolume. For nu-
relaxation time τ of the LJ liquid may be estimated mericaltimeintegrations,weusethefourth-orderRunge-
LJ
as τ = 0.067 with the definition G(τ )/G(0) = e 1. Kutta method, in which a single physical time step ∆t
LJ LJ −
The time scaleoftemporalvariationsinthe macroscopic is divided into four sub-steps. The time evolution of a
flowsmaybe muchlargerthanthestress-relaxationtime quantity φ, which is to be determined by the equation
of the present LJ liquid. Thus, for the macroscopic flow ∂φ/∂t=f(t,φ), is written as
behaviours of simple fluids, one can ignore the memory
∆t
effect in a local stress due to the history of previous flow φ =φn+ f(t ,φn), (9a)
velocities. That is, we can assume the local equilibrium ∗n+12 2 n
∆t
states at any instant in the macroscopic flows. φ =φn+ f(t ,φ ), (9b)
In this section, we present a model for multiscale sim- ∗n∗+21 2 n+21 ∗n+12
uorlaytieoffnescotsf.MDThaendbrCidFgDingforscshimempleeoflfuiMdsDwiatnhdouCtFmDemis- φ∗n+1 =φn+∆tf(tn+21,φ∗n∗+12), (9c)
∆t
bscaospedicoqnutahnetiltoiecsa.l eTquhielibmriuulmtiscaassleumsipmtuiolnatoiofntshearmeapcreor-- φn+1 =φn+ 6 f(tn,φn)+2f(tn+21,φ∗n+12)+
h
formedforthedrivencavityflowsofthesimpleLJliquid, 2f(tn+1/2,φ∗n∗+1/2)+f(tn+1,φ∗n+1) .
and the results are compared with those of the Newto-
(9id)
nian fluid. Special attention is given to the efficiency of
thepresentmultiscalemethodandtothenoisearisingin The time evolution of the fluid velocity v is computed
this method. by the above set of equations. On the other hand, the
pressurepisdeterminedsothatthefluidvelocitysatisfies
the incompressible condition (6) at each sub-step. The
A. Multiscale Model procedure at each sub-step is written as
Incompressibleflowswithuniformdensityρ andtem- p=p+ψ, (10a)
0
perature T0 are described by the following equations: v =v−τ∇ψ, (10b)
∂v e 1
α =0, (6) △ψ = ∇v, (10c)
∂x e τ
α
wherepisthepressureobtainedeattheprevioussub-step,
∂vα ∂vα 1 ∂Pαβ v is the velocity obtained by solving equation (9) at the
+v = +g , (7)
∂t β∂x ρ ∂x α presentsub-step,τ isthe timeincrementofthesub-step,
β 0 β e
and ψ is the correction term to obtain the divergence-
e
where x is the Cartesiancoordinate system, t the time,
α free velocities. The remaining three components of the
v thevelocity,ρ thedensity, P thestresstensor,and
α 0 αβ stress tensor, T , are to be computed directly by MD
αβ
g the external force per unit mass. Throughout this
α simulations. The details of the method are described in
work, the subscripts α, β, and γ represent the index in
the next subsection. Note that the calculation of T is
αβ
Cartesiancoordinates,i.e.,{α,β,γ}={x,y,z},andthe
carried out at each sub-step of Eq. (9).
summation convention is used. The stress tensor P is
αβ
written in the form
2. Local MD sampling
P =−pδ +T , (8)
αβ αβ αβ
where p is the pressure and δαβ is the Kronecker delta. We compute the local stresses by MD simulations ac-
We may assume that Tαβ is symmetric, Tαβ =Tβα, and cording to the local strain rates, rather than the local
traceless, Tαα=0, for isotropic simple fluids.28 To solve flow velocities themselves. A schematic diagram of the
theaboveequations,oneneedsaconstitutiverelationfor methodisdepictedinFig. 3. AttheCFDlevel,thelocal
the stress tensor Tαβ. In our multiscale method, instead strain rate tensor Eαβ is defined as
ofusinganyexplicitformulassuchastheNewtoniancon-
stitutive relation, Tαβ is computed directly by MD sim- E = 1 ∂vα + ∂vβ , (11)
ulations. αβ 2 ∂x ∂x
(cid:18) β α(cid:19)
5
staargrgaenrgeedment rotation of MD cell ∆t
CFD
y′ y′ y′ tMD
y z′lMDx′ z′ x′ z′ x′ ∆x
Transienttime
x excludedfrom ∆tMD
z sampling MD
(a) Space decomposition (b) Time evolution
FIG.3: Schematicdiagram ofthemultiscalesimulation methodforsimplefluids. (a)Staggered arrangementofthevelocityv
andstress tensorP on alattice-mesh grid of theCFD calculation. CFD simulations areperformed inareference coordinate
αβ
(x,y,z), while MD simulations are performed in a rotated coordinate (x′,y′,z′) so that the diagonal components of E′ all
αβ
becomezerousingtheproceduredescribedinSecIIB.TheCFDsystemisdiscretisedintocubicsubsystemswhosesidelength
is ∆x. Each subsystem is associated with an MD cell whose side length is lMD, with the Lees-Edward periodic boundary
conditionundersheardeformation. Notethatthree-dimensionalMDsimulationsareusedatthemicroscopiclevelevenforthe
problemsforwhichone-ortwo-dimensional analysisisappliedatthemacroscopic level. (b)Aschematictimeevolutionofour
multi-scale method. TheCFD simulation proceeds with atimestep of ∆t,andtheMD simulation iscarried out fora lapseof
time tMD only to sample thelocal stress Tα′β at each node point and time step of CFD.
wherethe incompressiblecondition, E =0,is tobe sat- Non-equilibriumMDsimulationsforsimpleshearflows
αα
isfied. We can now define a rotation matrix Θ by which in the rotated Cartesian coordinates are performed in
the strain-rate tensor E is transformed to many MD cells according to the local strain rate E ’s
αβ ′
defined at each lattice node of the CFD. Once a local
0 E E stress tensor P is obtained at the MD level, the local
x′y x′z α′β
E′ =ΘEΘT =Ey′x 0 Ey′z , (12) stressat eachlattice node Pαβ in the originalcoordinate
E E 0 systemis obtainedbycombiningthe pressurepobtained
z′x z′y
aprioribyCFDandatensorT obtainedbysubtracting
  α′β
wherethediagonalcomponentsallvanish. Thistransfor- the isotropic normal stress components from P as
α′β
mation enables us to perform MD simulations with the
Lees-Edwards periodic boundary condition. P =ΘT[−pI+T′]Θ=−pI+ΘTT′Θ, (16)
WenotethattheLees-Edwardsperiodicboundarycon-
where I is the unit tensor.
ditioncannotcreateaflowfieldforanarbitraryvelocity-
In the MD simulations, we solve the so-called SLLOD
gradienttensor∂v /∂x inanMDcellbutcanreproduce
α β equations of motion with the Gaussian iso-kinetic
a velocity profile for an arbitrary (symmetric) strain-
thermostat:30–32
rate tensor E using three components of the velocity-
αβ
gradient tensor, e.g., ∂vx/∂y, ∂vz/∂y, and ∂vx/∂z. For drj = pj +γ˙r e , (17a)
simple fluids, the antisymmetric part of the velocity- dt m yj x
gradient tensor, Ω = 1(∂v /∂x −∂v /x ), does not
αβ 2 α β β α
affect the stress tensor Tαβ. Thus, the present multi- dpj =f −γ˙p e −ζp , (17b)
scale method using the Lees-Edwards boundary condi- dt j yj x j
tion in each MD cell is applicable to the general (three-
where e is the unit vector in the x direction and the
dimensional) flows of simple fluids. x
index j represents the jth particle (j = 1,··· ,N). r
j
Theoff-diagonalstresstensorT iscomputedaccord-
α′β and p + mγ˙r e are the position and momentum of
ing to E and then passed to CFD after transform- j yj x
α′β the jth particle, respectively, f is the force acting on
j
ing back into the original coordinates, T . For two-
αβ the jth particle due to the LJ potential in Eq. (4), and
dimensional flows [∂/∂z=0 and v =0], Θ and E are ex-
z ′ γ˙ is the shear rate experienced by each MD cell, which
pressed as
is written as γ˙ =2E in the present multiscale scheme.
x′y
Notethat,intheSLLODequations,p /mrepresentsthe
j
cosθ sinθ
Θ= , (13) deviation of velocity of each particle from the mean flow
−sinθ cosθ
(cid:18) (cid:19) velocity γ˙ryjex in the MD cell. The friction coefficient
Ex′y =Ey′x =−Exxsin2θ+Exycos2θ, (14) ζ is determined to satisfy the constant temperature con-
dition dT/dt = 0 with T = |p |2/3mN. The friction
j j
where coefficient ζ is calculated as
P
θ = 1tan−1 −Exx . (15) ζ = (fj ·pj −γ˙pxjpyj)/ |pj|2. (18)
2 (cid:18) Exy(cid:19) Xj Xj
6
We integrate Eq. (17) with Eq. (18) using the
lMD tMD N L U Re ∆x/lMD ∆t/tMD
leapfrog algorithm33 with the Lees-Edwards sheared pe-
C1 6.84 4.68 256 219 0.46 59 1.0 1.0
riodicboundaryconditioninthey directionandperiodic
C2 6.84 9.36 256 438 0.91 235 2.0 2.0
boundary conditions in each cubic MD cell.
C3 6.84 3.11 256 876 1.83 941 4.0 4.0
The space integral of the instantaneous microscopic
stress tensor reads as C4 8.62 9.36 512 438 0.91 235 1.59 2.0
C5 8.62 18.7 512 438 0.91 235 1.59 1.0
N
1 dULJ(ξ)ξαξβ C6 10.86 15.7 1024 2780 0.58 941 8.0 8.0
VP = p p − (19)
α′β m αj βj dξ ξ
j=1 allpairs
X X TABLEI: Simulationparametersforcavityflows. lMD,tMD,
and N are the side length of MD cell, sampling duration of
where we rewrite the momentum of the jth particle,
MD simulation, and number of particles contained in each
p +mγ˙r e , as p . Here, V is the volume of the MD
j yj x j MD cell, respectively. L, U, and Re are the system size, ve-
cell, i.e., V = lM3D, and ξ is the relative vector rj −rj′ locityoftheupperwall,andReynoldsnumberonthesystem,
betweenthe two particles rj and rj′. The density is also respectively. ∆x/lMD and ∆t/tMD are the efficiency factors
fixed in each MD simulation by the number of particles of the present multiscale simulations.
andtheboxsize. Thus,wecalculatethelocalstressusing
so-called NVT-ensembles, P (ρ,T,E ). The macroscopic
′ ′
stress is averagedin steady states after the transient be-
haviour has vanished. In the present computations, the of the stress-relaxation function G(t) shown in Fig. 2,
transienttimeineachMDprocessissetasatenthofthe η = ∞G(t)dt. The computational domain is divided
0
CFD time step, ∆t/10. The initial states in each MD into 32×32 uniform lattices.
R
simulationarecreatedbyre-scalingthethermalvelocities The results of the multiscale simulations are shown in
of molecules to be fixed to the local temperatures with- Figs. 4 and 5. Figure 4 shows the steady-state velocity
outchangingthemolecularconfigurationsfromthoseob- profiles time-averaged over t/∆t=[900,1000]. It is clear
tained by the previous process. that our multiscale method can successfully reproduce
In the following, ∆t and ∆x represent the time-step the characteristic flow properties of cavity flows with
andthemeshsizeofCFDcalculations,andt andl different Reynolds numbers; the size of the vortex be-
MD MD
represent the sampling duration and the side length of comes larger as the Reynolds number increases. Figure
a MD cell, respectively. The two parameters ∆t/t 5 shows the snapshots of velocity profiles for the case
MD
and ∆x/l represent the efficiency of our multiscale of Re=980 at different time steps. Here, the results ob-
MD
simulations. In the present simulations, we fix the time tained by multiscale simulations with the efficiency fac-
step of the MD simulation ∆t as ∆t =0.005. The tors ∆x/l =∆t/t =8 are compared with the results
MD MD MD MD
sampling durations of the MD simulation t are set to obtained for the Newtonian fluid. The results show that
MD
be largerthanthe correlationtime ofthe temporalshear a small vortex first appears in the upper-right corner
stress for the bulk fluids, t ≫ τ . The time step of and then moves gradually toward the centre of the box,
MD LJ
CFD ∆t is set to be small enough for the CFD solutions with the size of the vortex increasingas time passes. Al-
to be stable. thoughfluctuationsareseenintheinstantaneousvelocity
profiles, our multiscale method can successfully repro-
duce the characteristic behaviours in the time-evolution
B. Driven cavity flows of driven cavity flow.
The comparisons of the spatial variations of local ve-
Thesimple LJliquidiscontainedinasquareboxwith locities, strain rates, and stresses obtained by the mul-
a side length L. At t=0, the upper wall starts to move tiscale simulation and those of the Newtonian fluid are
from left to right at a velocity V . A non-slip bound- shown in Figs. 6 and 7. It is seen that the deviations of
w
ary condition is applied at each wall: v =V and v =0 the local shear stresses obtained by the multiscale sim-
x w y
at y=L and v =v =0 at the other walls. At the up- ulations are larger than those of the local velocities and
x y
per left and right corners, v =V and v =0 are applied. strain rates. The instantaneous fluctuations of the local
x w y
Although the non-slip condition causes singularities in shearstressesarenotable,buttheyaresmoothedbytak-
the strain rates and stresses at the upper corners in the ing the time averages. The fractional RMS deviations of
continuum description, we use the conventional non-slip the local velocities and stresses obtained by the multi-
boundary condition even at the corners to compare the scale simulations from those obtained by the usual CFD
resultstothoseobtainedbytheusualCFDcalculation.34 are also shown in Table II. The comparisons between
We notethatmultiscalesimulationstoresolvethe singu- C2, C4, and C5 show the effect of the efficiency factors
lar forces at the corners can be found in Refs. 35,36. ∆t/t and ∆x/l on the RMS deviations with fixed
MD MD
The multiscalesimulationsareperformedfor the param- system parameters L, U, and Re. The RMS deviations
eters listed in Table I. Here, the Reynolds number is increase as the efficiency factors increase with fixed sys-
definedasρUL/η,andtheviscosityofthecorresponding tem parameters, e.g., L, U, and Re. On the other hand,
LJ liquid is η=1.7, which is calculated by the integral in the comparisons between C1 – C3, the RMS devia-
7
FIG. 4: The steady-state velocity profiles of the cavity flows for (a) Re=59 (C1 in Table I), (b) Re=235 (C2 in Table I), and
(c) Re=941 (C3 in Table I).The velocity profiles are time averaged overt/∆t=[900,1000].
RMS deviation thanthat of a full MD simulation. The largerthe ratios,
themoreefficientthesimulations;however,thestatistical
v Pxy
fluctuations also increase. In the following part, we will
C1 0.049 0.36
discuss the nature of the fluctuations inmore detail.37,38
C2 0.034 0.25
To handle the statistical noise explicitly, we rewrite Eq.
C3 0.024 0.39 (16) as
C4 0.024 0.17
C5 0.017 0.12 P =−pI+ΘT(T′+R′)Θ, (20)
∗
C6 0.022 0.36
wheretheoff-diagonalstresstensorT ,whichistobede-
′
termined by MD sampling, is decomposed into the non-
fluctuatingstressT andthefluctuatingrandomstressR
TABLE II: Fractional root-mean-square (RMS) de- ′ ′
due to the thermal∗noise. The magnitude of each com-
viations in the steady states, which are defined as
ponent of the random stress included in MD sampling,
q(QR=0Lvd)xa2nRdTfdotr(tQhe−sQheNaSr)2s/trLe2ssT(/QM=aPx|xQy)N.SH|,efroer, Tthreepvreelsoecnittys hcRooM′2rDdipnqaitweshaenredpdoanndotqforlelporwestehnetstuhmeminadteioxnincoCnvaertnetsioiann,
the sampling duration in the calculation of the RMS in the
should depend both on the size of the MD cell l and
steady state. MD
thelengthoftimet overwhichtheaverageistakenat
MD
the MD level; hR2 i=hR¯ (l ,t )2i, where R¯(l,t)
M′ Dpq pq MD MD
represents the random stress tensor averaged in a cubic
tions of the velocities slightly decrease in the order C1
with a side length l and over a time duration t.
– C3, although the efficiency factors increase in the or-
AttheCFDlevel,whichisdiscretisedwithameshsize
der C1–C3. This unexpected finding is explained by the
∆xandatime-step∆t,thephysicallycorrectmagnitude
fact that the local strain rates increase as the Reynolds
shouldbehR2 i=hR¯ (∆x,∆t)2i. Ifthecentral-limit
number becomes larger,andthus the localshearstresses C′ FDpq pq
theorem, hR¯ (l,t)2i∝ 1/l3t, is assumed, the following
are also large in comparison to the noise intensity. The pq
simple formula can be used:
comparison of C3 and C6 also shows that the fractional
RMS deviations decrease even though the efficiency fac- ∆x 3 ∆t
tors increase. It should be noted that, in the simulation hR2 i= hR2 i. (21)
M′ Dpq l t C′ FDpq
parameter C6, the number of particles contained in each (cid:18) MD(cid:19) (cid:18) MD(cid:19)
MDcellisfourtimesgreaterthaninC3andthesampling
This line of reasoning leads to the following expression
durationofthe MD simulationis aboutfive times longer
for the correctly fluctuating stress tensor P:
than in C3. The values of C3 and C6 in Table II can
be explained by the property of fluctuations discussed
3
l t
bcaenlowb.e Tpehrefsoermfaecdtswiintdhichaitgehtehffiatcimenuclytisfcaaclteorssimfourlatlaiornges P =−pI+ΘT T∗′+s ∆MxD ∆MDt RM′ DΘ (22)
(cid:18) (cid:19) (cid:18) (cid:19)
Reynolds numbers and large system sizes.
 
As mentioned above, the ratios ∆t/t and ∆x/l to be used in CFD instead of Eq. (16). This equation
MD MD
measure the efficiency of our multiscale simulations. For indicates that if we couldre-weightthe randomly fluctu-
example,inthecaseof∆t/t =∆x/l =8,thecompu- ating part R while the non-fluctuating part T remains
MD MD ′ ′
tational efficiency is approximately 82 ×8 times better untouched, hydrodynamic simulations includin∗g correct
8
FIG. 5: Time evolutions of the velocity profile for the cavity flow with Re=941. The left column shows the results of the
multiscalesimulation with∆x/lMD =∆t/tMD =8(C6inTableI)andtherightcolumnshowsthecorrespondingCFDresults.
thermalfluctuationscanbeperformedwithinthepresent transformation of T is defined as
x′y
framework. Incidentally,wecanexplaintheresultsofthe
fractional RMS deviations in Table II using the relation 2M 12M 1
(21) while noting that the quantity R is normalised by Π {k}= 1 − − Tˆ {x}exp(−ik·x), (23)
ρU2. x′y 4M2 x′y
nXx=0 nXy=0
where x = (n ∆x,n ∆x) is the position of each lat-
x y
tice node (n ,n ), k = (2πm /L,2πm /L) is the wave
x y x y
vector, n ,n ,m ,m are integers, M is the lattice
x y z y
We note that the important key toward the develop- number in each x- and y-axis, and Tˆ {x} is defined
x′y
ment of the multiscale simulation is the separationof T
′ as Tˆ {x}=T (x + ∆x/2,y + ∆x/2) for 0 ≤ x,y ≤
andR. We thus carriedoutthe spectralanalysisforthe∗ x′y x′y
fluctua′tions in the total stress tensor computed directly L, Tˆx′y{x}=Tx′y{2L − x,y} for L < x ≤ 2L, and
fromMDsimulations,T =T +R. ThediscreteFourier Tˆ {x}=T {x,2L−y} for L<y ≤2L.
′ ′ ′ x′y x′y
∗
9
FIG. 6: Comparisons between theinstantaneous profilesof velocity Uvˆα, strain rate(U/L)Eˆxy, and shearstress (ρU2)Pˆxy on
theline of x/L=0.5 obtained by themultiscale simulation with ∆x/lMD =∆t/tMD =8 (C6 in Table I) and those obtained by
CFD for the cavity flow for Re=941. The symbols show the results of the multiscale simulation and the solid lines those of
CFD. The dashed lines in (b) and (c) show thetime-averages over t/∆t=[900,1000].
FIG. 7: Comparisons of the instantaneous profiles of velocity Uvˆα, strain rate (U/L)Eˆxy, and shear stress (ρU2)Pˆxy on the
lineofy/L=0.5obtainedbythemultiscalesimulationwith∆x/lMD =∆t/tMD =8(C6inTableI)andthoseobtainedbyCFD
for the cavity flow for Re=941. See also thecaption in Fig. 6.
The power spectra h|Π {k}|2i calculated from our using the fluctuation-dissipation theorem in the case of
x′y
multiscale simulations are plotted in Fig. 8 (a) for the ∆x/l =∆t/t =1.
MD MD
case of ∆t/t = ∆x/l = 1, which corresponds to In Fig. 9, one sees how the fluctuations depend on
MD MD
the case of Fig. 4 (a) (C1 in Table I). The angle bracket the ratios ∆x/l and ∆t/t . Here, in comparison to
MD MD
h···i indicates the time average taken at steady state at the reference case (a) [∆x/l =∆t/t =2], the num-
MD MD
the CFD level. One can see that the overall structure is ber of particles used in the MD simulations are dou-
rathersimple. Thereexistarelativelylargepeakaround bled in the case of (b) [∆x/l =1.59, ∆t/t =2], and
MD MD
k = 0 and rather flat distributions throughout the k both the number of particles and the sampling duration
plane. Theformercorrespondstothe contributionsfrom used to take the time average are doubled in the case
the non-fluctuating part T , and the latter corresponds of (c) [∆x/l =1.59, ∆t/t =1]. The noise intensity
′ MD MD
to the contributions from∗the random stress R. The decreases with decreasing ratios ∆x/l and ∆t/t ,
′ MD MD
same quantity obtained by conventional fluctuating hy- consistent with the central-limit theorem Eq. (21); i.e.,
drodynamics using a constant Newtonian viscosity and the noise intensity in (b) is approximately half that in
a random stress, the intensity of which is determined by (a), and the intensity in (c) is about one quarter of that
the fluctuation-dissipation theorem37, is shown in Fig. 8 in (a).
(b) for comparison.39 Those two plots are surprisingly From the overall properties, we can confirm that the
similar , including the fluctuation part. This similarity fluctuations arising in the present multiscale simulations
means that our multiscale simulation generates fluctu- are white noise, which obeys the central-limit theorem,
ations quite consistent with fluctuating hydrodynamics writtenasEq. (21)or(22). Theresultsalsosuggestthat
10
(cid:520)(cid:91)(cid:92) (cid:95)(cid:21) (cid:520)(cid:91)(cid:92) (cid:95)(cid:21)
(cid:552)(cid:56)(cid:21) (cid:552)(cid:56)(cid:21)
noise part
10(cid:80)(cid:91)20 30 10(cid:80)(cid:91)20 30
2E-05 2E-05
30 30
20 20
0 10 (cid:80) 20 30 10(cid:80)(cid:92) 0 10 (cid:80) 20 30 10(cid:80)(cid:92)
(cid:91) (cid:91)
(cid:11)(cid:68)(cid:12) (cid:48)(cid:88)(cid:79)(cid:87)(cid:76)(cid:86)(cid:70)(cid:68)(cid:79)(cid:72) (cid:86)(cid:76)(cid:80)(cid:88)(cid:79)(cid:68)(cid:87)(cid:76)(cid:82)(cid:81) (cid:11)(cid:69)(cid:12) (cid:41)(cid:79)(cid:88)(cid:70)(cid:87)(cid:88)(cid:68)(cid:87)(cid:76)(cid:81)(cid:74) (cid:43)(cid:92)(cid:71)(cid:85)(cid:82)(cid:71)(cid:92)(cid:81)(cid:68)(cid:80)(cid:76)(cid:70)(cid:86)
FIG. 8: The fluctuations of T′ for the case of cavity flow with Re=59. The power spectrum h|Π′ {k}|2i for the present
xy xy
multi-scale model with ∆x/lMD = ∆t/tMD =1 is shown in (a); the corresponding result from the fluctuating hydrodynamics
is shown in (b) for comparison. Πx′y represents the discrete Fourier transform of Tx′y. mα is defined as mα=(L/2π)kα, where
k is the wavevector. The insets on each figure show the h|Πx′y|2i-mx plane.
FIG. 9: The fluctuations of T′ for the case of cavity flow with Re=235. The power spectra h|Π′ {k}|2i is plotted in (a) for
xy xy
the case of Fig. 4 (b). In (b), only the numberof particles is doubled; other parameters are unchanged from (a). In (c), both
the number of particles and the sampling time of Tx′y are doubled. Πx′y represents the discrete Fourier transform of Tx′y. mα
is defined as mα=(L/2π)kα, where k is thewave vector. The insets on each figureshow theh|Πx′y|2i-mx plane.
somewhite-noise-reductionalgorithm,suchasalow-pass i.e., the “memory effect”.
filter,couldbeusefulintheCFD-MDcouplingprocesses. In this section, we consider the flow behaviours of a
polymeric liquid confined in parallel plates; we assume
that the flow direction is restricted to being parallel to
the upper and lower plates and that the local macro-
III. ONE-DIMENSIONAL POLYMERIC FLOWS
scopic quantities, e.g., the flow velocities, strain rates,
and stresses, are uniform along the flow direction. This
As we have seen, the relaxation time of the local assumptionallowsus to calculate the macroscopicquan-
stressesinsimple fluids is veryshort;forthe LJ liquidin tities using the usual fixed-mesh system without tracing
the previous section, it is estimated as τ < 0.1 in the the advection of a fluid element that contains a mem-
LJ
LJ unit time (Fig. 2). Thus, we can predict the local ory of the configuration of polymer chains. We note
equilibrium state at each instant in the sampling of the that one must treat the convection of the memory along
localstressesforthe simple fluid (Fig. 3). For polymeric the streamlines inthe generaltwo- orthree-dimensional
liquids,however,therelaxationtimeismuchlongerthan flows. The extension to polymeric flows with the advec-
that of simple fluids, and thus it happens to be larger tion of the memory will be given in the next section.
than the characteristictime of the macroscopicflows. In
Fig. 2, we show the stress relaxation of a model poly-
mer melt, whichis discussedin this section. Inthis case, A. Multiscale modelling and the simulation
wecannotpredictthe localequilibriumstates withinthe method
time-step duration to resolve the macroscopic flow be-
haviours. Instead, we must consider the history of the We consider a polymer melt with uniform density ρ
0
slow dynamics of polymer chains in each fluid element, and temperature T between two parallel plates (Fig.
0