driven Mixing by Rayleigh-Taylor Instability in the Mesoscale
Modelled with Dissipative Particle Dynamics
Witold Dzwinel
AGH Institute of Computer Science, al. Mickiewicza 30, 30-059, Krakw, Poland
David A.Yuen1
Minnesota Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 554151227, USA
Abstract
In the mesoscale mixing dynamics involving immiscible fluids is truly an outstanding problem in
many fields, ranging from biology to geology, because of the multiscale nature, which causes
severe difficulties for conventional methods using partial differential equations. The existing
macroscopic models incorporating the two microstructural mechanisms of breakup and
coalescence do not have the necessary physical ingredients for feedback dynamics. We
demonstrate here that the approach of dissipative particle dynamics (DPD) does include the
feedback mechanism and thus can yield much deeper insight into the nature of immiscible
mixing. We have employed the DPD method for simulating numerically the highly nonlinear
aspects of the Rayleigh-Taylor (R-T) instability developed over the mesoscale for viscous,
immiscible, elastically compressible fluids. In the initial stages we encounter the spontaneous,
vertical oscillations in the incipient period of mixing. The long term dynamics are controlled by the
initial breakup and the subsequent coalescence of the microstructures and the termination of the
chaotic stage in the development of the R-T instability. In the regime with high capillary number
breakup plays a dominant role in the mixing, whereas in the low capillary number regime the flow
decelerates and coalescence takes over and causes a more rapid turnover. The speed of mixing
and the turnover depend on the immiscibility factor, which results from microscopic interactions
between the binary fluid components. Both the speed of mixing and the overturn dynamics
depend not only on the mascrocopic fluid properties, but also on the breakup and coalescent
patterns, and most importantly on the nonlinear interactions between the microstructural
dynamics and the large-scale flow .
Keywords:
Rayleigh-Taylor
mixing,
mesoscopic
fluids,
immiscible
fluids,
breakup
and
coalescence morphology, dissipative particle dynamics (DPD)
1
corresponding Author, Minnesota Supercomputing Institute, University of Minnesota, 1200 Washington Av. South,
Minnesota, 55415-1227, USA e-mail: [e-mail] fax: [phone]
1
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
INTRODUCTION
The problem of mixing and dispersion of heterogeneities in viscous and immiscibel fluids is found
universally in many fields in science and engineering [1-4]. In spite of the recent advances in
understanding the mixing of homogeneous fluids in the macroscale [1,3], realistic mixing
problems of heterogeneous multicomponent fluids are inherently difficult due to the complex,
multiple-scale nature of the flow fields, the intrinsic rheological complexity, and the cross-scaling
nature of the nonlinear physics, such as grain-boundary processes, chemical reactions and
geological processes with complex rheologies [2,4]. For these reasons, mixing problems have
traditionally been treated on a case-by-case basis. Classical modelling with partial-differential
equations [5] becomes intractable, if one wishes to examine all of the details simultaneously on a
cross-scale basis. As shown in [1-3], current models incorporate two competing processes in
mixing; breakup and coalescence. These two mechanisms involve one-sided interactions without
any feedback from the microstructural dynamics back to the global flow structure. This problem of
feedback of microstructural dynamics is very important in many applications, such as dispersion
of solid nano-particles through a polymer blend and dynamical mixing of crystals in magmatic
flows. In order to address some of the issues raised above, we will focus our attention on the
canonical problem of viscous mixing driven by the Rayleigh-Taylor instability in immiscible fluids.
Numerous factors, which influence the development of the Rayleigh-Taylor instability in a
simple macroscopic fluid, such as the surface tension, viscosity, compressibility, effects of
converging geometry, three-dimensional effects, the time dependence of the driving acceleration,
shoc ks, and variety of forms of heterogeneities, are being investigated for many years
theoretically, experimentally and by computer simulations employing computational fluid
dynamics techniques (CFD) [5-20]. The studies of coupling the breakup and coalescence
processes in the context of the R mixing in immiscible fluids in the mesoscale by employing
-T
continuum hydrodynamics equations are still difficult and not reliable.
In mixing problems dealing with immiscible fluids one needs to monitor closely the
interface separating the different components. One fundamental problem in simulating the motion
of sharp interfaces is a proper description of transition phenomena, such as merging and
reconnection of interfaces in wave dynamics. These tough numerical problems can be solved in
continuum fluid models by using the so-called front-tracking method [21-23]. In mesoscopic
systems in which thermal fluctuations cannot be neglected, the internal boundary structure
remains unknown, which means that the dynamical equations of the interface cannot be solved.
This point has stressed in numerical simulations of jets in the nanoscale [24].
The processes of breakup and coalescence in mixing can in fact be viewed as population
dynamics of microstructural structures [1,3], which is driven and maintained by the chaotic flow.
Such an idea has inspired us to employ the discrete-particle approach for coupling the breakup
and coalescencing processes present in the mixing of immiscible binary fluids.
2
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
The microscopic molecular dynamics (
MD) approach can be successfully employed in
6
2
large-scale simulations of realistic microscale flows over an area of 10 . They have revealed
the formation of bubbles and spikes for the R mixing [25-28] and breakup of nanojets created
-T
by injection of fluid through nanoscale convergent gold nozzles [24]. However, the MD
simultations take up too much computer time in extending the same technique for simulating
complex fluids over the meoscopic length scale of 100 to 1000 (angstroms).
The extension of discrete-particle method to larger spatial scales can be realized by
changing only the notion of the inter-particle interaction potential by treating a large-sized particle
as a cluster of computational molecules. This idea of up-scaling has already been adapted in the
dissipative particle dynamics (DPD) method [29-39]. This method is based similarly on the
principles of discrete particle mechanics connected by two-body potentials, like MD. The
dissipative particles can be viewed as "droplets" of MD atoms. The DPD method is intrinsically
mesoscopic in nature, because it can only resolve the center-of-mass motion of the droplet and
avoids the description of its internal state. It can potentially link the microscale at tens of angstrom
scale with the macroscopic flow at the millimeter level. The key factor, which dictates the use of
the DPD method in simulating the R-T mixing of immiscible fluids is its recent successful
application in simulating complex fluids and phase separation [33-35]. By using a two-level MDDPD technique, we have successfully simulated the fragmentation, agglomeration and ordering in
colloidal arrays [35-39].
In this paper we will employ DPD simulations to direct our attention on the characteristic
microstructures developed in the course of mixing of immiscible fluids. This study would entail
monitoring simultaneously both the break -up and coalescencing processes and also the mutual
interaction between the microstructures and the global flow.
First, we introduce the DPD model, the numerical method and the simulation conditions.
Then we investigate the mixing driven by the R-T instabilities in a rectangular box. Next, we
present and discuss the simulation results. We focus then on the temporal evolution of the
interface between two fluids and its dependence of the elastic compressibility, immiscibility and
kinematic viscosities of the two superimposed fluids. Finally, we will discuss both the advantages
and disadvantages of the DPD model in coupling in one single model the processes of breakup
and coalescence with the macroscopic fluids, at least, up to the millimeter scale.
THE DPD MODEL AND SIMULATION CONDITIONS
Numerical model
In the DPD model [29,33-35] the discrete particles move about within the confines of a
rectangular box with a height h and a length Lxy . Periodic boundary conditions are imposed along
the x-direction and reflecting boundary conditions are employed in the y-direction. We have
3
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
divided the box into two parts, with the upper part of the box filled up with heavy fluid particle (H)
and the lower part filled with lighter fluid particles (L). An external gravity field G pointing
downwards is present. The particles are defined by the mass Mi, position ri, and momentum pi.
We use the two-body, short-ranged DPD force FT= FC +FB+ FD as in [39]. This type of force
consists of conservative FC , dissipative FD and Brownian FB components. The value of FT=
FC =FB= FD =0 for rij>rc . Otherwise, we give the following definitions:
FC = 1 (rij ) e ij ,
FD = M 2 (rij ) (e ij o v ij ) e ij , F B =
ij
t
1 (rij ) e ij
(1)
where:
rc
1() and 2() - are the weight functions defined such that
n D m (r )dr = 1 for m=1,2.,
0
rij - the distance between particles i and j,
rc - a cut-off radius, for which 1(r) = 2(r)= 0,
nD - an average particle density in D-dimensional system (D dimension of the system),
e ij - a unit vector pointing from particle i to particle j,
- the scaling factor for the conservative part of collision operator,
- the scaling factor for the dissipative force,
- the scaling factor for the Brownian motion,
ij - a random variable with a zero mean and actually normalized variance.
We assume that the normalized weight functions 1(rij) and 2(rij) are linear (see Eqs.(5)) as it is
in [30-35]). According to the fluctuation-dissipation theorem they are chosen such that
2
2(rij)= 1(rij) [30-33].
The temporal evolution of the particle ensemble obeys the Newtonian equations of
motion. For integrating the equations of motion we employ the leap-frog algorithm in timesteping
n
n
for the particle positions ri , and the Adams-Bashforth scheme for the particle velocities vi and
n
momenta pi . For the two-component fluid, where k=g(i) and l=g(j) denote the types of particle i
and j (k,l {H,L}), the equations of motion in 2-D space can be represented in the following
discretized form.
n+
pi
1
2
n
= pi
1
2
kl ij
n
n
+ kl 1 rijn kl M kl 2 rijn eij o ~ ij +
vn
1 rijn eij t
t
j i
()
( )(
n+
rin+1
1
pi 2
= rin +
t
Mi
)
n+
pn =
i
pi
1
2
n
+ pi
2
()
(2)
1
2
(3)
4
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
~n =
vi
1
3
1
n
n
3 p i 2 pi 2 , ~ij = ~ n ~ n , rijn =
v n vi v j
2 Mi
()
1 r
n
ij
(
rin r n
j
)
2
, en =
ij
rin r n
j
rijn
(4)
if g (i ) = g ( j )
Mi
n
n2
3 rij
6 rij
2M M
n
(5)
= 2 1 , 2 rij = 2 1 , M ij =
i
j
if g (i ) g ( j )
rc
rc
rc n
rc n
Mi + M j
()
where:
t the timestep,
Mkl - the mass of DPD particle; for interactions between particles of different kind (k l) the center
of mass motion mean is computed.
In calculating the total force acting on each particle at a given timestep we have
employed an improved algorithm based on sorting out the linked-cells [40], which has O(N)
complexity (where N is the number of particles). The method is very efficient for serial
computations giving linear increase of CPU time with the number of particles and is more than
two times faster than the standard linked-cells method [41]. A typical simulation time for
5
5
N=1.8 10 particles in two dimensions in 1.210 timesteps for a particle density of 20 particles
within a cut-off radius rc is about 5 hours on a single R12000/300 processor of the SGI/Origin
2000 system.
Input data
In Table 1 we present the properties of the input data, which define physically the DPD particle
fluid model. The density k of the k th particle system is equal to:
k = S mk nD
(7)
where, mk is the mass of atom (or molecule) of fluid k , S - scaling factor of the spatio-temporal
scale simulated, Mk =S mk is the mass of DPD particle and nD is a molecular density. For 2
-D
2
particle system (D=2) n2 = where
= 3 n3 , n3 - particle density in 3 and - the average
-D,
distance between the neighboring particles. We assume that n2 is the same for both heavy (H)
and light (L) fluids and the density contrast H / L=5 is constant. It determines the Atwood number,
an important parameter governing the density jump in the R-T problem:
=
H L
H + L
(8)
5
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
to be equal to 2/3. In the R-T instabilities found in nature the density contrast ranges from O(1) for
geological systems to O(10) or greater in astrophysics.
Similar to the MD simulations, in [29, 33-34] dimensionless units are employed to define
the DPD particle system. However, in the former case, the temporal scale (and the timestep) is
determined by the temperature of the system related to the Lennard-Jones potential well depth,
which together with the size of atom couples spatial and temporal scales in MD simulation. For
this reason, the dimensionless units define precisely spatio-temporal scale for MD particle
system. In DPD model the interparticle force is employed instead of potential function (see Eq.(1))
[29-39]. Thus the reference level of potential energy, similar to the depth of potential well for the
L-J interactions, is lacking. The size of the DPD particle is also not precisely defined and should
be matched to the spatio-temporal scale S under interest. Moreover, the Brownian part is
dependent on the timestep t (see Eq.(1)). Therefore, the values of coefficients of DPD collision
operator, thus the physical properties of the system, will strongly depend on the scaling factor S
for a given range of spatial-temporal scales.
To control the actual scale of simulation, we have matched (see [35]) the DPD particle
system to the saturated liquid Ar thermodynamic properties by using kinetic theory equations [2932]. We employ this spatio-temporal scale as a reference point to our earlier 2 MD simulations
-D
of the R-T mixing in the microscale [26-28]. From the physical properties of H (heavy) and L
(light) fluids such as: the temperature T of the system, kinematic viscosities H, L and sound
velocities c H and cL, we can calculate the DPD interparticle force coefficients kl, kl, kl (see
Eq.(1,2)) by employing the following equations from kinetic theory [29,31] and in the continuum
approximations:
2
kl = 2 kl k BT M kl
Pkl =
n D kl r
2 D
vkT =
kl =
,
(9)
1
P = ck2
k
2
k BT
Mk
(10)
(11)
2
1 kl n r
2 D( D + 2)
(12)
rc
where the average moment values in Eqs.(10,12)
r
m
= n r m (r )dr .
0
6
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
We define also the following values, which describe the dynamical r gimes of the DPD model
e
[42]:
= rc / Dv T , rc = s , fr =
L xy
(13)
rc
where Lxy is the length of computational box in x and y directions. We assume that =100, to
obtain realistic kinematic viscosities for the particle fluids system. However, the real viscosities of
the L and H fluids can be then about 70% higher than these presented in Table 1 [31]. As shown
in [42], for <15, f r>30 and s >2.5 the transport coefficients values are predicted by the kinetic
theory with relatively high accuracy. Otherwise, we can detect quantitatively the differences
between the theoretical predictions and the numerical simulations.
Table1. Input data
a) Macroscopic data
H
H-L interface
L
Atomic mass m
40
8
DPD part. mass M
1000
200
Density [kg/m ]
1400
280
Temperature T [K]
84
84
84
2.06 10 (for =100)
3.59 10 (for =100)
4.64 10 (for =100)
a) high_c [N/m]
0.233
0.233-0.500
0.233
b) low_c [N/m]
0.0233
0.0233-0.055
0.0233
a) high_c [m/N]
31.5
38.0
b) low_c [m/N]
57.3
66.0
3
2
Viscosity [ m /s]
-7
-7
-7
P i n 2-D (Eq.(10)):
Compressibility
Gravitational accel.
11
11
10
10
2
[ m/s ]
b) Microscopic data
-9
[ 10 m ]
1.06
1.06
1.06
s =rc [in ]
2.5
2.5
2.5
[dimensionless]
12,5 50 100
12,5 50 100
12,5 50 100
h [nm]
224
224
Lxy/ h
2:1 or 8:1
2:1 or 8:1
Lxy/ h [in rc units]
84:42 or 336:42
84:42 or 336:42
Number of particles
22,500-90,000
22,500-90,000
No. part. in rc sphere
20
20
7
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
No. of timesteps
-13
t [ 10
s]
120,000
120,000
1 or 2
1 or 2
The kinetic theory equations have been developed in the limit without any conservative
forces. Therefore, the transport coefficients computed from the theory can be used as a first
approximation. For example, the kinematic viscosity is underestimated theoretically but its
approximate dependence on the density is similar in both cases [31]. For a better matching, one
can apply an iterative procedure as in [43]. The contrast in dynamic viscosities p= L/H is chosen
to be close to 1 (0.6-2.2). As suggested in [1], for such the viscosity ratio wide spectrum of
breakup patterns are observed in laboratory experiments.
The simulations were carried out for DPD particle fluids with low and high sound
velocities c (low_c and high_c systems, respectively). Bec ause the pressure Pk in fluids H and L
should be the same [44], the difference in c comes from the density contrast between H and L
fluids.
Let us define the value of
= PHL-P k,
(14)
which is equal to the difference between the partial pressures computed by using Eq.(10) for k l
(P HL) and for k=l (P k =PH =P L). The value of - called here immiscibility factor governs the
immiscibility of both the L and H fluids. For 0 the fluids are miscible. Otherwise, they separate
one from another [33-35,39]. We have conducted R mixing simulations for several values of
-T
for both low_c and high_c DPD particle systems (see table1).
As shown in [35], the DPD particle system defined in Table 1 (S 25) reveals liquid
ordering in spatial scales larger than the simple L fluid (see Fig.1). While for greater S we can
-J
observe the Kirkwood-Alder fluid-solid transition [45] (see, the RDF plot for S=100 in Fig.1). For a
lower density n2, lower sound velocities and shorter cut-off radius rc , the value of S , in which DPD
fluid does not freeze, can be greater. Because the complex fluid has a longer length-scale
structure than the ordinary molecular scales, we may assume that the DPD particle system may
simulate a colloidal-like suspension with nanoscale colloidal particles, about 2nm. As shown in
[37, 39], defining the colloidal and solvent types of particles and assuming the short-ranged
attractive forces between them, we can obtain a more sophisticated ordering in colloidal fluids by
including the effects due to the creation of micelles and colloidal arrays.
Due to the limited range of the spatial-temporal scales of the system, the gravitational
acceleration G (or an other external shock) should be very large for observing the mixing process
for reasonable number of timesteps. The value of G assumed here is an order of magnitude
smaller than the value of G we used for the Rayleigh-Taylor mixing in the microscale (S=1) [268
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
9
2
28]. Therefore, we may expect that for S=10 mixing will take place for G values of order 10 m/s
2
- the range of acceleration commonly used in mixers. The model requires some modifications at
larger spatial scales, than those in our DPD model. For example, instead of using a linear weight
function, we can apply a more realistic types of (r), e.g., the first derivative of DLVO potential
used for colloids [45-46] or a Gaussian-like function (see [47]).
In spite of the problems dealing with matching the DPD model to the larger spatial scales,
we show below that in the present form this model can be used as a viable first approximation for
investigating the dynamic properties of mesoscopic fluids. From a computational point of view,
DPD approach has an advantage in that it is less demanding than using MD. In the simulation of
the system defined in Table 1 for a fixed spatio-temporal scale (e.g., 1 m in 10 ns) we find
astonishingly that MD needs roughly 10 times more memory and 100 times more CPU time than
DPD. These ratios can be even higher for greater values of S.
SIMULATION RESULTS
As shown in [26-28], there are basically four stages in the temporal evolution of the microscale RT instability in miscible Lennard-Jones fluids:
1.
The incipient regime, during which the system oscillates vertically after the gravity field
has been turned on. The horizontal perturbations are created due to the thermal
fluctuations.
2.
The perturbations, which correspond to the most unstable wavelength m , start growing
exponentially.
3.
The large scale regime enters in, when the effects from the initial conditions have been
2
erased. The thickness of the mixing layer increases as t .
4.
The final chaotic regime with many length-scales participating.
For immiscible and viscous R mixing of very thin fluid layers, the final chaotic regime may not
-T
be resolved, due to the small-scale coalescent process involving droplets. Below we report the
corresponding stages of immiscible mixing in mesoscale on the basis of simulation results
obtained by using the DPD model.
Incipient regime
In macroscopic continuum models, the R-T instability is initiated by introducing incipient 2
-D
harmonic perturbations (e.g. [8,9]). In the linear regime disturbances larger than a critical
9
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
wavelength crit grow exponentially [6] in the initial stages of the R mixing. For some types of
-T
initial perturbation, dominating wavelengths can depart from the characteristic one. We note that
in 3-D circumstances the fastest growing wavelength may be different [59]. This could lead to
characteristic mixing patterns in the final stages such as diapiric spacing close to the
characteristic wavelength [9].
On the other hand, in the world of discrete particles the R-T instability emerges
spontaneously [26-28] without any preference for wavelength selection. In the microscale, over
the nanometer length scale, both thermal fluctuations and the fluid elastic compressibility cause
principally the process of mixing. Because of the limited depth h of the computational box, initial
compression of lighter fluid induces the vertical oscillations of the particle system. Assuming that
the oscillations are harmonic, the period of these oscillations can be roughly estimated to be:
h
= ~
c
where
~
c
(15)
is the average sound velocity of the system. The energy of the one-dimensional
oscillations is transferred by the amplification of the thermal fluctuations.
This situation resembles the initiation of fingering observed in experiments from thin liquid
films flowing down an incline [36,48], where 2-D periodic adjacent wave fronts precede the 3
-D
instabilities. As show laboratory experiments [48], when entrance flow rate is perturbed at
adequate frequency by applying small sinusoidal pressure variations, after some start-up time
the straight 2 wave fronts brake spontaneously resulting in 3 wave production and then in
-D
-D
chaotic flow. The pressure fluctuations enable the energy of 2-D oscillation be transferred in
horizontal direction. We show that because of the large thermal fluctuations in the modeling of
thin films flowing down a slope, the 2-D waves are broken shortly after the initiation.
Together with the thermal fluctuations, the elastic compressibility constitutes the very
important factor impacting on the start-up time of the R fingering. For the high_c system, i.e.,
-T
the fluid with lesser elastic compressibility, where
~ =971m/s and h=224nm (h=42 in r units),
cH
c
the period of oscillation 0.72 ns, that is, about 3
600 timesteps. As shown in Figs.2a,b, the
value of for the high_c system is about 3 times shorter than for the low_c system. This comes
from the Eq.(15), because the sound velocity
than
~ =307m/s is approximately 3 (10) times less
cL
~ . Likewise for a computational box 4 times higher, increases proportionally according
cH
to Eq.(15) (compare Fig.2a and Fig.2c). This means that for the macroscopic systems with an
infinite depth h, the oscillations do not appear as . For finite values of h, due to the lack of a
physical mechanism for converting the vertical oscillations into horizontal perturbations in the
continuum models, they will not be able to excite the R-T instability.
10
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
As shown in Fig.3, the start-up time is distinctly shorter for more compressible fluids, i.e.,
for the low_c system. For this system, the initial one-dimensional oscillations of larger amplitudes
and longer oscillation period appear (see Eq.(15) and Fig.2) than for less compressible fluid
high_c . Therefore, there is a higher probability for the creation of the vertical perturbations, which
correspond to the most unstable wavelength m.
Linear regime
During the evolution of the R instability, vertical oscillations last for a long time, thus coupling
-T
nonlinearly with the large thermal fluctuations at the density interface. It is indeed difficult to
identify the linear regime in Fig.3, in which the perturbations corresponding to the most unstable
wavelength start to grow exponentially. In this regime, the capillary number is small. We define
here the capillary number to be the ratio of viscous forces H V/ v to capillary forces / l [49] ,
where V/ v is the characteristic shear rate, is the surface tension and l is the striation
thickness, respectively. The density interface wrinkles and small ripples grow into fingers. Their
size and separation depend, among other things [6-20], on the surface tension, viscosity and
elastic compressibility. Below we will estimate approximately whether some theoretical
predictions from the macroscopic model are valid for compressible, viscous DPD fluids in the
linear regime of the R-T mixing.
For incompressible, viscous fluids with the surface tension =0, the system becomes
unstable for disturbances at all wavelengths. For a macroscopic fluid the mode of maximum
instability appears at the wavelength [6-7]:
2
m 4
G
where:
=
H + L
H + L
1/ 3
(16)
is the average kinematic viscosity, H , L and H, L are respectively the
dynamic viscosities and densities of heavy and light fluids. Let us define additionally
Nm =
L xy
(17)
m
which is approximately equal to the number o fingers at the fluids interface at the beginning of
f
the R-T mixing. When the surface tension is taken into account, unstable modes occur only for
> crit where [6]:
crit = 2
?
[G ( H
L )]
(18)
11
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
and crit is independent on viscosity. For inviscid fluids m= 3crit . For viscous fluids m increases
with the surface tension and with fluid viscosity [6].
The surface tension at the density interface for DPD system depends on the
immiscibility of the two superimposed fluids. In the DPD system the immiscibility is defined by
value (see Eq.(14)). For =P HL-P H 0 (see Eqs.(10,14)), the fluids H and L are miscible, while for
>0 they are immiscible. From Fig.3, we can see that the speed of mixing decreases rapidly with
value both for the high_c and for the low_c systems. Therefore, we may expect the surface
tension at the interface between two immiscible fluids to be a strong function of . In order to
find () dependence we have employed the Laplace scaling law in the same way as in [33-35].
As shown in Fig.4, where we display two plots of
() for two different kinematic viscosities
(given in units) of the particle system, the surface tension is clearly a non-linear monotonically
increasing function of and . In particular, in the miscible case, (i.e., =0) the surface tension
is equal to 0 (see Fig.4).
Employing Eqs.(16-17) to the values from Table 1 for the high_c system with =0 and
with box aspect ratios f =Lxy /h= 2 and f =8, we obtain that Nm1=1.8 and 7.2, while for small surface
-11
tension 110
N (=0.03 N/m, see Fig.4) it comes out from Eqs.(18,17) that Nm2=2.2 and 8.9
for the respective aspect ratios. The values of Nm2 will be even smaller in the viscous case.
Because the Nm values are very small (m value is equal to the thickness of the L and H fluids)
and boundary effects influence the flow, this situation is very different from the macroscopic
models, in which the total thickness of the two fluid layers is taken to be infinite.
Surprisingly, the number of fingers (Figs.5,6) resulting from DPD simulations is very close
to the Nm values computed from the formulae derived from the continuum hydrodynamical
equations. The finite-thickness correction for mh reduces the effective Atwood number by only
less than 2% [12]. As one would expect, the number of fingers is greater in the more
compressible (low_c) case in which the surface tension is an order of magnitude smaller than the
high_c system [35] (see Figs.5,6).
Regimes with large-scale features and the breakup of microstructures
In the turbulent R instability, the sustained acceleration causes the dominant spatial scales to
-T
evolve self-similarly and the initial scales are forgotten. The mixing layer depth hT increases with
k
time as hTA TGt ( - the Atwood number (see Eq.(8)), G - gravitational acceleration, A
T
proportionality constant) as shown by many numerical [7,8,11] and experimental results [18-20].
For the two end-members:
1.
two fluid layers of equal depths,
12
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
2.
very thin heavy fluid layer placed at the surface of lighter fluid,
the values of k =2 and k=1 are obtained respectively. In [28] we show that for two superimposed
microscopic Lennard-Jones fluids of equal depth the value of k =2. Moreover, the proportionality
constant A B for bubble mixing layer contribution is approximately the same as these obtained in
laboratory experiments and macroscopic simulations. A similar scenario can be found in DPD
fluids. As shown in Fig.3a, the growth of interface length with time for two fluid layers of equal
2
depth can be well approximated by a + b t +ct polynomial in the period of time from the beginning
of the mushroom structures formation to the breakup of stems. For the cases without the erosion
2
by mushrooms the ( >0) the t dependence can be attributed to the propagation speed of the
mixing front. We can conclude from DPD simulation results depicted in Fig.3b that in the case of
thinner region with heavy fluid, the time exponent k changes from 2 to 1. Thus the same behavior
is found for both the macroscopic and mesoscopic DPD particle system.
The capillary number increases with time due to the increase of the shearing rate, which
also influences the viscosity. This shows the non-linear dependence of the viscosity on the
shearing rate. Because the box is very thin, growing fingers and the bubbles formed after fingers
breakup (Figs.5,6) saturate almost completely the lighter fluid from the bottom. As shown in Fig.
3, mixing is more vigorous for more compressible system and slows down with increasing surface
tension (and ). For fully miscible DPD fluids ( =0), due to the local decrease of viscosity and its
dependence on local shearing on the density interface, the fingers and mushroom structures of
lighter fluid erode forming plume like patterns (see Fig.5). For small values of the surface tension
(=0.03 N/m) in the low_c system, secondary fingers and thin needles appear instead of plume
6
(see Fig.5c, and Fig.6a). In [39] we inspected this phenomena for larger particle system (2 10
particles) and more compressible fluids. It turns out that the secondary fingers appear after
decompression of the lighter fluid, i.e., when the less viscous fluid pushes the more viscous fluid.
This dynamical effect is very similar to the Saffman-Taylor instability [50]. For the high_c system
and low_c fluid with greater surface tension, the bubble erosion occurs only locally in selected
places with low viscosity caused by the greatest rate of shear (see Fig.6b,c and Fig.7a).
The erosion patterns presented in Figs.5,6a change the flow radically in comparison with
situation depicted in Figs.6b,7. The kinetic energy of the light particles overpowers that for the
heavy phase (see Fig.8a). During miscible mixing ( =0) the fingers disappear and light particles
disperse completely in the heavy fluid (see Fig.5b). The density interface vanishes and chaotic
motion of particles can be observed. Other type of coupling between flow and breakup occurs for
the low_c system with small . The clustering of the thin fibers initiated by the secondary fingers
and needles into thicker threads shown in Fig.9, dissipate the energy from the system. The
coalescencing process becomes to be the dominant type of flow resulting in rapid overturn (see
Figs.6a).
13
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
We observe the fast mixing of the thin mesoscopic fluid layers at the density (light -heavy)
interface at the onset, which is caused by the erosion of the fingers and mushroom structures.
This is to be contrasted with the flow patterns shown in Figs.7,10, in which relatively stable
bubbles appear. The process of detachment of mushroom structures is responsible for the
change in flow type from growing fingers to the dynamics of free bubbles. For the high viscosity
cases (=100) the stems of mushroom structures are saturated. They become thinner and
approach a zero thickness without any breakup of the stem (see Figs.6,7a, 9,10). By decreasing
the viscosity of L and H liquids eight times ( =12.5), thus resulting also in more than twofold
increase of the surface tension for =0.03 N/m (see Fig.4), we find a radical change in
mechanism of stem disintegration. The four principal mechanisms, the same as those responsible
for droplets breakup in macroscale [51], can be observed in DPD simulation of the R instability.
-T
As shown in [1,51], moderately extended drops for capillary number close to a critical value
breakup by a necking mechanism due to a concentration of low viscosity from the high-rate of
shear. Due to the elongation flow (see Fig.11 - velocity field) in mushroom stems and p= L/ H 1
assumed, according to [1], the necking mechanism cause the breakup in Fig.7b. In the two
rectangles in Fig.7b (1 and 2, marked by thick line) one can see the two phases of this process.
After breakup (rectangle 1) two daughter droplets are created (rectangle 2) and the stem remnant
contracts producing the spherical droplet (rectangle 3).
Another mechanism responsible for the stem breakup is shown in the same figure in
dashed rectangle 4 in Fig.7b. The mushroom structure detaches thus producing from the stem tip
a cluster of small bubbles, which can overtake the large bubble (see pictures 1 and 3 in Fig.11, 3
is the reverse of 2 image). This process is similar to the tipstreaming mechanism, in which small
drops break off from the tips of moderately extended pointed drop [51]. The conditions at which
tipstreaming occurs are not well known [1]. From Fig.7b one can see that the stem is distinctly
larger and thicker than those in the necking case. Moreover, the remnant pieces do not relax
back so quickly to spherical droplet as in the case of necking.
A supercritically extended drop in the presence of low shear rate and surface tension
causes a breakup of the end-pinching type. When the drop is stretched to a highly extended
thread the microstructure arises from short-wavelength capillary instabilities . The thread becomes
unstable to small fluctuations and will eventually disintegrate into a number of drops consisting of
secondary droplets between the larger primary drops. Because the end-pinching and capillary
instabilities are only observed when [1]:
1.
the length of filament is more than 15 times the initial radius of the drop,
2.
the flow is halted,
we cannot found these instabilities in thin film flows and vigorous flow situations, as is in Fig.6,7.
In Fig.12 we show the zoomed-in fragments of snapshots from DPD simulations of the R-T
instability in a deeper computational box and with a heavy fluid layer (gray structures in Fig.12)
14
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
ten times thinner than the light phase. As shown in Fig.12, the three mechanisms: necking, endpitching and capillary instabilities can be observed on different portions of the same extended
thread. The same situation was encountered previously by Tjahjadi and Ottino [52] in laboratory
experiments.
Mixing can be severely influenced by viscosity stratification in the system, as shown by
the work of Ten et al. [53]. A large enough viscosity contrast would inhibit mixing [53,54]. Lower
contrast in viscosity p= L/ H results in much better mixing due to critical erosion of mushroom
structures (Fig.13a). Whereas, non-monotonic dependence of interfacial length as a function of
average viscosity of L (light-heavy) fluids for high_c system (see Fig.13b), reflects the change
-H
in the flow character for mushroom patterns.
Coalescence
Due to the extremely thin width of the box, after breakup of stems the flow speed after reaching a
maximum slows down (see Fig.8). The capillary number decreases and the two fluids overturn.
As shown in Fig.13, the density interface begins to shrink, due to the overturn and appearance of
coalescing microstructures.
In the course of mixing of immiscible fluids the coalescence process is competitive with
the tendency toward breaking-up of the microstructures. In the case of zero gravity, when the
immiscible fluids are initially mixed and is sufficiently small, the fluids separate and average
domain size R growths with time obeying algebraic scaling low R(t)=t. In 2-D space the value of
is equal to in the beginning and 2/3 for longer simulation times. As a result, the two fluids
separate completely after same time [33-35]. Otherwise (large value), at the onset, immiscible
fluids separate out into stable droplets of various sizes, thus creating emulsion-like substance
[35,39] (see Fig.14a).
Coalescence consists of three phases: collision, drainage of the thin liquid between the
two drops and the eventual rupture of the thin film [1]. The time for coalescence depends mostly
on the drainage time and the collision frequency depending on the dispersed phase volume
fraction.
In the case of no erosion, i.e., high viscosity and high surface tension (large , see
Fig.10), the overturning of the two fluids evolves smoothly without any microscopic breakup
morphology. However, the thin film of the heavy fluid persists for a long time. This is because the
drainage time decreases with the drop size [1] thus larger drops are less likely coalescing than
smaller droplets. The similar behavior is observed for the low viscosity and high surface tension
case (see Fig.7b), in which the breakup patterns occur but the bulk of matter is confined in the
bubbles, thus resulting in more than 90% overturn within a very short period of time. Because of
the coalescence of the larger structures and the agglomeration of the remaining droplets to the
15
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
bulk evolve slowly, the overturning takes very long time. It is even longer in the case of a higher
viscosity of the lighter fluid (see Fig.13a), which reduces greatly the mobility of the interface.
In the case of erosion the coalescing microstructures can be seen only for >0. As shown
in Fig.6a, in which the erosion is high, it occurs also for larger shearing rate resulting in clustering
of thin fibers into thicker threads. However, coalescence is more likely for a low capillary number
at the end of mixing process. In the places of the lower shearing rate the threads coalesce
trapping heavy fluids in "holes" in the bulk of light fluid (Fig.5c, Fig.6b). Simultaneously, the
droplets of light fluid coalesce in the bulk of the heavy fluid (Fig.6b, the last snapshot). Both
processes are of much finer granularity as compared to the cases without erosion and much
faster because the coalescing droplets are orders of magnitude smaller. Large erosion and then
fast coalescence result in rapid stop of vertical flow and e
mulsification of one or both fluids [1].
Thus a very long overturn time is observed. For the extreme cases, e.g., large differences in
density and elastic compressibility contrasts of the two fluids and a very large external field, the
creation of stable foam-like structures (see Fig.14b) may completely stop the overturning process.
CONCLUSIONS
In this work we have employed the dissipative particle dynamics approach to attack the problem
of immiscible Rayleigh-Taylor mixing at the mesosopic scale. The significant advances of DPD
can be delineated as follows:
1.
Conceptually the DPD method is simple and is based on fundamental physical principles of
Newtonian mechanics.
2.
Unlike the lattice-gas models [55,56], DPD particles evolve over continuum space in real
time, thus allowing for realistic visualization and understanding.
The method can produce quantitative results, when constructed within the context of cross-scale
systems [43, 57-58] in conjunction with microscopic (MD) or/and macroscopic particle models
(e.g., smoothed particle dynamics - SPH [47]).
The use of DPD in the simulation of complex fluids in the mesoscale has many distinct
advantages, such as
1.
Fluid granularity effects can be investigated, e.g., thermal fluctuations influence on incipient
of mixing,
2.
The microstructures are coupled with the macroscopic flow.
3.
Other processes, e.g. microscopic forces between dispersed phase and solvent, can be
included (see a two-level model presented in [37-39]).
4.
Sophisticated boundary conditions with complicated multi-phase interfaces can be easily
simulated.
16
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
The interactions between DPD particles influence strongly the speed of mixing. The surface
tension is the effect arising from the immiscibility of two superimposed fluids and is an
increasing function of . We show that in the case of high surface tension mainly larger
mushroom structures appear and carry the bulk of the fluid. The small microstructures do not
exert any influence on the flow. For lower viscosity, breakup patterns of mushroom stems appear.
One can distinguish in DPD simulation all four types of droplet breakup [1]: necking, tipstreaming,
end-pinching and capillary instabilities. This microscopic morphology changes the flow type from
rising mushrooms to bubble dynamics.
For small values of and viscosity contrast p= L/ H 1 the R-T mixing occurs due to the
erosion of larger mushroom structures. This morphological breakup, other than droplet breakup
mechanisms described in [1], produces fine grained mixing, which results in fast coalescence of
small droplets. The coalescing droplets dissipate energy from the system and decelerate the flow,
thus influencing the flow structure.
The DPD approach suffers also some drawbacks in the simulation of the R instability.
-T
As in the other discrete-particle methods, DPD takes a lot of computer time. The 2-D simulations
presented here were carried out for small particle systems, which correspond to the fluid samples
of 1 m in size. Less than 20 ns can be simulated for very thin liquid layers. For such systems the
gravitational acceleration needs to be very large.
The system can be also matched to the actual spatio-temporal scales by rescaling, i.e.
increasing scaling factor S. However, we show in [35] that matching DPD system to the spatiotemporal scale under interest is difficult. For relatively small scaling factor S and for given
physical properties of the liquid, the DPD system "freezes". We can easily "melt" it by assuming
unrealistically small sound velocities in the particle system. For non-equilibrium systems like the
R-T mixing, such an assumption may produce erroneous results. Therefore, we may surmise that
the DPD model involves the change of the shape of (r) weight function. For larger scales S we
suppose that it should take the form of a Gaussian function, the same as it is used for the SPH
method [47]. This is reasonable; since the DPD method can be treated as a special, mesoscopic
case of the SPH technique [32]. We have shown here that the DPD method can be successfully
employed for simulating Rayleigh-Taylor mixing in mesoscopic complex fluids with microscopic
ordering for scales greater than 10 angstroms.
Acknowledgments
Thanks are due to Dr Arkady Ten and Dr Dan Kroll from the Minnesota Supercomputing Institute
for inspiring discussions. Support for this work was provided by the Energy Research Laboratory
Technology Research Program of the Office of Energy Research of the U.S. Department of
17
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
Energy under subcontract from the Pacific Northwest National Laboratory and partly by the Polish
Committee for Scientific Research (KBN) Grant No. 8T11C00615.
REFERENCES
1.
Ottino, J.,M., De Roussel, P., Hansen, S., Khakhar, D.,V., Mixing and dispersion of
Viscous Liquids and Powdered Solids, Advances in Chemical Engineering, 25, 105-204,
2000.
2.
Ten, A., Yuen, D.A., Podladchikov, Yu.Yu.,
Larsen,
T.B., Pachepsky, E. and
A.V.Malvesky, Fractal features in mixing of non-Newtonian and Newtonian mantle
convection, Earth Planet. Sci. Lett., 146, 401-414, 1997.
3.
Ottino, J.M., The Kinematics of Mixing: Stretching, Chaos and Transport, 364 pp.,
Cambridge Univ. Press, 1989.
4.
Metcalfe, G., Bina, C.R.
and
J.
M.
Ottino,
Kinematic considerations for mantle
mixing, Geophys. Res. Lett, 22, 743-746, 1995.
5.
Schiesser, W.E. , Computational Mathematics in Engineering and Applied Science,
ODEs, DAES, and PDEs, 587pp, CRC press, 1994.
6.
Chandrasekhar, S., Hydrodynamic and Magnetohydrodynamic Stability, Oxford Press,
(1961).
7.
Sharp, D., H., An overview of Rayleigh-Taylor instability, Physica D, 12, 3-18, 1984.
8.
Youngs, D.L., Numerical simulation of turbulent mixing by Rayleigh-Taylor Instability,
Physica D, 12, 45-58, 1984.
9.
Schmeling, H., On the relation between initial conditions and the late stages of RayleighTaylor instabilities, Tectonophysics, 133, 65-80, 1987.
10. Mikaelian, K.O, Turbulent mixing generated by Rayleigh-Taylor and Richtmyer-Meshkov
instabilities, Physica D, 36, 343-357, 1989.
11. Deng, Y., Glimm, J., Sharp, D., H., Mixing and chaotic microstructure, Los Alamos
Science, 21, 124-132, 1993
12. Mikaelian, K.O, Rayleigh-Taylor instability in finite-thickness fluids with viscosity and
surface tension, Phys.Rev. E, 54, 4, 3676-3681, 1996.
13. Ortoleva, P.J., Geochemical Self-Organization, 412pp. Oxford Univ. Press, 1994.
18
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
14. 1996 Division of Plasma Physics 38th A
nnual Meeting, 11-15 November 1996, Denver,
Colorado, http://flux.aps.org/meetings/dppfl96.html
15. 1998
Division
of
Plasma
Physics
Meeting
Program,
November
16-20,
1998
New Orleans, LA http://www.aps.org/BAPSDPP98/
16. Mirin, A.A., Cohen, R.H., Curtis, B.C., Dannevik, W.P., Dimits, A.,M., Duchaineau, M.,A.,
Eliason, D.E., Schikore, D.R., Anderson, S.E., Porter, D.,H., Woodward, P.,R., Shieh,
L.,J., Very High Resolution Simulation of Compressible Turbulence on the IBM SP
System, Proc. of Supercomputing 99, Portland, Oregon, USA, Nov. 13-19, 1999
17. Oparin, A., Abarzhi, S., Three-dimensional bubbles in Rayleigh-Taylor instability, Physics
of Fluids, 11/11, 3306-3011, 1999
18. Read,K.,I., Experimental Investigation of Turbulent Mixing by Rayleigh-Taylor Instability,
Physica D, 12, 45-58, 1984
19. Voropayev, S.I., Afanasyev, Y.,D., van Heijst, G.,J.,F., Experiments on the evolution of
gravitational instability of an overturned initially stably stratified fluid, Phys. Fluids A, 10,
2461-2466, 1993
20. Dimonte, G. and M. Schneider, Density ratio dependence of Rayleigh-Taylor mixing for
sustained and impulsive acceleration histories, Phys. Fluids, 12/ 2, 304-321, 2000.
21. Glimm, J., Menikoff, R., Sharp, D.H. and Q. Zhang, A numerical study of bubble
interactions in Rayleigh-Taylor instability for compressible fluids, Phys. Fluids, A , 2,
2046-2054, 1990.
22. Glimm, J., Zhang, Q. and D.H. Sharp, The renormalization group dynamics of chaotic
mixing of unstable interfaces, Phys. Fluids,A 3, 1333-1335, 1991.
23. Unverdi, S.O. and G.
Tryggvason,
A
front-tracking
method
for viscous,
incompressible, multi-fluid flows, J. Comput. Phys., 100, 25-37, 1992.
24. Mosler, M., and Landman, U., Formation, Stability, and Breakup of Nanojets, Science,
289, 1165-1169, 2000.
25. Dzwinel, W., Alda, W., Kitowski, J., Moscinski, J., Wcislo R., and Yuen, D., Macro-scale
simulations using MD method, Molecular Simulation, 15, 343-360, 1995.
26. Moscinski, J., Alda, W., Bubak, M., Dzwinel, W., Kitowski, J., Pogoda, M., Yuen, D.,
Molecular Dynamics Simulations of Rayleigh-Taylor Instability, in: Annual Review of
Computational Physics, V , 97-136, 1997.
27. Alda, W., Dzwinel, W., Kitowski, J., Moscinski, J., Pogoda, M., Yuen, D.A, Complex FluidDynamical Phenomena Modeled by Large-Scale Molecular Dynamics Simulations,
Computers in Physics, 12, 6, 595-600, 1998.
19
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
28. Dzwinel, W., Alda, W., Pogoda, M., Yuen, D.A., Turbulent mixing in the microscale,
Physica D, 137, 157-171, 2000.
29. Hoogerbrugge,
P.J.,
Koelman,
J.M.V.A.,
Simulating
Microscopic
Hydrodynamic
Phenomena with Dissipative Particle Dynamics, Europhysics Letters, 19, 3, 155-160,
1992.
30. Espaol, P., Serrano, M., Zuniga, I., Coarse-graining of a fluid and its relation with
dissipative particle dynamics and smoothed particle dynamics, International Journal of
Modern Physics C, 8, 899, 1996.
31. Marsh, C., Backx, G., Ernst, M.H., Static and dynamic properties of dissipative particle
dynamics, Physical Review E, 56, 1976, 1997.
32. Espaol, P., Fluid particle model, Physical Review E, 57/3, 2930-2948, 1998.
33. Coveney, P., V., Novik, K., E., Computer simulations of domain growth and phase
separation
in
two-dimensional binary immiscible fluids using dissipative particle
dynamics, Physical Review E, 54(5), 5134-5141, 1996.
34. Novik, K.E., Coveney, P., V., Using Dissipative Particle Dynamics to Model Binary
Immiscible Fluids, International Journal of Modern Physics C, 8, 909-915, 1997.
35. Dzwinel, W., Yuen, D.,A., Matching Macroscopic Properties of Binary Fluid to Dissipative
Particle Dynamics Interactions, International Journal of Modern Physics C, 11/1, 1
-25,
2000.
36. Dzwinel, W., Yuen, D.A., Dissipative particle dynamics of the thin-film evoluation in
mesoscale, Molecular Simulation, 22, 369-395, 1999.
37. Dzwinel, W., Yuen, D., A., A Multi-level, Discrete Particle Model in Simulating Ordered
Colloidal Structures, Journal of Colloid and Interface Science, 225, 179-190, 2000.
38. Dzwinel, W., Yuen, D.,A., A two-level, discrete particle approach for large-scale
simulation of colloidal aggregates, International Journal of Modern Physics C, 11/5, 10371062, 2000.
39. Dzwinel,W., Yuen, D.,A., Large-Scale Simulations of Colloidal Dynamics With Fluid
Particle Model (FPM) and Dissipative Particle Dynamics (DPD), European Physical
Journal E-direct, (submitted October, 2000), www.msi.umn.edu/~dwitek/e-paper.html
40. Dzwinel, W., Bargiel, M., Kitowski, J., and Moscinski, J., Link -List and Method of Lights in
Molecular Dynamics Simulation Search for the Best Method of Forces Evaluation in
Sequential MD Codes, Molecular Simulation, 4, 229-239, 1989.
41. Hockney,R.,W., Eastwood, J.,W., Computer Simulation Using Particles, McGraw-Hill Inc.
1981
42. Espaol, P., Serrano, M., Dynamical regimes in DPD, Phys. Rev. E, 59/6, 6340-7, 1999.
20
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
43. Hadjiconstantinou,
N. ,G.,
Petera,
A.,
T.,
Heterogeneous
atomistic-continuum
representation for dense fluid systems, International Journal of Modern Physics C, 8, 4, ,
967-976, 1997.
44. Cahn, J.W. and J.E. Hillard, Free energy of nonuniform system. I. Interfacial free energy,
J. Chem. Phys., 28, 258-267, 1958.
45. Gast, A.,P., Russel, W.,B., Simple Ordering in Complex Fluids, Physics Today, 24-30,
December 1998.
46. Gelbart, W.M., Ben-Shaul, A., The New Science of Complex Fluids, J. Phys. Chem. 100,
13169-13189 1996.
47. Libersky L.D., Petschek, A.G. Carney, T.C. Hipp, J.R..Allahdadi, F.A, High Strain
Lagrangian Hydrodynamics, J. Comp. Phys., 109/1, 67-73, 1993.
48. Oron, A., Davis, H., S., Bankoff, S., G., Long-scale evolution of thin liquid films, Reviews
of Modern Physics, 69/3, 931-980, 1997.
49. Daoud, M. and Williams, C.E. (Eds), Soft Matter Physics, Springer Verlag, Berlin, 320
pp., (1999).
50. Saffman, G.T., Taylor, G.I., The penetration of a fluid into a porous medium or Hele-Shaw
cell containing a more viscous liquid, Proc. R. Soc. London, Ser, A., 245, 312, (1958)
51. Stone, H,A., Dynamics of drop deformation and breakup in viscous fluids, Ann.Revs.Fluid
Mech., 26, 1994, 65-102.
52. Tjahjadi, M., and Ottino, J.,M., Stretching and breakup of droplets in chaotic flows,
J.Fluid.Mech., 232, 1991, 191-219.
53. Ten, A.A., Podladchikov, Yu.Yu., Yuen, D.A., Larsen, T.B. and A.V. Malevsky,
Comparison of mixing properties in convection with the particle line method, Geophys.
Res. Lett., 25/16, 3205-3208, 1998.
54. Ten A.A., Yuen, D.A., Podladchikov Yu.,Yu., Visualization and Analysis of Mixing
Dynamical Properties in Convecting Systems with Different Rheologies., Electronic
Geosciences, 1999 4:1.
55. Zaleski, S. and Ph. Julien, Numerical simulation of Rayleigh-Taylor instability for single
and multiple salt diapairs, Tectonophys., 206, 55-69, 1992.
56. Rothman, D.,H., Zaleski, S., Lattice-Gas Cellular Automata: Simple Models of Complex
Hydrodynamics, Cambridge University Pres Cambridge, UK, (1997).
57. Dzwinel, W., Alda W., Kitowski, J., Yuen, D.A., Using discrete particles as a natural
solver in simulating multiple-scale phenomena, Molecular Simulation, 25, 6, 2000, 361384.
21
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
58. Dzwinel, W., Alda, W., Yuen, D.A., Cross-Scale Numerical Simulations Using Discrete
Particle Models, Molecular Simulation, 22, 397-418, 1999.
59. Kaus, B.J.P. and Yu. Yu. Podladchikov, Forward and reverse modeling of the threedimensional viscous Rayleigh-Taylor instability, Geophys. Res. Lett., in press, 2001.
22
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
RDF
MD L-J liquid
3
DPD S=25
DPD S=100
RDF
2
1
0
0.0
0.5
1.0
1.5
2.0
2.5
DISTANCE
Fig.1. Radial distribution functions for DPD fluid for two different values of S. Distance is given in units.
The R DF for the Lennard-Jones fluid (dashed line) is presented to show the contrast in the spatial scale
between L -J and DPD fluid for S=25.
23
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
b
a
c
0.24
PRES SURE [ N/ m]
0.23
0.22
0.21
0.20
0
40000
80000
120000
160000
TIME (*2.0E-13)
Fig.2. The variations of thermodynamic pressure (computed from the virial law [33-35]) with time for three
DPD simulations describing the R-T instability. Upper line (dark) corresponds to the L fluid, bottom one
(gray) to the H fluid. a) high_c fluid, b) low_c fluid, sound velocity c is about 3 times smaller than that for a),
c) high_c fluid but the box height H is 4 times greater than that for a) and b). The thick line in Figs a) and b)
is the moving average fit of the pressure.
24
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
a)
INC REA SE OF IN TERF AC E LEN GTH ( *2. 24E- 9 m)
5000
Delta=0.0
Sonic Speed
High
4000
Low
sqr(t) fit
3000
2000
Delta =0.013 N/m
1000
Delta =0.033 N/m
0
10000
20000
30000
40000
50000
60000
TIMESTEP
b)
2
10
Width of m ixi ng la y er
9
Thickness ratio (light fluid)
8
7
0.80 (power 1.06)
6
0.75 (power 1.64)
5
0.50 (power 2.02)
4
3
2
1
100
2
3
4
5
6 789
1000
2
3
4
5
T ime s tep
Fig.3. a) The temporal evolution of interface length between H and L fluids for high_c a nd low_c 2-D DPD
particle systems for various values of the immiscibility factor . The t2 regime is approximated by 2nd order
polynomials. b) The mixing layer growth with time for 3-D DPD particle system for various fluid thickness
ratios. The numbers on the left in the legend show the light fluid thickness in correspondence to the unit box
height. The numbers in parenth eses correspond to the fit of the time exponent k.
25
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
10
Viscosity
SU RF AC E TENSION ( *E-11 N )
8
lower ( Om ega =12.5)
higher (Omega = 100)
6
4
2
0
0.00
0.01
0.10
1.00
IMMISCIBLITY FACTOR (N/m)
Fig.4. The surface tension dependence on the immiscibility factor for the high_c fluid. For lower kinematic
viscosity ( =12.5) the surface tension is higher. The lines are drawn as guides for visual purposes.
26
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
a)
10 ns
6 ns
b)
4 ns
10 ns
C)
6 ns
10 ns
Fig.5. The snapshots from DPD simulation of the R-T mixing of two DPD fluids L and H a) high_c fluid (=0,
=100) b) low_c fluid (=0, =100) c) low_c fluid (=0.03, =100). The box dimensions in rc units are
84 84.
27
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
28
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
-9
T=3.2 10 sec.
b)
a)
-9
T=4.8 10 sec.
-9
T=6.4 10 sec.
-9
T=7.6 10 sec.
-8
T=1.0 10 sec.
T=1.2 10-8 sec.
c)
Fig.6. The snapshots from DPD simulation of the R-T mixing of two immiscible fluids L and H (=0.033 N/m,
=100) a) low_c and b) high_c systems. We emphasize the initiation of bubbles erosion by using two
rectangles. In c) the velocity field for the fragment (in black rectangle) of the last snapshot from b) is shown.
29
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
The highest shearing rate corresponds with the places of maximum bubble erosion. The box dimensions in
rc units are 84336.
30
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
-9
T=3.2 10 sec.
-9
T=6.4 10 sec.
a)
b)
-9
T=8.0 10 sec.
T=4.8 10-9 sec.
4
1
2
-8
T=1.0 10 sec.
-9
T=6.4 10 sec.
3
T=1.2 10-8 sec.
-9
T =8 10 sec.
Fig.7. The snapshots from DPD simulation of the R-T mixing of two immiscible fluids L and H a) low_c fluid
(=0.041 N/m, =100) b) high_c fluid (=0.033 N/m, =12.5).
2 00
b)
a)
K IN ETIC E NER GY (K )
1 60
L
H
H
L
H
1 20
80
L
40
0
2 0000
4000 0
60 000
80000
1000 00
TIMESTEP (*2.0E-13 sec.)
Fig.8. The evolution of the kinetic energy of DPD particle systems with time a) low_c, =0 b) low_c, =0.04
N/m c) high_c, =0.3 N/m.
1
2
3
31
Dzwinel, W. and D.A. Yuen, "Mixing Driven by Rayleigh -Taylor Instability in the Mesoscale Modeled with Dissipative Particle
Dynamics", International J. of Modern Physics C, Vol. 12, No. 1, 91-118, 2001.
Fig.9. The zoomed-in snapshots from simulations presented in Fig.6a. The clustering of the thin fibers
initiated by the secondary fingers and needles into thicker threads
10 ns
35 ns
15 ns
Fig.10. The snapshots from DPD simulation of the R -T mixing for a large immiscibility factor ( =0.31 N/m).
1
2
3
4
Fig.11. The zoomed-in snapshots from simulations presented in Fig.7b The mushroom stem severs and
produces a cluster of small bubbles, which can overtake the large bubble (see 1 and 3, 3 is the reverse of 2
image). The velocity field is shown for situation from picture 1.
32
end- pinching
necking
necking
capillary instabilities
Fig.12. Three breakup microstructures observed on different portions of the same extended thread obtained
from DPD simulation of turnover of thin heavy layer in the bulk of lighter fluid.
500
INC R EA SE OF I NTER FAC E LE NGTH ( *2.2 4E- 9 m)
IN C RE ASE OF IN TERFA C E LEN GTH ( *2. 24E -9 m)
2500
a)
2000
Viscosity ratio
0.6
2.2
1500
1000
500
b)
Omega
400
1 00
50
1 2.5
300
200
100
0
0
0
40000
80000
TI MESTEP
120000
0
40000
80000
120000
TIMESTEP
Fig.13. The temporal evolution of interface length between H and L flui ds for DPD system for different p and
contrasts.
33
Fig.14. The large differences in density and compressibility contrasts for the two DPD fluids, results in
emulsification (a) or creation of stable foam like structures (b).
34