Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 6.1.4 - Dispersion Analysis of Diffusion and Stiffness Filters for the One Dimensional Two Fluid Model
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00146
 Material Information
Title: 6.1.4 - Dispersion Analysis of Diffusion and Stiffness Filters for the One Dimensional Two Fluid Model Instabilities
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Fullmer, W.
Prabhudharwadkar, D.
Vaidheeswaran, A.
Ransom, V.
Lopez de Bertodano, M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: two-fluid model
ill-posed
Kelvin-Helmholtz
stratified wavy flow
 Notes
Abstract: The Inviscid Kelvin-Helmholtz (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 ill-posed 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 ill-posed while the latter is well posed. Linear stability analyses (i.e., dispersion analysis and Von Neuman analysis) and non-linear 1D and 2D CFD TFM simulations are performed. It is observed that linear theory describes the initial stage of wave formation but that non-linear effects become dominant later on. The 1D results are explained by the non-linear 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00146
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 614-Lopez-ICMF2010.pdf

Full Text

Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Linear and Non-linear Analysis of the Kelvin-Helmholtz 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: two-fluid model, ill-posed, Kelvin-Helmholtz, stratified wavy flow

Abstract

The Inviscid Kelvin-Helmholtz (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 ill-posed 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 ill-posed while the latter is well posed. Linear stability analyses (i.e.,
dispersion analysis and Von Neuman analysis) and non-linear 1D and 2D CFD TFM simulations are performed. It is observed
that linear theory describes the initial stage of wave formation but that non-linear effects become dominant later on. The 1D
results are explained by the non-linear 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 Kelvin-Helmholtz (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 ill-posed 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 one-dimensional TFM becomes ill-posed
at the transition. However Ramshaw and Trapp (1976)
found that adding surface tension makes the model
well-posed 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
one-dimensional TFM mass and momentum equations to
obtain a well-posed 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 two-fluid 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 well-posed
while the other is ill-posed. Linear stability analyses (i.e.,
dispersion analysis and Von Neumann analysis) and





Paper No


non-linear 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 (kgm-3)
o surface tension (Nm-1)
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 one-dimensional two-fluid 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 two-fluid 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(1-a) p +-(1-a)pvf =0 (2)
t ax


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

V BDv
apg + apv =
at apv ax (3)
P()
-a -apgg,+Mg
Dx
av av,
(1-a)pf D +(1-a)pfvf x=
(4)
-(1- a) f- (1-a)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 ill-posed. 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 quasi-linear 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(k-ox) (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 non-trivial 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 de-stabilize 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
K-H stability condition is not satisfied. The basic model
corresponds to the simplest two-fluid model which
neglects all second order terms, surface tension,
hydrostatic pressure forces and interfacial shear. Clearly
this model is ill-posed (in the strictest sense) because oA -


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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
ill-posed. 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
two-fluid code includes this term where the coefficient
includes an arbitrary term chosen to make the model
well-posed (Bestion 1990). In the present stratified case
the K-H stability condition is not satisfied so the
hydrostatic term is not enough to make the model
well-posed.
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 two-fluid 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 ill-posed 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 well-posed 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 two-thirds
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 fi-g, [v -f ] At


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

n +1 n
(l-)Pf) -vf

+ I _--n n[ n At
i i Ax
+((l-a)p) vf [v -v. =


-(l- -_ n+l ]
x (17)

+ gyH (1- a )p, (a, L 1- a
-+(1-ap gA1- t+f,[- v 1 At

-- I At


-aH(1-a)[a n1-3a ,-+3, 1-a,] 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 m-1 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 multi-equation 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 30-June 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 Non-linear Analysis

The Two-Fluid 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
ill-posed, 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 non-linear 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
shock-type 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 30-June 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 Von-Neumann 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 non-linear 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 Kelvin-Helmholtz
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 quasi-2D. 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 counter-current 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
Eulerian-Eulerian two-fluid model was chosen to compare
the results directly with the 1D model. The two-fluid
model equations were solved using the phase-coupled
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 no-slip 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 30-June 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
counter-currently. 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 30-June 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 30-June 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 non-linear 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, two-fluid 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 Navier-Stokes 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, 647-656
(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 two-field approach, Int. J. Heat and Fluid Flow, 29,


Paper No






Paper No


460-478 (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
two-fluid model including artificial diffusion. IMA J of
Applied Math., 73, 651-667 (2008)

Ishii, M. & Hibiki, T., Thermo-fluid dynamics of
two-phase flow, Springer (2006)

Kreiss, K.-O. and Ystrom, J. A note on viscous
conservation laws with complex characteristics. BITNumer.
Math., 46, S55-S59 (2006) [Also see: Kreiss, K.-O. and
Ystrom, J. Parabolic problems which are ill-posed in the
zero dissipation limit, Math. Comput. Modeling, 35,
1271-1295 (2002)]

Krishnamurthy, R. and Ransom, V. H. A non-linear
stability study of the RELAP5/MOD3 two-phase model.
Presented at Japan-US Seminar on Two-Phase Flow,
Berkeley, CA (1992)

Liao, J., Mei, R. and Klausner, J. F A study on numerical
instability of inviscid two-fluid model near ill-posedness
condition. ASME Summer Heat Transfer Conference, San
Francisco, CA (2005)


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Mortensen, G. A. and Trapp, J. A. A computational model
for two-phase flow in one-dimensional channels using a
discrete continuous formulation. EGG-NRE-1103 (1993)

Pokharna, H., Mori, M., and Ransom, V. H. Regularization
of two-phase flow models: a comparison of numerical and
differential approaches. J. Comput. Phys., 134, 282-295
(1997)

Ramshaw, J. D. and Trapp, J. A. Characteristics, stability
and short-wavelength phenomena in two-phase flow
equation systems. Nucl. Sci. Eng., 66, 93 (1978)

Ransom, V. H. Summary of research on numerical methods
for two-fluid modeling of two-phase flow. Unpublished.
(2000)

Thorpe, J. A. Experiments on the instability of stratified
shear flow: immiscible fluids. J. Fluid Mech. 39 25-48
(1969)

Trapp, J. A. and Mortensen, G. A. A discrete particle model
for bubble-slug two-phase flows, J. Comput. Phys., 107,
367-377 (1992).

vonNeumann, J. and Ritmyer, R. D. A Method for the
Numerical Calculation of Hydrodynamic Shocks, J.
Applied Phys., 21, 232-237 (1949)





7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 s-H
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(1-a)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(1-a)psin- -] At
2 Ax


fAt + ap


-fAt


-fIAt


fAt+(1-a)pf


i2asn kAx] At
i2asin --I
12 Ax

i2(1-a)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
(1-a)


f1-v,(1-e-ikA)-- (1-a)pf 0
\ Ax)


*k x 3kAx At
ioH(1-a) 6sin -2sin At




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs