7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Development of spatial instability waves in laminar twophase mixing layers
L. Cheung* and T.A. Zaki*
Department of Mechanical Engineering, Imperial College London, London, SW7 2AZ, UK
1.cheung@imperial.ac.uk and t.zaki@imperial.ac.uk
Keywords: twophase shear flow, interfacial instabilities, Parabolized Stability Equations
Abstract
In this work we examine the growth and evolution of interfacial instabilities in the context of spatially developing,
laminar, twophase mixing layers. Previous studies of interfacial instabilities using temporal, parallelflow linear
stability analyses have successfully identified the initial mechanisms responsible for the waves' growth. These studies
have also demonstrated the presence of mode competition (Yecko et al. 2002), where the dominant unstable mode
can vary between internal and interfacial modes, depending on the flow parameters. However, in the current work, we
find that the dynamics of twophase mixing layers can be sensitive to additional factors, such as the spreading of the
mean flow, finite amplitude effects, and nonlinear interactions. To investigate these effects, we develop a twophase
formulation based on the nonlinear Parabolized Stability Equations (PSE). This formulation includes nonparallel
effects, nonlinear modal interactions, a coupled mean flow correction, and finite amplitude deformation of the
interface. Using this approach, we investigate mode competition, and find that the streamwise evolution of the flow
can significantly alter the predictions of locallyparallel, linear stability analyses. Nonlinear calculations of two and
threedimensional instability waves also show that nonlinear effects can both limit the growth of instability waves,
and simultaneously destabilize highfrequency, linearlystable modes.
Introduction
The spatial stability in twophase shear flows is a
paramount concern to many technological and industrial
applications, such as capillary jet breakup, wall films in
airfoil deicing systems, and twophase pipe and chan
nel flows. The majority of previous works have stud
ied the growth and development of the interface through
either Direct Numerical Simulations (DNS) or weakly
nonlinear, temporal stability analysis, but many nonlin
ear mechanisms have yet to be extensively investigated.
In the current work, we employ a nonlinear formulation
based on the Parabolized Stability Equations (PSE) in
order to study the evolution of disturbance waves in spa
tially developing, twofluid, laminar mixing layers. With
this framework we are able to show how nonlinearity
and mean flow modifications can significantly alter the
growth of instability modes.
The initial instability mechanisms in twofluid shear
flows have been previously explained through parallel
flow, Linear Stability Theory (LST). The work of Miles
(1957, 1960) identified the critical layer mechanism,
while Yih (1967) uncovered the long wavelength inter
facial instability due to viscosity stratification. Short
wavelength instabilities have also been investigated by
Hooper and Boyd (1983); Hinch (1984); Renardy and
Joseph (1985a,b). Within the same temporal LST frame
work, Yecko et al. (2002) discuss the possibility of mode
competition between different types of instabilities in
twophase mixing layers. The potential for energy am
plification due to nonmodal effects was also investi
gated by Yecko and Zaleski (2005) using the method
ology developed by South and Hooper (1999) and van
Noorden et al. (1998) for bounded twofluid flows.
To explore the evolution of these twophase instabil
ity waves beyond the early stage development regime,
a number of weakly nonlinear theories for twophase
flows have been proposed (Hooper and Grimshaw 1985;
Renardy 1989; Renardy and Renardy 1993; King and
McCready 2000). The majority of these formulations
are based upon the StuartLandau, GinzburgLandau, or
KuramotoSivashinsky equations. For marginally unsta
ble instability waves near the bifurcation point, Stuart
Landau and GinzburgLandau theories allow an insta
bility wave amplitude equation to be developed which
accounts for a limited set of nonlinear interactions.
King and McCready built upon the StuartLandau the
ory of Sangalli et al. (1995, 1997) to include cubic
and quadratic modal interactions in twofluid, bounded
flows. As a result, they observed instability waves
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
reaching steadystate amplitudes due to energy transfer
to higher frequency linearly stable modes. For small
wavenumber, longwave disturbances, the Kuramoto
Sivashinsky equation allows the temporal evolution of
the wave amplitude to be captured. Formulations based
on this theory have been used to study the effect of har
monic forcing in channel flows (Renardy and Renardy
1993), and uncovered the existence of a mean flow shift.
Nonlinear instabilities in twophase flows can also be
modelled through Direct Numerical Simulations (DNS)
of the NavierStokes equations. Simulations by Lom
bardi et al. (1996); Fulgosi et al. (2003); Alexakis et al.
I'"114); Boeck et al. (2007); Li et al. (1998); Li and
Renardy (2000) have shown that DNS can track the
evolution of the disturbance wave from its initial gen
eration to ligament and droplet formation. Different
methods have also been developed to accurately capture
complex deformations of the interface, ranging from
volumeoffluid (Gueyffier et al. 1999), level set (Suss
man et al. 1994; Herrmann 2008), or diffuse interface
(Ding et al. 2007) methods. Challenges remain, how
ever, when scaling DNS to realworld twophase prob
lems (Gorokhovski and Herrmann 2008). Due to their
computational cost, high fidelity direct calculations are
used relatively sparingly in the study of interfacial per
turbations.
Therefore, stability theory has continued to provide
valuable insight into the evolution and dynamical mech
anisms of twofluid shear flows. With the appropriate
formulation, such insight can be obtained in a much
more computationally efficient manner than direct cal
culations of the NavierStokes equations. Hence, the
aim of this work is to develop an efficient and accu
rate nonlinear instability wave formulation to study fi
nite amplitude disturbances in spatially developing, two
phase shear flows. Our approach extends the PSE frame
work introduced by (Bertolotti et al. 1992), and which
has since been applied to complex geometry flows (Her
bert 1997), aeroacoustics (Malik and Chang 2000), and
combustion (Day et al. 2001). The nonlinear twophase
PSE formulation described herein provides a means to
track the evolution of instability waves and the interface.
It also consistently incorporates the effects of nonparal
lelism, nonlinearity, and mean flow distortions.
Our work seeks to address several open questions re
garding the effects of nonparallelism and nonlinearity
in the context of laminar twofluid mixing layers. For
instance, the phenomena of mode competition between
liquid and interfacial instabilities has been studied by
Yecko et al. (2002) using parallel linear stability analy
sis. The nature of this competition is unknown in a spa
tially developing flow and is herein investigated. Also,
by comparing the results from both linear and nonlin
ear stability computations, we demonstrate the effects of
f (( t)
U2, P2' 12
(a) Physical configuration
)7
f($)
  :
