Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Linear and Nonlinear Analysis of the KelvinHelmholtz Instability
with the One Dimensional Two Fluid Model
William Fullmert, Deoras Prabhudharwadkar, Avinash Vaidheeswaran,
Victor Ransom and Martin Lopez de Bertodanot
School of Nuclear Engineering, Purdue University, 400 Central Dr.,
West Lafayette, IN 47907
twfullmer@purdue.edu, *bertodan@purdue.edu
Keywords: twofluid model, illposed, KelvinHelmholtz, stratified wavy flow
Abstract
The Inviscid KelvinHelmholtz (IKH) instability in horizontal stratified flow has been the subject of considerable research
within the framework of the one dimensional Two Fluid Model (1D TFM). It is of particular interest to the analysis of the
stability of the model because it initiates the high frequency waves that are amplified as a result of the illposed nature of the
model.
The objective of the present work is to analyze the basic 1D TFM and how it is affected by inclusion of transverse gravity (i.e.
hydrostatic force) and surface tension. The basic case is illposed while the latter is well posed. Linear stability analyses (i.e.,
dispersion analysis and Von Neuman analysis) and nonlinear 1D and 2D CFD TFM simulations are performed. It is observed
that linear theory describes the initial stage of wave formation but that nonlinear effects become dominant later on. The 1D
results are explained by the nonlinear analysis of Kreiss and Yalom (2008) that showed that the behavior of a system of
Burgers equations with a small value of artificial viscosity produces dissipative shocks that limit the amplitude of the waves.
The results are compared with the experimental data of Thorpe (1969). The data show that the waves reach an asymptotic
amplitude. The 2D CFD simulations show details of the breaking waves and the vortices that develop. The 1D TFM with
surface tension cannot simulate breaking waves or vortices but shocks are obtained instead which combine with the first order
upwind difference scheme to limit the amplitude of the waves by numeric viscous dissipation.
The inclusion of surface tension results in wavelengths that are similar to the experimental values, but the amplitude of the
waves is too large. Additional constitutive relations to account for the breaking of waves should improve the accuracy of the
model.
Introduction
The Inviscid KelvinHelmholtz (IKH) instability in
horizontal stratified flow has been the subject of
considerable research within the framework of the one
dimensional Two Fluid Model (1D TFM). It is of particular
interest to the analysis of the stability of the model because
it is the mechanism that produces high frequency waves
that are amplified as a result of the illposed nature of the
model.
Taitel and Dukler (1976) used linear stability theory to
predict the flow regime transition from stratified to slug
flow when the IKH instability occurs. In this simple
interpretation the onedimensional TFM becomes illposed
at the transition. However Ramshaw and Trapp (1976)
found that adding surface tension makes the model
wellposed but the fastest growing waves obtained from
linear stability theory, i.e., capillary waves, are very small.
Holmas et. al. (2008) added artificial viscosity to the
onedimensional TFM mass and momentum equations to
obtain a wellposed model that had a prescribed fastest
growing wavelength. This is equivalent to the standard
industrial practice to use the upwind scheme with
computational cells of arbitrary size to obtain significant
numerical viscosity. However, Holmas et. al. (2008) did
not offer a physical basis for the artificial diffusion in the
mass conservation equations which is inconsistent with the
TFM formulation and they did not compare their results
with experimental data.
It has been shown by Kreiss and Ystrom (2006) that the
nonlinear behavior of a system of Burgers equations
similar to the 1D twofluid model with a small value of
artificial viscosity results in dissipative shocks that limit
the amplitude of the waves. This theoretical result is
significant because it shows that the amplitude of the
waves in the TFM with the upwind scheme is limited even
when the mesh size is small. The shock steepening and
formation is analogous to the wavelength cascading results
reported by Krishnamurthy and Ransom (1992). They
analyzed a toroidal stratified flow where a long wavelength
disturbance was transferred to shorter wavelengths, until
the energy was dissipated by numerical viscosity. They
hypothesized that this process was similar to turbulent
eddies cascading down to the Kolomogrov length scale.
Bartosiewicz et. al. (2008) compared a 2D TFM CFD
model with the IKH data of Thorpe (1969). In this case the
waves actually break and the amplitude reaches an
asymptotic value. The breaking waves in the 2D model
occur instead of the shocks that appear in the 1D model.
The objective of the present work is to analyze the 1D
TFM including gravity (i.e. hydrostatic force) both with
and without surface tension. The first case is wellposed
while the other is illposed. Linear stability analyses (i.e.,
dispersion analysis and Von Neumann analysis) and
Paper No
nonlinear simulations with the TFIT 1D computer model
of Ransom (2000) and with 2D CFD simulations using
FLUENT are performed. The results are compared with the
data of Thorpe (1969).
Nomenclature
interfacial drag per unit volume (Nm 3)
gravitational constant (ms 2)
channel height (m)
liquid level (m)
imaginary unit
wavenumber (m1 )
pressure (Nm 2)
velocity (ms 1)
Greek letters
a void fraction ()
A finite difference
E artificial continuity viscosity
S independent vector
X wavelength (m)
g viscosity (Pas)
p density (kgm3)
o surface tension (Nm1)
c frequency (s')
5 eigenvalue
Subsripts
f heavy fluid
g light fluid
I imaginary
i interface
j control volume edge
m control volume center
x axial component
y transverse component
Supersripts
n time step
T transpose
t complex conjugate
Dispersion Analysis
In this paper the onedimensional twofluid model will be
used to analyze stratified channel flow of two immiscible
fluids. The model used in the subsequent analysis, Eq. (1)
 (4), is equivalent to the area averaged twofluid model of
Ishii and Hibiki (2006) with the following simplifications:
mass transfer is generally neglected; wall shear, axial
(viscous) shear and turbulent stresses are assumed to be
insignificant; all covariance terms are ignored. An energy
or temperature equation need not be modeled as all
conditions considered will be isothermal. The equations of
mass and momentum are
ap, +apv = 0 (1)
Dt x
S(1a) p +(1a)pvf =0 (2)
t ax
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
V BDv
apg + apv =
at apv ax (3)
P()
a apgg,+Mg
Dx
av av,
(1a)pf D +(1a)pfvf x=
(4)
(1 a) f (1a)pf g,+Mf
ax
The standard subscripts g and f will be used here where g
is identifies the lighter fluid and f identifies the heavier
fluid. The momentum sources, Mik contain the generalized
drag terms of steady drag, virtual mass, Basset, lift and
turbulent diffusion forces. For this analysis we will only
consider the steady drag force using a drag coefficient
model,
g =CDPf vg v f (5)
for the lighter phase. The equal and opposite drag force is
applied to the heavier phase. The average phasic pressures
can be given in terms of the interfacial pressure and a
hydrostatic pressure. Additionally the interface pressure
jump can be related by means of the surface tension force
a2h
Pig P = x (6)
for a liquid level h(x)= (l a)H where H is the
channel height. The two average phase pressures can then
be written as a function of the gas side interfacial pressure
and the void fraction.
Pg = Pig TPgyaH
P = Pig +i gyH (1 )+aH x(8
The second terms on the RHS of equations (7) and (8) are
the hydrostatic forces for a rectangular channel.
To simplify the analysis both phases will be considered
incompressible. This is an acceptable approximation
because it can be shown that for the conditions considered
here, the acoustic roots always exist and are not
responsible for the model being illposed. Therefore
expanding the derivatives in Eq. (1) (4), neglecting
density gradients and introducing the hydrostatic and
surface tension forces given by Eq. (7) and (8), the original
system of partial differential equations can be cast into
vector equation form
Dq a D 2 9 3 9
A +B + E +H +C=0 (9)
Dt Dx Dx Dx3
where 9== [a,VgVf, Pig is the vector of the
independent variables. The coefficient matrices and the
source vector C are given in the appendix. The second
order term is included in Eq. (9) although no physical
second order terms appear in the governing equations. This
artificial diffusion term is added for comparison of the
model of Holmas et al (2008) to the von Neumann analysis
of the discrete equations which contains inherent
Paper No
numerical viscosity. This term is not meant to represent
physical phenomena as it contains both velocity diffusion
in the momentum equations and void diffusion in the
continuity equations. The third order coefficient matrix
contains only one entry for surface tension force
contribution to the liquid momentum equation.
Equation (9) is a quasilinear partial differential equation
that renders itself to a dispersion relationship of Fourier
component representation (Ramshaw & Trapp 1975). To
perform the dispersion analysis an infinitesimal
perturbation, (p is introduced about an initial reference
state. The reference state, which is satisfied automatically,
is subtracted out and the equation is linearized by
neglecting products of disturbances. Assuming either a
steady state initial condition or that the perturbation
wavelength is much smaller than the characteristic
lengthscale of the reference state (Pokharna et al. 1997),
the resulting perturbation equation is then
A +B +E
at Bx Bx2
a3 a(10)
+H 3 + p= 0
ax aqp
The perturbation is assumed to have the form of a traveling
wave
S' = SPoei(kox) (11)
where (p90 is the initial amplitude, k is the wave number
and co= wR + io is the frequency. The speed of propagation
is given by fR/k and the growth (or decay) rate is given by
ao. Substituting Eq. (10) into Eq. (11) yields a
homogeneous system of a equations,
i(A+ikB +(ik)2 E+(ik)3 H+ &=0 (1
which must satisfy
det A kB ik2E +k3H +i =0 (1:
I 9p
for a nontrivial solution to exist.
The roots of Eq. (13) were determined using Mathematica
with the following physical parameters: pf = 1000 kg/m3,
pg = 780 kg/m3, f = 0.1 m/s, v = 0.1 0.5 m/s, a= 0.5,
a= 0.04 N/m, H = 0.03 m and the flow to be horizontal, ie
gx = 0. The properties have been chosen based on the
experimental setup of Thorpe (1969) which will be used to
benchmark some of the numerical results. The use of a
relatively dense secondary fluid, rather than a gas, makes it
possible for to destabilize the flow at low relative
velocities. It also has the effect of decreasing the critical
wavelength and increasing the growth rate.
Figure 1 shows the growth rate obtained by solving Eq. 13
for the imaginary component of the growth rate, Im((o) =
ro, for the wavelengths of practical importance, A = 2r/k.
All curves are for vf = 0.1 and vg = 0.5 m/s for which the
KH stability condition is not satisfied. The basic model
corresponds to the simplest twofluid model which
neglects all second order terms, surface tension,
hydrostatic pressure forces and interfacial shear. Clearly
this model is illposed (in the strictest sense) because oA 
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
co as A 0. Interfacial drag is tested as a (mathematical)
method of stabilization, as shown in Figure 1 in red for CD
= 0.1. However, this changes the growth rates very little,
the long wavelength growth rates are damped slightly but
the result approaches the basic model as the wavelength
approaches zero. This is due to the algebraic modeling of
the drag, meaning it will not affect the characteristic
eigenvalues (Pokharna et al. 1997).
Im(w)[l/s]
 Basic
 Drag
Hydrostatic
 Surface Tension
 Diffusion
I , .[m]
0.02 0.04 0.06 0.08 0.10
Figure 1: Growth factors from dispersion analyses.
When only the hydrostatic pressure contribution is
considered, i.e. the first order void gradient terms, a slight
stabilizing effect is seen over the entire domain even at
very small wavelengths. However, the model remains
illposed. In dispersed flow the coefficient of the void
gradient is modeled by an algebraic expression dependant
S on the square of the relative velocity. Park et al, (1998)
have shown that for monodispursed bubbly flow, the range
of stability can be increased by increasing the coefficient
of this term from 14 to unity. In fact the CATHARE
twofluid code includes this term where the coefficient
includes an arbitrary term chosen to make the model
wellposed (Bestion 1990). In the present stratified case
the KH stability condition is not satisfied so the
hydrostatic term is not enough to make the model
wellposed.
When the surface tension is included the model becomes
well posed. Inclusion of surface tension results in a cutoff
wavelength (around 2 cm) where the growth vanishes.
Also seen is a critical wavelength where the growth rate is
maximum. This is a well known result (Ramshaw & Trapp
1978) but it has been previously criticized because, for
physical values of surface tension, numerical
discretizations on the order of millimeters are required in
order to observe the cutoff. This in turn causes further
concern from some who suggest that the 1D model (which
has been averaged over an area presumably much larger
than the millimeter scale) should not be used for short
wavelength phenomena. Although the following numerical
analysis is indeed "pushing the boundary" of the validity
of the 1D twofluid model, insightful conclusions can still
be drawn from the results.
The final curve in Figure 1 corresponds to the basic model
plus the diffusion matrix E with diffusion in both mass and
Paper No
momentum equations of both phases. The viscosities of the
continuity equations are assumed equal magnitude to the
momentum viscosities of 0.001 Pa s. If only the realistic
diffusion of the momentum equations is considered, it does
not result in a cutoff wavelength. The added viscosity to
the entire system produces a result similar to representing
the differential equations as finite difference equations
with numerical viscosity. As will be shown in the
following section, even illposed differential models do not
have infinite growth as A 0 when they are represented
numerically. However this should not be understood as
rendering the model wellposed because the cutoff
wavelength is dependant on the mesh (grid spacing) and
will result in a decreasing critical wavelength as the mesh
is refined.
Von Neumann Analysis
The von Neuamnn analysis of the finite difference model is
equivalent to the dispersion analysis of the differential
model where the growth rate as a function of wave number
can be determined. The mass and continuity equations,
including hydrostatic and surface tension effects, are
discretized in Eqs. (14) (17) for a staggered grid where
the flow properties of void fraction, pressure and density
are stored at the volume centers, m, and the phase
velocities are stored at the cell boundaries, j = m /2. The
numerical method also uses first order upwinding where
Eqs. (14) (17) have been evaluated for cocurrent positive
flow. The momentum convection is first order upwind and
the third order surface tension term has been evaluated
with a forward/centered scheme (also called twothirds
upwind) at the j 2 position. The overbars indicate the
average properties of the j 2 positions.
n nl n+ In
+Vf jf (15)
 a l ) pn At (16) 0
Wp 1 ap A
\ s l m '
+ w v v v =
At
_a  (16)fn
Ax
() gn n+1 tv v
^^Ai Ax
)L H ,1 a .
p ,g fig, [v f ] At
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
n +1 n
(l)Pf) vf
+ I _n n[ n At
i i Ax
+((la)p) vf [v v. =
(l _ n+l ]
x (17)
+ gyH (1 a )p, (a, L 1 a
+(1ap gA1 t+f,[ v 1 At
 I At
aH(1a)[a n13a ,+3, 1a,] 2A
] Ax
As with the dispersion analysis, both phases are
approximated as incompressible, eliminating all density
difference terms. Then the continuity equations are
centered about a mass cell center, m, and the momentum
equations centered about a cell boundary, j. The products
of independent variables are again linearized, for example:
n+l n+l
,Vg ~j+l9: 1
n+l a n+l
a a v > (18)
m +1/2 m1 g l/2
1(Vg,+11/2 g,m 1/2 )n g (a a l)
The von Neumann analysis is preformed on the linearized
equations by considering the evolution of a single Fourier
mode, from old time to new time, by decomposing the
independent variables into p(x,t) = q(t)e After
some rearrangement, a common e il and e can be
canceled from the mass and momentum equations,
respectively. Moving the new time variables to the left
hand side and the old time variables to the right results in a
^n+l n
linear system of the form Np = 9p Multiplying
both sides by the inverse of the new time matrix finally
leads to the recursion relation for the new time variables in
terms of a so called growth factor or amplification matrix,
G.
^n+l ^n 1 ^n
9p = G =N 19p (19)
The stability requirement of a single equation is Igl < 1 for
all wavelengths to have a stable, bounded solution [5].
This is extended to the multiequation system requiring
that the modulus of the eigenvalues of the amplification
matrix, f, be less than or equal to unity, f e < 1.
The eigenvalues were found with Mathematica using the
conditions listed in the previous section with Ax = 1mm
and At = 0.5 ms. The eigenvalues are adjusted to for
comparison with the growth factors of the dispersion
analysis by the relation V = e The adjusted
expression for the basic model and for the model with
hydrostatic pressure and surface tension is shown in Figure
2 for wavelengths from 2Ax to 10 cm. The growth at
wavelengths less than 2Ax has been neglected because it
can not be represented with finite spacing of Ax.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
corresponds to the most dangerous wavelength of the basic
model in Figure 2. It is also smaller than the cutoff
 Basic
 Hydrostatic + Surface Tension
0 l iI I '[m ]
0.00 0.02 0.04 0.06 0.08 0.10
Figure 2: Growth rates from von Neumann analyses.
The basic model in Figure 2 is similar to the differential
model with diffusion shown in Figure 1. The growth rate
increases exponentially as the wavelength decreases until
numerical damping becomes significant, then decreases
very rapidly to a minimum at 2Ax. A cutoff wavelength
exists near 5Ax, similar to the Holmas et al. model of the
dispersion analysis. On the other hand when surface
tension and hydrostatic pressure are included, the discrete
analysis is in excellent agreement with the differential
analysis. There is a cutoff wavelength with physical
significance around 2 cm and a critical wavelength,
although slightly skewed, around 3 cm. The major
conclusions of the linear studies will be compared to
nonlinear analysis by numerical calculations with the
scheme described by Eq. (14) (17).
Numerical Nonlinear Analysis
The TwoFluid Interfacial Temperature (TFIT) code is used
for the 1D analysis. TFIT was developed by Ransom
(2000) from the DISCON2 code (Trapp and Mortensen
1992, 1993). The finite difference scheme for the mass and
momentum equations is given in the previous section. The
code also solves energy equations for the phasic
temperatures, interface mass transfer and interface
temperature even though the problems presented here are
adiabatic.
A sample problem was set up to demonstrate numerical
aspects outlined by the preceding linear analysis. To
simplify the problem the wall and interfacial drag
coefficients were set to zero since the dispersion analysis
has shown that their effect is not very significant. The flow
and material properties are the same as given previously as
well as the mesh spacing and time step. A rectangular
channel with a length of 1 m has initial void fraction of 0.5
and phase velocities of 0.5 and 0.1 m/s. A long wavelength
perturbation is imposed on the void fraction at x = 0.3 m
with amplitude of 0.05 and wavelength 10 cm (well above
the critical wavelength). A short wavelength perturbation
of 1 cm with amplitude of 0.005 is added on the entire
domain. It is important to note that the 1 cm perturbation
,,! t,, i ,.... .... ... .;
i I Vi
Distance [m]
Figure 3: Evolution of long and short wavelength liquid
level perturbations for the basic model (left) and the model
with surface tension and hydrostatic pressure (right) at: (a)
initial condition, (b) t = 0.05s, (c) t = 0.10s and (d) t =
0.15s.
wavelength predicted by the linear analyses when surface
tension and hydrostatic pressure are considered. The
results for the first 0.15 s (300 iterations) are shown in
Figure 3.
The basic model on the left hand side of Figure 3 is
illposed, as expected. In the second frame it is clear that
the short wavelength perturbation is growing faster than
the long wavelength perturbation. When surface tension
and hydrostatic pressure are included however, the short
wavelengths decay and after 300 iterations the interface is
practically smooth. This is a highly desirable result. The
long wavelength perturbation on the other hand grows and
begins to develop a steep, shock like formation.
Waveforms of this type have been studied by Kreiss and
Ystrom (2006) for a system of Burger's equations
resembling the two fluid model. Their nonlinear analysis
demonstrates that the initially unstable process is bounded
when such shock type structures develop. Continued
evolution of this case would show that although more
shocktype wave forms develop from the exiting wave, the
amplitude of each is limited. This is due to the highly
dissipative interaction of the shocks and numerical
viscosity. In Figure 3(d) is it clear that the initial
wavelength of 10 cm has been reduced by the near
doubling of the wave. The same case was run on three
Paper No
Ln( t)/At [1/s]
120 r
*BMX' `Ti ffl'lJ``w I 'WM'J`tg `,iJ wm# ` 'XJJi'N ``,
i,,l,,,,,l
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
different meshes Ax = 2.0, 1.0 and 0.5 mm, with the same
courant number. The initial growth rate of the large wave
and the development of the shock are nearly identical for
each case. Figure 4 shows the wave at t = 0.2 s for the
three grids. The amplitude and the slope of the fully
developed shock are equal for each case. The development
of the adjacent waves is slightly different for each case as
the smaller spacing allows for more harmonic excitation.
However, as evident by the initial disturbance, each case
will result in waves with shocks of similar amplitude and
slope.
0.8
7 0.6
0 0.4
0.2
0.3 0.35 0.4 0.45 0.5 0.55
Distance [m]
Figure 4: Evolution of a 10 cm wave at 0.2 s.
To further validate the model, the growth rate for the first
200 iterations is compared with the VonNeumann analysis
shown in Figure 2. A set of single wavelength perturbation
cases were run that have initial conditions similar to Figure
3, without the short wavelength perturbation superimposed.
The wavelengths were varied from 2 to 10 cm and the
amplitude was measured for the first 400 iterations. In
every case the amplitude was saturated by the end of the
simulation. The amplitudes of the first 200 iterations were
(poorly) fitted with exponential curves to estimate the
growth rates. The fit of the curves are somewhat irrelevant
because of the nonlinear behavior. What is important is
the relative magnitude of the growth to determine the
cutoff and fastest growing wavelengths.
Figure 5 shows that the critical wavelength is closer to 4
cm than the 3 cm obtained by the linear analyses. This is
interpreted as a result of larger numerical diffusion present
in TFIT than was modeled by the von Neumann analysis.
The 2 cm wavelength is an approximation of the cutoff
wavelength; actually this wave initially decayed slightly,
but spread and by 100 iterations began to grow at a
wavelength of approximately 3 cm. The spreading or
diffusion of the wave is a property of the surface tension.
The 10 cm single wave perturbation was also used to find
the critical or destabilizing velocity of the lighter phase.
Both the dispersion and von Neumann analyses predict that
the critical velocity is 30 cm/s. The nonlinear result is
closer to 32 cm/s, but very similar to the linear
approximation.
000 0.01 0.02 0.03
4,*'4
0.04 0.05 0.06 0.07 0.08 0.09 0.10
Wavelength, X [m]
Figure 5: Approximate growth rates of single wave
disturbances after 200 iterations.
Experimental Benchmark
The results of the previous sections are now benchmarked
against experimental data. The experiment is that of
Thorpe (1969) who studied the KelvinHelmholtz
instability of two immiscible liquids; the heavier fluid
being water and the lighter a mixture of tetrachloride and
commercial paraffin with the properties given previously
and pf = 0.001 and pg = 0.0015 Pa s. The rectangular glass
channel was 1.83 m long with height 3 cm and width 10
cm, making the flow quasi2D. From a steady state initial
level of 1.5 cm (a= 0.5), the channel was suddenly tilted
an angle of 0 = 4.1 degrees causing a countercurrent flow,
settling type problem. The flow of each phase was
accelerated by the axial gravity source. By neglecting the
closed end effects and viscosity Thorpe (1969) estimated
that the velocity increases linearly during the first three
seconds as:
A_ phk sin(O)g
vk =+ A t, k=f g (20)
Phf +Pf hg
which can be used with the linear analyses to determine a
critical time for the onset of wave generation. In the
previous linear analysis velocity values were selected well
beyond the IKH criterion. When the relative velocity is
decreased to the point of neutral stability, the growth rate
decreases at all wavelengths, critical wavelength increases
and the previously distinguishable peak is broadened.
The critical relative velocity was determined for multiple
wavelengths as shown in Figure 6. The asymptotic critical
relative velocity is 19 cm/s which occurred at 1.1 s.
However the growth rates for these long wavelength
instabilities are very near unity. The growth rate becomes
significant at time 1.3 to 1.5 s. For the data shown below,
the instability was estimated to occur at 1.88 s. However
this includes approximately 0.25 s to compensate for half
of the time to tilt the tube.
Paper No
Paper No
035 
0.25
0.2
0s
0.1 
o1
0 0.05 0.1 0.15 0.2 0.25 0.3
Wavelength [m]
Figure 6: Critical velocity and time of IKH instability
The TFIT code was again used for the 1D simulations.
Again surface tension was included and wall drag was
neglected, interfacial drag was accounted for but has little
effect.
Fluent was used for the 2D simulations. The
EulerianEulerian twofluid model was chosen to compare
the results directly with the 1D model. The twofluid
model equations were solved using the phasecoupled
SIMPLE algorithm with a first order upwind scheme for
both momentum and volume fraction transport equations.
Uniform, rectangular meshed grids were used with a length
scale of 0.5 and 0.25 mm. The solution was found to be
mesh size independent and the reported results utilize the
finer mesh with a time step of 0.5 ms. For these
simulations interfacial drag is neglected, surface tension is
not included and no turbulence modeling is used, i.e. only
laminar stresses. All walls used the noslip boundary
condition. The initial condition was prescribed as perfectly
separated fluids with zero initial velocities. At the initial
time, the gravitational acceleration component was
modified to obtain the desired final tilted position of the
experiment, 4.1 from the normal. Surface tension can be
neglected for the 2D case because its principal effect is on
the shape of the waves. Bartosiewicz et al. (2008) has
shown that including surface tension in a 2D simulation
stiffens the interface which results in delayed wave
development, a slight increase in the critical wavelength a
reduction in wave amplitude and a reduction in the amount
of roll up wave breaking.
Discussion
The results of the 1D and 2D models are compared to the
data in Figure 8 for four characteristic wave development
stages. Note that these stages do not occur at the same time
for each case. Similar to the sample case in Figure 3, the
1D model exhibits wave growth to a point where sharp
void gradients similar to shocks appear on the windward
side of the waves. These highly dissipative regions are
responsible for controlling the growth in the 1D model (see
Figure 4). The shock steepening and formation has been
analyzed by Kreiss (2006). Clearly the 2D results which
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
show wave breaking and vortices are more physical. These
vortices were identified experimentally by Banner and
Phillips (1974) who observed their formation on the
forward face of the waves as the water tumbles forward
without necessarily having discontinuities in the slope.
Initially all of the vorticity is concentrated in a thin layer
surrounding the interface as the phases accelerate
countercurrently. As the KH instability develops and the
waves break the vortex sheet stretches and folds resulting
in a significant increase in viscous dissipation. Figure 7
shows contours of vorticity in the vicinity of a breaking
wave at 3 s. Here the vortex sheet is fully split with red and
blue signifying clockwise and counterclockwise rotation.
Figure 7: Vorticity surrounding a breaking wave in the 2D
simulation at 3.0s
In the 1D results there is a noticeable shift in the
predominant wavelength on the interface from the onset of
the disturbance to the fully developed waveforms. A Fast
Fourier Transformation (FFT) analysis of the interface
shows that shortly after the interface has become
completely wavy, the most predominant wavelength is near
4 cm. Over a period of about a second, the peak shifts to 6
cm. These measurements are similar to dominant
wavelengths of the 2D results: 1.8 3.2 cm and the
experimental data: 2.5 4.5 cm.
Another significant discrepancy between the 1D model and
the data is delay in the onset of the instability. The TFIT
interface only becomes wavy in the center of the domain
about a full second after the data and the 2D simulation.
This is in contrast to the calculations of the previous
section, where the TFIT results were in agreement with the
dispersion analysis to less than 1%. Although TFIT can
predict neutral stability of finite waves well, the inherently
stable numerical method has difficulty generating a finite
amplitude disturbance from an initially smooth interface.
The waves that appeared in the center of the domain were
originated by waves that were generated at the ends of the
domain. In reality there are short wavelength disturbances
throughout the experimental test section that are caused by
the environment, vibrations during tilting, etc.
To test this hypothesis TFIT was run with a 0.3 mm
amplitude, 10 cm long wavy perturbation. Figure 9 shows
a Fourier analysis of the interface in the center of the
domain. An FFT of the void fraction is integrated over all
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
t = 2.057 s
.meal r n,
t = 2.157 s
t = 1.8 s
t = 2.293 s
V' 1 "
t = 2.5 s
t = 2.411 s
'T .
Distance [m]
Figure 8: Waves resulting from tilted stratified flow predicted by the ID (TFIT) and 2D (Fluent) models compared to the data of Thorpe (1969 permission Cambridge University Press)
Paper No
t = 2.75 s
t = 1.5 s
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
wavenumbers at different times. In the
initial disturbance, the integral FFT sho\
after about one second which is due an
number of waves. Later the integral FFT
dramatically which corresponds to wave gr
0 05 1 15 2 25
Tim [s
Figure 9: Integral of FFT of the interfac
perturbation and sample wave forms at 2.0
This trend is extrapolated to approximate a
s. After two seconds the integral FFT
behavior. Figure 10 shows two typical wa
and 2.1 s. It is observed that when the shock
develop, the dissipation begins to reduce
the waves. At the same time the sur
dispersing the waves, which allows new v
in the center of the old structures. Th
increase in amplitude until dissipation halt
the cycle is repeated. This cycle will cc
initial perturbation wavelength, 10 cm,
final critical wavelength of approximately
1 1.02 1.04 1.06 1.08 1.
Distance [m]
Figure 10: Developed waves from an initial
presence of an
vs some growth
Conclusions
increase in the In this work we have used linear analyses of the
increases rather mathematical (dispersion) and numerical (von Neumann)
owth. differential equations for the two fluid model. The linear
analyses were used to test the frequency domain
characteristics of the model. The most important result is
that the inclusion of surface tension third order void
derivative term causes the growth factor to exhibit a cutoff
wavelength. This was verified by nonlinear analysis using
the TFIT code to simulate simple cases with initial void
fraction perturbations. The results show growth of the long
wavelength wave, while the small amplitude, short
wavelength disturbances below the cutoff frequency were
damped out.
The IKH experimental data of Thorpe (1969) was used as a
benchmark to test several aspects of the TFIT 1D model.
Additionally, twofluid 2D simulations of the experiment
were performed with the commercial CFD code Fluent.
The data show that the waves reach an asymptotic
amplitude and a predominant wavelength. This may be
3 35 4 explained by the wave braking process described by
Banner and Phillips (1974). Whereas the 2D two fluid
model reproduces the breaking phenomenon, the 1D model
e with an initial is not capable of simulating wave breaking, but
and 2. s. compensates by generating dissipative shocks if an upwind
scheme is used. The wavelength of both simulations were
in 'onset' at 1.88 similar to the data and to the linear analysis.
show oscillatory It was shown that the wave amplitude was overpredicted
ave forms at 2.0 with the 1D model. Furthermore it was shown that the
k type structures dissipation mechanism of the model was responsible for
the amplitude of stopping the wave growth. Physically the dissipation is
face tension is caused by wave breaking and the formation of vortices,
saves to develop which the 2D model is capable of simulating (i.e.,
ese new waves NavierStokes equations). However, the 1D model (i.e.,
s the growth and Burgers equations) relies on numerical dissipation at sharp
ntinue until the void gradient regions resembling shocks. These waveforms
has reached the were shown to be stable as predicted by the analysis of
6 cm. Kreiss and Ystrim (2006) as long as first order upwinding
is used. However, Liao et al. (2005) have shown that
second order upwinding schemes are much more unstable
that the simpler first order upwinding. So the advantage of
using higher order numerical schemes for stratified 1D
flow may be precluded if the two fluid model does not
have constitutive relations that can describe the effect of
the breaking waves and the vortices.
Acknowledgements
i This material is based upon work supported under a
Department of Energy Nuclear Energy University Program
S 1 Graduate Fellowship.
References
Banner M. L., Phillips, O. M., On the incipient breaking of
small scale waves. J. Fluid Mech., 65, part 4, 647656
(1974)
1 1.12 1.14 Bartosiewicz, Y., Lavieville, J. and Seynhaeve, J.M. A
first assessment of the NEPTUNE CFD code: Instabilities
al disturbance. in a stratified flow comparison between the VOF method
and a twofield approach, Int. J. Heat and Fluid Flow, 29,
Paper No
Paper No
460478 (2008)
Bestion, D. The physical closure laws in CATHARE Code.
Nucl. Eng. Design, 124, p229 (1990)
Holmas, H., Sira, T., Nordsveen, M., Langtangen, H. P.,
and Schulkes, R. Analysis of a 1D incompressible
twofluid model including artificial diffusion. IMA J of
Applied Math., 73, 651667 (2008)
Ishii, M. & Hibiki, T., Thermofluid dynamics of
twophase flow, Springer (2006)
Kreiss, K.O. and Ystrom, J. A note on viscous
conservation laws with complex characteristics. BITNumer.
Math., 46, S55S59 (2006) [Also see: Kreiss, K.O. and
Ystrom, J. Parabolic problems which are illposed in the
zero dissipation limit, Math. Comput. Modeling, 35,
12711295 (2002)]
Krishnamurthy, R. and Ransom, V. H. A nonlinear
stability study of the RELAP5/MOD3 twophase model.
Presented at JapanUS Seminar on TwoPhase Flow,
Berkeley, CA (1992)
Liao, J., Mei, R. and Klausner, J. F A study on numerical
instability of inviscid twofluid model near illposedness
condition. ASME Summer Heat Transfer Conference, San
Francisco, CA (2005)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Mortensen, G. A. and Trapp, J. A. A computational model
for twophase flow in onedimensional channels using a
discrete continuous formulation. EGGNRE1103 (1993)
Pokharna, H., Mori, M., and Ransom, V. H. Regularization
of twophase flow models: a comparison of numerical and
differential approaches. J. Comput. Phys., 134, 282295
(1997)
Ramshaw, J. D. and Trapp, J. A. Characteristics, stability
and shortwavelength phenomena in twophase flow
equation systems. Nucl. Sci. Eng., 66, 93 (1978)
Ransom, V. H. Summary of research on numerical methods
for twofluid modeling of twophase flow. Unpublished.
(2000)
Thorpe, J. A. Experiments on the instability of stratified
shear flow: immiscible fluids. J. Fluid Mech. 39 2548
(1969)
Trapp, J. A. and Mortensen, G. A. A discrete particle model
for bubbleslug twophase flows, J. Comput. Phys., 107,
367377 (1992).
vonNeumann, J. and Ritmyer, R. D. A Method for the
Numerical Calculation of Hydrodynamic Shocks, J.
Applied Phys., 21, 232237 (1949)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Appendix
Matrices of the dispersion analysis:
1 0 0 0
1 0 0 0
A=
0 ap, 0 0'
0 0 (1 a)p 0
, 0 0 0
E 0 0 0
0 g, 0 0
0 0 li 0
fjg (Vg v1)
gxfg H
gP sH
H
Matrices of the von Neumann analysis:
0
0
fpf, (Vg V
H
f Hp (Vg V
H
Vg
vf
B gyHapg
L g,H (1 a)p,
0 0 0
0 0 0
H=
0 0 0
a(1a)H 0 0
0
0
)
0
0
0 0
v ap, 0
0 vf (1 a)pf
0
0
0
0
.i2an kAx At
i2asin 
2 Ax
0
( l)inkAx At
i2(1 a) sin ]
1 2 Ax
igyHapg sin kA At
igH(n 1 2 Ax
igyH(1a)psin ] At
2 Ax
fAt + ap
fAt
fIAt
fAt+(1a)pf
i2asn kAx] At
i2asin I
12 Ax
i2(1a)sinkx] At 
S2 Ax
Sv (1 ek)At
Ax
vi (1 e" ) 1
Ax
0
0
v (1 e ikhA)A
0
Paper No
0
0
a
(1a)
f1v,(1eikA) (1a)pf 0
\ Ax)
*k x 3kAx At
ioH(1a) 6sin 2sin At