,w 
r7=0
(b) Mapped configuration
Figure 1: Schematic of a twofluid, incompressible mix
ing layer in the (a) physical coordinate system and (b)
transformed coordinates.
modal interactions and finite interface deformations on
the growth of instability waves, as well as the modifi
cation of the mean flow due to finite amplitude distur
bances.
Methodology
Governing Equations
In this current study, we consider twodimensional,
incompressible laminar mixing layers of the type shown
in figure 1. An upper stream of fluid 1, with density pi,
viscosity p1, and velocity U1, flows over a lower stream
of fluid 2, with properties P2, P2, and U2. All variables
used in this study have been nondimensionalized using
U1 as the velocity scale, the initial vorticity thickness
60 as the length scale, and the material properties of the
upper stream. Thus, pi = 1 1, and the lower stream
density and viscosity have the values p, p2/1p and
Pr = 2/H1, respectively.
The NavierStokes equations governing the velocities
u(x, t) and pressure P of such a mixing layer can be
written as
V u 0
Du P V2U
P DLu VP + Re Vu
DtRe
where the Reynolds number Re = iU15o/pi. At the
interfacial location y f(x, z, t), the velocities u and
stress tensor 7 must satisfy the jump conditions
Eu f = 0, Eu it 0 (2a)
IT. IPI = P ancfi, IT. t = 0 (2b)
where f and t are the unit normal and tangential vectors
to the interface, and *I* denotes a jump across the fluid
boundary (Eq qly=f+ ql =f). In equation (2b),
a is defined as the surface tension coefficient and K the
interfacial curvature, and the Weber number can be writ
ten as We = piU So/a6 In addition, the location of the
interface itself is determined by the evolution equation
For mixing layers which experience no change in in
terfacial topology, and where f(x, z, t) remains single
valued, equations (1) and (2) completely determine the
solution to the problem. However, a direct solution of
the NavierStokes equations for spatially evolving flows
can be computationally expensive for the spatially devel
oping flows under consideration. Instead, this investiga
tion focuses on modelling the interfacial behavior using
linear and nonlinear instability waves. The problem for
mulation is then separated into two components: (1) de
termination of the timeindependent mean flow, and (2)
solving for the time dependent perturbations.
For this to occur, we assume that the flow variables
p [u v w P]T can be decomposed into mean (y) and
perturbation (y) components
p = p(x, y) + ((x, y,z, t). (3)
As explained below, the mean flow variables p are deter
mined from the forced boundary layer equations, while
the perturbation components p are found through the
PSE formulation. The solution of the two components
are then coupled via nonlinear forcing terms. Through
the use of a mapped coordinate system, large amplitude
interface deformations are also handled consistently be
tween the mean flow and the perturbation solutions.
Mean Flow Equations
The choice of the mean flow profile is critical when
considering the effects of the spreading mixing layer,
nonlinear interactions, and mean interface distortions.
Simple choices, such as analytic profiles or the self
similar mixing layer, fail to account for the mean flow
distortion arising from finite amplitude perturbations.
As discussed in this section, a more robust method to
calculate the mean flow is through a direct solution of
the forced boundary layer equations.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The governing equations for the mean flow compo
nents p can be derived by substituting the decomposi
tion (3) into the NavierStokes equations (1). Applying
a time average to the equations eliminates terms which
are first order in p. The equations can be further sim
plified by using the standard boundary layer approxima
tions, resulting in:
OU OV
+ = 0o
OiU i9U 02U
PU +pVy Re )y
PUxPVDy Re Dy2
(4a)
Fx,BL (4b)
Df (4c)
Dt
The nonlinear term Fx,BL on the right hand side of
(4b) arises from quadratic and higher order terms in
y, which can be of finite amplitude and can possess
a nonzero mean. In the free stream, the velocities
U(y oc) U1 and U(y oc) U2. The initial
inlet profile for U(y) and the upper and lower boundary
conditions for V are found from the selfsimilar solution
of the twofluid mixing layer, as shown in White (1991).
A consistent set of interfacial conditions can be found
for the mean flow variables by taking the time average
of equations (2).
U O U t 0
T1. f= P o f acrf, 17 t = 0 (5b)
Equations (5) are to be enforced at the mean interface
location y = f(x, z), which is also time independent.
All variables, such as 7 and K, are also assumed to sepa
rable into mean and fluctuating components. The slowly
evolving nature of the spatial mixing layer implies the
mean curvature K remains small, and the mean normal
of the interface remains close to n (0, 1, 0) both of
these assumptions can be verified a posteriori.
After applying the condition that IP = 0, and the
boundary layer approximations on V, the final results
for the mean interfacial conditions at y f are
Uo =0, oV =0
' 0, 0. (6b)
Equations (4) with conditions (6) are advanced down
stream using a CrankNicolson scheme with fourth or
der central differencing in y. The downstream marching
of the 9 variables is also synchronized with the solu
tion of the 9 variables, in order to exchange information
through the nonlinear forcing term F,,BL 
The evolution of the streamwise velocity profile U(y)
for a typical mixing layer, with parameters U2/UI =
U(x,y) 1 (a)
0.8 /
0.6 /
0.4
0.2
0
2 0 2
4
0(x) (b)
3
2
1
20
x
Figure 2: (a) The mean velocity profile U(x, y) at dif
ferent streamwise positions: x 0 ( ), x = 25
(), and x 50 ( ). (b) The momentum thick
ness O(x). The mixing layer parameters are given by
U2/Ui 0.03, p, 10, and p 2.
0.03, p, = 10, and p, = 2, is illustrated in figure 2. The
width of the mixing layer, as defined by the momentum
thickness
U(y) U2 dy
is also shown. The typical spreading rate for mixing lay
ers used in this study was small, generally on the order
of dO/dx r 0(102).
Perturbation Equations
By inserting the using the decomposition (3) in
the NavierStokes equations and removing the time
independent base flow components o, the equations gov
erning the perturbationvariables [iu v w4 p] can be
written as
V u f 0, (7a)
S+ (V) U + (U.V V) i (7b)
aOt
+P V2fu
Re
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and the perturbation interface equation is
Df
Dt (7c)
Dt
Similarly, the perturbation interface conditions can be
expressed as
ui fi = 0, fi t] t 0 (8a)
In i IP5 ft kain, rI t] 0 (8b)
In contrast to equations (2), these interfacial conditions
must be enforced at the true, time dependent interface
location y f(x, z, t). Enforcing these conditions at
the true interface location is essential to capturing the
effect of large amplitude, nonlinear deformations of the
interface. However, imposing the conditions at the time
dependent wave height f (x, z, t) presents additional dif
ficulties when attempting to use a frequency domain rep
resentation for the disturbance components p. These dif
ficulties can be overcome through the use of a mapped
coordinate system, as explained in the following section.
Mapped Coordinate System
Prior efforts at modelling instability waves in two
phase shear flows have generally applied interfacial con
ditions about a linearized location at y 0 While
this allows the interfacial conditions for each individ
ual mode to be enforced at a fixed location, the actual
position of the interface is timevarying and dependent
on the local velocities. However, by using a coordinate
transformation (9) to map the physical domain to a com
putational domain where the interface remains static,
we can avoid these difficulties in enforcing the interface
conditions.
T t (9)
y ,x
In the mapped coordinates (E,r, ), the interface is per
manently located at r f( () (see figure 1), and the
coordinate mapping function H satisfies the condition
H(Q, r = f, C, T) = f(, z, t). Therefore, the instanta
neous interface location y = f + f is always mapped
to r f. Note that this coordinate transform does
not allow for a multivalued interface deformation 
f(x, z, t) and H((, r, C, T) are required to be single
valued. Hence, although large amplitude interface de
formations can occur, the coordinate mapping precludes
the twophase system from cresting or undergoing vor
tex rollup.
From the coordinate transformation given by (9), the
derivatives in physical space are transformed according
to the following expressions:
oo U1y)U2 1
Pfi Vuf,
ax a
0t 0T + ag
0y 09 0
9y 0 + 7
az a( (a
where the coefficients are given by
dH OH
1 Q )1 Qr
1 
g1 +t 1, g1 1 aH
a? 1 O 1 1 + aT
Using the derivative transforms (10) in equations (7)
allows us to write the perturbation equations in a coor
dinate system following the interfacial wave, and also
eliminates the need to carry a separate equation (7c) to
track the location of the interface.
PSE Formulation
The spatially developing, twophase mixing layer pre
viously described is susceptible to convectively unstable
interfacial disturbances. Once forced at the inlet, these
disturbances travel downstream and can grow exponen
tially before saturating. For a mean flows which evolve
sufficiently slowly, it is possible to reduce the pertur
bation equations (7) to a parabolized system (Bertolotti
et al. 1992). This parabolized system of equations can
then be used to track the evolution of the interfacial dis
turbances downstream from the inlet. The PSE approach
has been previously applied to singlephase mixing lay
ers, jets, and boundary layers, and is extended here to the
twophase problem. The methodology presented below
is analogous to the PSE formulation used for compress
ible flows by Cheung and Lele (2009).
In deriving the governing equations for the pertur
bation quantities q(4, y, (, T), we assume the following
normal mode representation in the mapped coordinates
M,N
= mn($,)Amn() exp {i/3C im},
m,n
(10)
where the amplitude factor A, () is written as
An mn exp \1 acmn (')d'} (11)
This normal mode assumption makes use of shape func
tions Pmn and wavenumbers amn~() which are func
tions of the streamwise coordinate , as well as the span
wise wavenumber f3 and temporal frequency w m. In
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
addition, an initial amplitude ,n, is given for each shape
function ,mn.
The coordinate mapping function H((, r, (, T) can be
expanded in a similar manner, using the modal expan
sion
M,N
HQE,r,TI)= I Emnhmn(y) X
m,n 0
exp i amn(') d$' + iz iWmT }
The frequencies cm and wavenumbers amn, f, are cho
sen to correspond to those in equation (10), and the func
tions hmn (r) are chosen such that 0h/dy = 0 at r f,
and h(r) 0 as r oo. For the purposes of this
formulation, hmn (y) are taken to be a set of smooth
functions in r, which decay exponentially away from the
interface at the same rate as the linear stability eigen
modes. To satisfy the kinematic condition (7c) at the
interface, we also require that
h. (= f(
Umn (f)
i(amnU(f) Wim)
where 0mn is the vertical velocity shape function given
by the PSE solution.
Although ,mn and am, are allowed to evolve in the
streamwise coordinate, we assume that they vary slowly
relative to the typical instability wavelength. This trans
lates to the following assumptions on the PSE variables
a2 a2
0 mn() < 0(1), and a mn() < 0(1),
(12)
such that the second order derivatives in can be ne
glected. Additionally, the mean flow is also assumed to
vary slowly, meaning that
1 0
Re a< 0(1),
Re aE
a 20
and < 0(1),
8(2
and these derivatives are neglected in the analysis as
well.
Inserting equations (1011) and the assumptions (12
13) into the perturbation equations (7) yields the govern
ing PSE system in the mapped coordinates. Expressed
in operator form, the Nonlinear PSE (NPSE) system is
mn{ I .mn [mn/Awn .., (14)
where Fmn denotes the nonlinear terms arising from
higher order products of the disturbances p and the map
ping function H(Q, r, (, T). The linear operator mn can
be written as
mn
iwmG + iamnA + B
09Tj
02 a
+C + i3D E+M
+ tm ,
using the matrices A, B, C, D, E, G, M and N de
fined in the appendix. The coupling between the pertur
bation equations and the mean flow is supplied through
the zerofrequency forcing terms, i.e., FX,BL = F,oo.
The Linear PSE (LPSE) formulation can be recovered
by setting the right hand side of (14) to zero and elimi
nating the forcing term Fn,BL from equation (4b).
The PSE system (14) is supplemented by an addi
tional normalization condition on .mn. This normal
ization condition removes the streamwise ambiguity in
herent in the representation (10). Because both 6mn
and am
physical
To close
gral noi
0o
Addition
function
be expr
where k
free str
and at
shape f
allel flo
The
tion eqi
(4) is s
Day et
equation
tral diff
direction
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
08
06
04
02
S\ Liquid Mode
Interfacial Mode
0 02 04 06 08
0 02 04 06 08
12 14 16
Figure 3: LST computed growth rates a, for the liquid
and interfacial modes of a mixing layer with U2/UI
0.03, pr 10, and p, 2, at the streamwise positions
x 0 (solid line ), x 5 (dashed line ),
x = 10 (dashdot line ).
Mode competition
Share dependent on , streamwise changes to the To illustrate the important effects present in a nonpar
d variable p can be absorbed into either variable. allel, slowly diverging mean flow, we first consider the
e the PSE system, we impose the following inte streamwise evolution of growth rates for different mode
rm: types in a twophase mixing layer. Prior studies by
Yecko et al. (2002) have emphasized the distinction be
(+ ,,, , d7 ) \ tween various viscous modes in twophase mixing lay
m D oDt 0 ers. They found that the most unstable mode can be ei
(16) their an TollmienSchlichting (TS) type liquid mode or
nally, the interfacial conditions (8) on the shape an interfacial mode, depending on the parameters of the
as p are imposed at the location j f, and can flow. The results of locally parallel, linear stability anal
essed in the frequency domain as yses have shown that at low frequencies, the TS liquid
modes generally dominate over the interfacial modes,
0, .. 0, 0, (1a) while the reverse is true at higher frequencies.
However, local parallelflow analyses only provide an
incomplete picture regarding the mode competition be
... ., tween liquid and interfacial modes. The growth rates
Pmn 2. W (17b) obtained by LST are valid only for a particular stream
2 J WeM wise position, whereas in reality, the spreading of the
'( i,,... ) 0, mean flow will cause these growth rates to evolve down
Siavmn) 0, (17c) stream. As shown in figure 3, the growth rate of the most
S 2 unstable frequency decreases quickly with downstream
+ P3, and M = mnU(f) w. In the distance for liquid modes. Over the same distance, LST
eam, the boundary conditions on mn, are predicts that the decrease to the maximum growth rate is
;, 0, as oc, (18) less severe for interfacial modes, which makes the out
come of mode competition unclear as the mixing layer
the inlet 0, the initial conditions for the evolves. A better description of mode competition can
unctions and wavenumbers are provided by par be supplied through a PSE analysis of the liquid and in
w Linear Stability Theory (LST). terfacial mode behavior, which will capture the effects
computational methodology for solving perturba from mean flow spreading and nonparallelism.
nations (14) and (16), and mean flow equations In the current study, we consider the growth rates of
similar to the numerical algorithms described in the interfacial and liquid modes of a mixing layer with
al. (2001) and Cheung and Lele (2009). The PSE the parameters U2/U1 = 0.03, pr 10, pr = 2,
ns (14) were discretized using fourth order cen We = 300, and Re=150, whose dispersion relation is
'erencing in r, and advanced in the streamwise given in figure 3. By keeping the wave amplitudes rel
n using the first order implicit Euler method. atively small and limiting the domain to the initial de
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
aix (.3 (a)
0.2
0.1
0 10 20 30
0 10 20 30
G(x)
LU
10 / 0 1
0
0 10 20 30
x
10 20 30
Figure 4: Mode competition for the mixing layer shown
in figure 3. (a) Growth rate ai(x) and (b) ampli
tude gain Glo(x) for the interfacial mode (dashed line
) and liquid mode (solid line ) at frequency
S= 0.75. Thick lines correspond to predictions from
the locally parallel LST, thin lines from the LPSE sim
ulations. The circles (o) correspond to the initial condi
tions provided by LST
velopment region of the mixing layer (where x < 30),
the linear formulation suffices for the current example.
In addition, moderately unstable instability waves are
used in this part of the study. Although the most unsta
ble, lowfrequency liquid mode could be examined using
LPSE, the relatively high gain and fast growth rate intro
duces the possibility of finite amplitude effects. Thus,
we delay the discussion of nonlinearity and the behav
ior further downstream until the following section, and
concentrate on mode competition below.
First, we consider the evolution of an interfacial mode
and a liquid mode which are equally unstable (near the
frequency c 0.75 in figure 3) at the inlet, according
to locally parallel LST predictions. Once in the LPSE
simulations, however, the initial guesses for the growth
rate a, are adjusted to accommodate for the nonpar
allel effects and the presence of a streamwise varying
mean flow, resulting in the interfacial mode being al
most twice as unstable. As the LPSE modes evolve fur
ther downstream, we find that the growth rates for both
Figure 5: Mode competition for the mixing layer shown
in figure 3. (a) Growth rate ai(x) and (b) amplitude
gain Glo(x) for the interfacial mode (dashed line ),
liquid mode (solid line ) at frequency c = 0.50.
Thick lines correspond to predictions from the locally
parallel LST, thin lines from the LPSE simulations. The
circles (o) correspond to the initial conditions provided
by LST.
modes slowly decrease, but whereas the liquid mode
becomes stable near x 10, the interfacial remains
unstable (see figure 4a). The rapid decay and satu
ration of the liquid mode results leads to gain factor
Gmn(x) = \A,(x) /A,,(0) which is several orders
of magnitude smaller compared to the interfacial mode
(figure 4b).
In a second example we again compare the behavior
of the interfacial and liquid modes, but at a lower fre
quency (c = 0.50) where LST predicts the liquid mode
to dominate over the interfacial mode according to the
dispersion relation shown in figure 3. The correspond
ing PSE predictions is shown in figure 5. We can see
that the growth rates for both modes are adjusted once
in the LPSE simulation, and the liquid mode does re
main slightly dominant until approximately x 3. The
liquid mode continues to decay beyond this point, how
ever, and the interfacial mode eventually overtakes the
liquid mode downstream.
When comparing the streamwise evolution to the lo
G(x) (b)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
cally parallel LST results, we find that the presence of
nonparallel generally depresses the growth rate for liq
uid modes, and increases the growth rate for interfacial
mode. This is seen in both the low frequency and high
frequency examples mentioned above. Taking these ob
servations into account, we conclude that while the pre
dictions of locally parallel LST may be valid at a partic
ular location of the mixing layer, capturing nonparallel
effects and the streamwise evolution of the mean flow is
essential to determining the downstream behavior of an
interfacial wave.
2
1 5 (a) (b)
05
Y0
1 5
2L
160
0 25 (c)
02
015
01
005
165 170 175 180 185
05 1
0
100 110 120 130 140 150 160 170 180 190
Figure 6: (a) Instantaneous perturbation vorticity con
tours for the interfacial mode forcing in the p, = 10,
p, = 2 mixing layer. (b) The streamwise velocity shape
function Ilao(y) at x = 189. (c) The maximum inter
facial wave height max{f (x, t)}, as computed by the
nonlinear ( ) and linear () simulations.
Linear and Nonlinear Computations
In the following section we compare and contrast the
effects of nonlinearity on both liquid and interfacial
modes in twophase laminar mixing layers. The results
of two mixing layers are considered. The first mixing
layer is similar to case shown in the previous section,
(p, = 10, p, = 2, We = 300) with a higher veloc
ity ratio U2/U1 = 0.30 in order to hasten the onset
of nonlinear effects. The Reynolds number based on
the velocity shear ReAu pi(Ui U2)6o/p1 is de
fined to be ReA = 100, which causes nonparallel ef
fects to become evident early on in the mixing layer. At
15 (a)
y
* J ) % ') Ij )
100 105 110
115 120 125 130
X
02 04
Itiol
0.4
0.2 
0.1
40 60 80 100 120
Figure 7: (a) Instantaneous perturbation vorticity con
tours for the liquid mode forcing in the p, 10, p, f 2
mixing layer. (b) The streamwise velocity shape func
tion I io(y)l at x 131. (c) The maximum interfacial
wave height max {f(x, t)}, as computed by the nonlin
ear ( ) and linear () simulations.
higher Reynolds numbers, the basic nonlinear mecha
nisms remain unchanged, but the extent of the linear re
gion of the flow is increased, resulting in higher ampli
tude gain before the onset of nonlinearity. This case is
used in the study of nonlinearity in the development of
twodimensional interfacial and liquid modes. The di
rect calculations of Boeck et al. (2007) also studied low
Reynolds number mixing layers with similar parameters
(p, = 10, We z 36 360) and noted that the interfa
cial deformation in such cases eventually lead to more
complex ligature and droplet formation.
For spatially evolving mixing layers, two mechanisms
can lead to the saturation of the low frequency funda
mental mode. First, the spread of the mean flow will
naturally dampen the growth of the mean flow, eventu
ally causing it to become neutral and saturate. However,
the fundamental mode can also quickly transition to a
neutral or decaying wave as a result of the aforemen
tioned interaction with higher harmonics. A compari
son between the linear simulations (which capture the
first mechanism) and nonlinear simulations (which cap
ture both) will help to distinguish the effects between the
two.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
60 80 100120140160180
x
10 '
50
100 150
Figure 8: Streamwise development of the integrated
modal energy E,,n (x) for the interfacial modes in the
p, 10, p, 2 mixing layer. The different line styles
correspond to nonlinear ( ), linear (), and LST
( ) computations.
Two Dimensional Mixing Layers
We first examine the effects of nonlinearity on the in
terfacial wave growth in the twodimensional, p, = 10,
p, 2 mixing layer. For the interfacial mode, the in
stantaneous interfacial location f((x, t) and perturbation
vorticity for the mixing layer is plotted in figure 6a. Im
mediately adjacent to the vorticity contours we also plot
the streamwise velocity shape function alo at the posi
tion x 189 and corresponding instant in time (figure
6b). Note that the discontinuity in the shape function
occurs at the physical position of the interface and not
at the nominal interface height y 0 In figure 6c we
compare the growth of the interface between the non
linear and linear simulations. During the initial region
x < 140, the growth of the interface shows little differ
ence between the nonlinear and linear cases. However,
further downstream the enhanced growth of the nonlin
ear interfacial modes is visible.
An even more dramatic difference is seen in the mix
ing layer forced by liquid instability modes, where the
presence of nonlinearity acts to suppress the growth of
the interface. A corresponding set of figures for the liq
uid modes is shown in figure 7, using the same layout as
105 o
E20
10 10o
0 50 100
105
(b)
10 10 /
1015
1020
0
50 100
x
Figure 9: Same as in figure 8, for the liquid modes in
the p, 10, p, 2 mixing layer
figure 6. Note that the spatial extent of the streamwise
shape function alo is larger in the liquid case than the
interfacial case. As seen in figure 7c, the initial growth
of the interface is largely linear for x <60, but eventu
ally nonlinear effects limit the maximum amplitude of
the interface to f(x, t) z 0.2. Purely linear simulations
predict that the liquid modes continue to grow exponen
tially and will ultimately reach unphysically large inter
face deformations. As a result of this particular exam
ple, we can see how nonlinearity plays a major role in
the saturation mechanism of instability waves.
A better understanding of the interfacial wave behav
ior in this twophase mixing layer can be found by ex
amining the modal energy development of individual in
stability modes. For the twophase mixing layers in this
study, we define the integrated modal energy for a par
ticular shape function ',, m ,,n (, q)A,, (() as
E.(x)= p ( '... 2*. .2+ w' ) dn.
(19)
Figure 8 presents the modal energy behavior for the
fundamental and first harmonic interfacial modes in the
p, = 10, p, = 2 mixing layer, as calculated by the
NPSE, LPSE, and parallel flow LST. The largest differ
ences are seen in the behavior of the fundamental mode
in simulations with and without nonparallel effects. As
mentioned in the earlier discussion on mode compe
tition, the presence of nonparallel effects and stream
wise changes greatly enhances the growth of the funda
mental mode, and as a consequence, the LST computa
tions vastly underpredict the energy content in Eo10 when
compared to the PSE computations.
The inclusion of the nonlinear terms leads to two vis
ible effects in the development of the interfacial modes.
The growth of the fundamental mode is slightly en
hanced, and nonlinearity also excites higher harmon
ics which were previously stable. For example, during
the linear region of the flow, the second harmonic E30
decays exponentially, but is subsequently excited after
x > 150 (figure 8b).
The effect of the nonlinear interactions on higher fre
quency instability modes can also be seen through a sim
ilar analysis of the modal energy for the liquid modes in
the same mixing layer. In the initial development region
of the mixing layer, very little difference is seen between
the three formulations for the fundamental liquid mode
(figure 9a). However, further downstream of the inlet
the effects of the nonparallel flow become apparent as
the LST calculations of the modal energy diverge from
the PSE calculations. As the instability modes reach fi
nite amplitudes, the distinction between the linear and
nonlinear formulations also becomes visible. In the ab
sence of nonlinearity, the fundamental mode continues
to grow exponentially downstream, whereas nonlinear
interactions should cause the mode to saturate at ener
gies near E,,F 3 x 10 3. Some of this energy is
transferred to higher harmonics of the system, such as
the 920 mode, which grows several orders of magnitude
as the fundamental mode reaches saturation.
Similar to the interfacial case, nonlinear effects also
play a major role in the streamwise development of pre
viously unimportant instability modes. Although the
first harmonic ~o2 is weakly unstable, it quickly be
comes neutrally stable as the mean flow spreads, be
fore being excited by the fundamental mode as men
tioned above. For the 930 and 440 modes, which are ini
tially stable at the inlet, the nonlinear interactions with
other unstable modes result in their immediate excita
tion and eventual destabilization (see figure 9b). While
linear theory predicts that these modes simply decay
and are negligible downstream, the results of the non
linear theory suggest that they serve as an energy sink
and could reach amplitudes comparable to the lower fre
quency modes. Therefore, the nonlinear amplification of
the 920 or 940 mode downstream of the inlet region may
be responsible for the creation of smaller scale ripples
on top of an interfacial wave.
Three Dimensional Mixing Layers
In this section we examine the development of three
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
dimensional instability waves using the second two
phase mixing layer. The parameters of this mixing layer
are given by pr = 2, p = 1, We = 300, using the
same ReAu as the previous mixing layer. The discrete
forcing frequencies for the liquid modes encompass both
low frequency, linear unstable modes and higher fre
quency, linearly stable modes. In this particular mix
ing layer, the liquid modes are the dominant instability
modes which perturb the interface, and also includes a
sufficiently broad spectrum of unstable spanwise har
monics which can interact with other two and three
dimensional modes.
The inclusion of these spanwise harmonics pro
vides an opportunity to examine the growth of three
dimensional structures in twophase mixing layers. In
previous DNS studies of singlephase, compressible and
incompressible mixing layers (Rogers and Moser 1992),
oblique modes were observed to play a role in the de
velopment of streamwise rib vortices and the spanwise
'kinking' of vortex rollers. Although the current mix
ing layers under consideration do not display the same
vortical structures, threedimensional effects have been
shown to be important to the formation of liquid fingers
and droplets in the two fluid case (Lozano et al. 2001).
1 5 
05
1 i
1 5
200 250 300 350
Figure 10: Instantaneous interfacial wave height f(x, t)
at z 0 in the pr 2, f 1 mixing layer, with pri
mary forcing at the 911 liquid mode. The different line
styles correspond to nonlinear ( ) and linear ()
computations.
In using nonlinear PSE to simulate the pr 2,
p 1 mixing layer with liquid modes, two spanwise
forcing schemes were employed: one in which the pri
mary forcing mode was 911 (with spanwise wavenumber
S= 0.10), and the other with 912 (and 3 0.20) as the
primary forcing mode. Similar to the findings of Pierre
humbert and Widnall (1982), we observe that the initial
growth rates of the instability waves were relatively in
sensitive to the spanwise wavenumber at low 3, but the
oblique modes are responsible for the development of
different spanwise structures downstream.
The comparison of the interfacial wave heights com
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 100 200 300
x
0 100 200 300
0 100 200 300
x
0 100 200 300
0 100 200 300
x
0 100 200 300
Figure 11: Modal energy E,, for liquid modes in the p, 2, 1 mixing layer, with (a)(c): primary forcing
from the 911 liquid mode, and (d)(f): primary forcing from the 912 mode. The different line styles correspond to
nonlinear ( ) and linear () calculations, with the thick line ( ) representing the primary forcing mode.
puted by the nonlinear and linear simulations is consis
tent with the results of the twodimensional mixing lay
ers in the previous section. Because the linear predic
tions do not properly capture the saturation mechanism
for instability waves, they generally tend to overestimate
the growth of the interface compared to the equivalent
nonlinear computations (figure 10). From the results of
the threedimensional mixing layers, we find that sub
stantial energy transfer occurs between not only between
modes with different harmonics, but also between modes
with different spanwise wavenumbers. In figures 1 lb
and 1If, we see that the primary energy transfer route
occurs in the excitation of an immediate harmonic by
the primary forcing mode (i.e, between 911 and 921, and
between 912 and 922). However, in all simulations, the
growth of the twodimensional modes were still affected
by the presence of the oblique modes. In figure 1 la and
1 Id, we find that the growth of the Plo mode is notice
ably increased throughout the streamwise domain, while
the >20 mode gains additional energy near the satura
tion point of the primary forcing mode. Additionally,
in figure 11, it should be noted that all linear predic
tions are obtained using the linear PSE. More prominent
differences between the two forcing schemes can be ob
served in the threedimensional interfacial wave patterns
that develop. As illustrated in figure 12, forcing due to
the 012 mode leads to greater spanwise variation, as ex
pected. This forcing method produced larger height dif
ferences between adjacent wave crests in the spanwise
direction, and enforces the idea that higher frequency
oblique modes are important to the development of liq
uid fingers. An examination of the threedimensional
patterns produced by the linear simulations yielded un
physically large amplitude variations in the spanwise z
direction.
Discussion
From the results of the previous section, we have seen
that nonlinear effects can both limit the growth of the
instability modes, but also excite previously stable high
frequency modes. In the following discussion, we use
the capabilities of this formulation to examine individual
nonlinear effects, including the influence of the mean
flow correction, and the effect of the initial amplitude on
the saturation mechanism.
In the context of singlefluid mixing layers, previous
studies have demonstrated that appropriately account
ing for the mean flow correction is essential in cap
turing the streamwise development of instability waves
105
1010
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 12: Spanwise plot of instantaneous interfacial wave height f(x, z, t) for liquid modes in the pr
mixing layer. (a) Primary forcing from 11 mode. (b) Primary forcing from 912 mode.
2, p = 1
150 200 250 300
Figure 13: Modal kinetic energy Em for liquid modes
in the pr = 2, r = 1 mixing layer. The different line
styles correspond to nonlinear calculations (), lin
ear calculations (), and nonlinear calculations with
no mean flow correction ().
and the vortical dynamics in shear flows (Cheung and
Lele 2009). In the context of the twofluid mixing layer,
we can isolate the effect of the corrected mean flow by
comparing simulations both with and without the mean
flow forcing term FX,BL in equation (4b). For the two
dimensional pr 2, pr 1 mixing layer, the results
can be seen in figure 13. In the situation where only
modetomode nonlinear interactions were allowed, the
eventual saturation of the fundamental mode occurred
at higher amplitudes, and the first harmonic mode ab
sorbed larger amounts of energy, when compared against
the fully coupled solution with the mean flow correction.
Correspondingly, the absence of the mean flow correc
tion leads to larger wave amplitudes for the interface.
x 103
2
4
6 We
8
100 150 200 250 300
Figure 14: The mean interface location f(x) for liquid
modes in the pr = 2, = 1 mixing layer, with increas
ing curves corresponding to We=200, 300, 600, 1000.
In addition, we find that the inclusion of the mean
flow correction is also responsible for deformations to
the mean interface location f(x). Typical linear theo
ries, using a laminar mean flow for the mixing layer,
generally specify a zero vertical velocity at the center
line, causing the mean interface to remain at y = 0.
However, in the nonparallel, spatially evolving mixing
layer, the upper and lower streams experience different
convective decelerations, and nonlinear interactions can
introduce a nonzero mean vertical velocity at the center
line, which causes the mean interface to deform (figure
14). As the Weber number increases, modifications to
f(x) appear earlier upstream, and even at higher lev
els of surface tension, the mean interface location is still
displaced. This result is qualitatively consistent with the
earlier work of Renardy and Renardy (1993), who also
observed a longwave modulation of the interface height
in twolayer CouettePoiseuille flow. Their formulation,
however, did not incorporate a mean flow correction
and also enforced interfacial boundary conditions at the
nominal interface location.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Conclusions
In this study, we have extended the traditional parallel
flow, linear stability analysis of twofluid mixing layers
to regimes where spatial evolution and nonlinear effects
are nonnegligible. To accomplish this, we developed
a formulation based on the Parabolized Stability Equa
tions, which allows us to investigate nonlinear mode and
mean flow interactions, the mean flow correction, and
finite amplitude effects on the development of the two
fluid interface.
Our results have indicated that the spatially evolution
of the mixing layer can significantly alter the outcome
of mode competition between liquid modes and interfa
cial modes. Some interfacial modes, which are initially
less unstable than liquid modes, can substantially am
plify due to the spreading of the mean flow grow and
overtake liquid modes at the same frequency. This im
plies that local LST predictions of the dominant mode
can fail to provide an accurate picture of the develop
ment of instability modes downstream.
Several nonlinear effects were also highlighted
through comparisons between linear and nonlinear cal
culations of spreading mean flows. In order to obtain the
correct amplitude growth for instability waves, both the
mean flow correction and the the nonlinear modal inter
actions must be accounted for. We also find that energy
transfer from the fundamental mode to higher harmon
ics preserves or amplifies smaller scale features in the
interfacial wave. Additionally, energy is also transferred
between the instability waves and the mean flow, alter
ing the mean interface location.
Acknowledgements
This research was supported by a Marie Curie Inter
national Incoming Fellowship within the 7th European
Community Framework Programme, as well as the Im
perial Marshall Sherfield Postdoctoral Fellowship.
Appendix
For the linear PSE operator given in equation (16) acting
on the vector ,,n = [,, ,,mm,, m, T, the matri
ces A, B, C, D, E, G, M,and N are given by
1000
1 0 0 0
C 0 1 0 0
Re 0 0 1 0
0 0 0 0
auE
P a
0
N 
Re
pU M
0
0
1
,a u
, + E,
.aw
0
0
pU M
0
0
0 + El
0
0
0
pU 
0
0
0
0
where E, = (a2 + 32) and M = 2., ,.
References
A. Alexakis, A. C. Calder, L. J. Dursi, R. Rosner, J. W.
Truran, B. Fryxell, M. Zingale, F. X. Timmes, K. O1
son, and P. Ricker. On the nonlinear evolution of wind
driven gravity waves. Physics of Fluids, 16(9):3256
3268, 2004.
F P Bertolotti, Th. Herbert, and P R. Spalart. Linear and
nonlinear stability of the Blasius boundary layer. Jour
nal ofFluid Mechanics, 242:441474, 1992.
Thomas Boeck, Jie Li, Enrique LopezPages, Philip
Yecko, and Stephane Zaleski. Ligament formation in
sheared liquidgas layers. Theoretical and Computa
tional Fluid Dynamics, 21(1):5976, 2007.
L. C. Cheung and S. K. Lele. Linear and nonlinear pro
cesses in twodimensional mixing layer dynamics and
sound radiation. Journal ofFluid Mechanics, 625:321
351, 2009.
M. J. Day, N. N. Mansour, and W. C. Reynolds. Nonlin
ear stability and structure of compressible reacting mix
ing layers. Journal of Fluid Mechanics, 446:375408,
2001.
Hang Ding, Peter D.M. Spelt, and Chang Shu. Diffuse
interface model for incompressible twophase flows with
large density ratios. Journal of Computational Physics,
226(2):20782095, 2007. ISSN 00219991.
M. Fulgosi, D. Lakehal, S. Banerjee, and V DeAngelis.
Direct numerical simulation of turbulence in a sheared
airwater flow with a deformable interface. Journal of
Fluid Mechanics, 482:319345, 2003.
Mikhael Gorokhovski and Marcus Herrmann. Modeling
primary atomization. Annual Review ofFluid Mechan
ics, 40(1):343366, 2008.
Denis Gueyffier, Jie Li, Ali Nadim, Ruben Scardovelli,
and St6phane Zaleski. Volumeoffluid interface track
ing with smoothed surface stress methods for three
dimensional flows. J. Comput. Phys, 152:423456,
1999.
T. Herbert. Parabolized stability equations. Annual Re
views ofFluidMechanics, 29:245283, 1997.
M. Herrmann. A balanced force refined level set grid
method for twophase flows on unstructured flow solver
grids. Journal of Computational Physics, 227(4):2674
2706, 2008. ISSN 00219991.
E. J. Hinch. A note on the mechanism of the instability
at the interface between two shearing fluids. Journal of
FluidMechanics, 144:463465, 1984.
A. P. Hooper and W. G. C. Boyd. Shearflow instability
at the interface between two viscous fluids. Journal of
FluidMechanics, 128:507528, 1983.
A. P. Hooper and R. Grimshaw. Nonlinear instability
at the interface between two viscous fluids. Physics of
Fluids, 28(1):3745, 1985.
Michael R. King and Mark J. McCready. Weakly non
linear simulation of planar stratified flows. Physics of
Fluids, 12(1):92102,2000.
Jie Li and Yuriko Renardy. Numerical study of flows of
two immiscible liquids at low reynolds number. SIAM
Review, 42(3):417439, 2000.
Jie Li, Yuriko Y Renardy, and Michael Renardy. A nu
merical study of periodic disturbances on twolayer cou
ette flow. Physics ofFluids, 10(12):30563071, 1998.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Paolo Lombardi, Valerio De Angelis, and Sanjoy Baner
jee. Direct numerical simulation of nearinterface turbu
lence in coupled gasliquid flow. Physics ofFluids, 8(6):
16431665, 1996.
Antonio Lozano, Felix Barreras, Guillermo Hauke, and
Cesar Dopazo. Longitudinal instabilities in an air
blasted liquid sheet. Journal of Fluid Mechanics, 437:
143173,2001.
M. R. Malik and C. L. Chang. Nonparallel and nonlinear
stability of supersonic jet flow. Computers and Fluids,
29(3):32765,2000.
J. W. Miles. The hydrodynamic stability of a thin film
of liquid in uniform shearing motion. Journal of Fluid
Mechanics, 8:593610, 1960.
John W. Miles. On the generation of surface waves by
shear flows. Journal ofFluid Mechanics, 3(02):185
204, 1957.
R. T. Pierrehumbert and S. E. Widnall. The two and
threedimensional instabilities of a spatially periodic
shear layer. Journal of Fluid Mechanics, 114:5982,
1982.
Michael Renardy and Yuriko Renardy. Derivation of am
plitude equations and analysis of sideband instabilities in
twolayer flows. Physics ofFluidsA: Fluid Dynamics, 5
(11):27382762, 1993.
Yuriko Renardy. Weakly nonlinear behavior of peri
odic disturbances in twolayer CouettePoiseuille flow.
Physics ofFluidsA: Fluid Dynamics, 1(10): 16661676,
1989.
Yuriko Renardy and Daniel D. Joseph. Couette flow
of two fluids between concentric cylinders. Journal of
Fluid Mechanics, 150:381394, 1985a.
Yuriko Renardy and Daniel D. Joseph. Oscillatory in
stability in a B6nard problem of two fluids. Physics of
Fluids, 28(3):788793, 1985b.
M. M. Rogers and R. D. Moser. The threedimensional
evolution of a plane mixing layer: the KelvinHelmholtz
rollup. Journal ofFluidMechanics, 243:183226, 1992.
M. Sangalli, C. T. Gallagher, D. T. Leighton, H.C.
Chang, and M. J. McCready. Finiteamplitude waves
at the interface between fluids with different viscosity:
Theory and experiments. Phys. Rev. Lett., 75(1):7780,
Jul 1995.
Massimo Sangalli, Mark J. McCready, and HsuehChia
Chang. Stabilization mechanisms of short waves in strat
ified gasliquid flow. Physics of Fluids, 9(4):919939,
1997.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
M. J. South and A. P. Hooper. Linear growth in twofluid
plane Poiseuille flow. Journal ofFluidMechanics, 381:
121139, 1999.
Mark Sussman, Peter Smereka, and Stanley Osher. A
level set approach for computing solutions to incom
pressible twophase flow. J Comput. Phys., 114(1):146
159, 1994. ISSN 00219991.
T. L. van Noorden, P. A. M. Boomkamp, M. C. Knaap,
and T. M. M. Verheggen. Transient growth in parallel
twophase flow: Analogies and differences with single
phase flow. Physics ofFluids, 10(8):20992101, 1998.
F. M. White. Viscous FluidFlow. McGrawHill, 1991.
P. Yecko and S. Zaleski. Transient growth in twophase
mixing layers. Journal ofFluid Mechanics, 528:4352,
2005.
Philip Yecko, St6phane Zaleski, and JoseMaria Fullana.
Viscous modes in twophase mixing layers. Physics of
Fluids, 14(12):41154122, 2002.
C. S. Yih. Instability due to viscosity stratification. Jour
nal ofFluidMechanics, 27:337352, 1967.
