THE NUMERICAL SIMULATION OF BUBBLE COALESCENCE
by
Valerie Ellen Walker
B.S., Colorado State University, 1984
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Mechanical Engineering
1991
I/U
This thesis for the Master of Science
degree by
Valerie Ellen Walker
has been approved for the
Department of
Mechanical Engineering
by
CrtXflA
J. Kenneth E. Ortega
Date
Walker, Valerie Ellen (M.S., Mechanical Engineering)
The Numerical Simulation of Bubble Coalescence
Thesis directed by Professor John A. Trapp
A computer model was developed to accurately simulate the
process of bubble coalescence in a bubbly flow regime. The
numerical model included proper drag relationships for the various
bubble shapes and regimes occurring in bubbly flow. A model of the
wake structure behind spherically capped bubbles was developed which
included models for both primary and secondary wake structures.
The numerical model was incorporated into an existing one
dimensional, mechanical, two phase flow program. The revised
program having the capability of tracking individual bubbles
throughout the flow area, with the ability to simulate bubbly flow
through the onset of slug flow.
Calculated values of bubble terminal velocities in water
i
compared well with measured data. Simulations of bubble coalescence
also compared well with measured bubble positions during coalescence
and coalescence times. The comparison were made for sugar water and
glycerin coalescences, yielding favorable results from close
distances (a few bubble radii) to separation distances up to 80 cm.
The form and content of the abstract are approved. I recommend its
publication.
Signed
CONTENTS
Figures......................................................vi
Acknowledgements .......................................... vii
Nomenclature ............................................. viii
CHAPTER
1. INTRODUCTION ............................................... 1
Background............................................... I
Scope..................................................... 3
2. PROGRAM STRUCTURE .........................................4
Governing Equations ....................................... 4
Program Scheme ............................................ 7
3. NEW CODE FEATURES......................................... 10
Bubble Hydrodynamics ..................................... 10
Drag Coefficients and Transition Criteria ................ 16
Spherical Regime ...................................... 16
Ellipsoid Regime ...................................... 20
Spherical Cap Regime....................................23
Bubble Geometry .......................................... 23
Bubble Coalescence ....................................... 25
Wake Structure........................................... 27
Primary Wake......................................... 27
Secondary Wake..........................................31
Drag on a Bubble in a Wake.................................35
Bubble Merging ........................................... 37
4. RESULTS AND DISCUSSIONS
39
Terminal Velocity ...................................... 39
Coalescence Times ........................................ 42
Original Model ....................................... 43
Final Model............................................45
Flow Visualization.........................................55
Discussion............................................... 58
REFERENCES....................................................... 61
APPENDIX ..........................................................62
A. Code Listing..................................................63
v
FIGURES
Figure
2.1 Program flow chart.........................................8
3.1. Bubble shapes.............................................11
3.2. Bubble regimes.............................................. 13
3.3. Reynolds number vs. drag coefficient......................17
3.4. Wake structure behind a spherical cap bubble (Bessler
and Littman) ............................................ 28
3.5. Primary wake shapes.......................................29
3.6. Primary wake velocity...................................... 32
3.7. Wake structure behind a spherical cap bubble
(code model).............................................36
4.1. Terminal velocity of bubbles in water
(Peebles and Garber 1953) 40
4.2. Terminal velocity of bubbles in water
(Haberman and Morton 1953) ................................ 41
4.3. Bubble positions during coalescence original model ... 44
4.4. Coalescence of same size bubbles original model .... 46
4.5.. Coalescence of a larger trailing bubble original model 47
4.6. Bubble positions during coalescence final model .... 49
4.7. Coalescence of same size bubbles final model............51
4.8. Coalescence of a larger trailing bubble final model . . 52
4.9. Glycerin coalescence . ..................................54
4.10. Water coalescence.......................................... 56
4.11. Numerical simulation of bubble coalescence in a pipe . . . 57
ACKNOWLEDGMENTS
I wish to thank John Trapp for the thesis topic and for the use
of the existing computer, program. I would also Tike to thank Gary
Young for his support and confidence throughout the project.
NOMENCLATURE
Variables
p density
V Volume
u velocity
U terminal velocity
t time
x horizontal distance or position
y vertical distance or position
g gravity
f wall friction coefficients
fp drag component
fu added mass component
a void fraction
c speed of sound
r radius (equivalent of a sphere)
D total drag
CD drag coefficient
du equivalent diameter
P pressure
B gd2P/uM
A area
a surface tension
fj, dynamic viscosity
n(s) percent of bubble in grid cell
LHS left hand side
RHS right hand side
a width of bubble
b height of bubble
L' distance from end of primary wake to following bubble
u velocity with u,v,w as the i,j,k components
Subscripts
b discrete (bubble) phase
pw primary wake
wc wake center line
sw secondary wake
f free stream
p primary
s secondary
Superscripts
n previous time step
n+1 new time step
CHAPTER 1
INTRODUCTION
Background
Two phase systems are essential in many engineering
applications as well as in nature. Rainfall, fog, and fermentation
are all considered as two phase systems. Industrial applications
utilizing two phase flows are numerous. Mechanical applications
include combustion, fluidized beds, nuclear reactor cooling, foams,
jets and sprays. Distillation, absorption and floatation are also
considered to be two phase systems.
Many two phase flows, such as steam and water, may exist in
various flow regimes. Bubbly flow is characterized by discrete
bubbles in a continuous fluid. Slug flow contains very large,
discrete bubbles which nearly fill the cross sectional flow area.
In annular flow a thin liquid film flows along the walls of the pipe
or vessel and steam flows along the core of the pipe. Drop flow is
characterized by drops of liquid which flow along the pipe and are
surrounded by a continuous steam phase. In addition to these flow
regimes, the flow regimes may be superimposed, resulting in
bubblyslug or annulardrop flow.
Of particular importance is the bubbly flow regime. In this
regime void fractions can range from that of a single discrete
bubble to foam or frothy flow, up to the slug flow regime.
Engineering applications in utilizing bubble flow include
bubblecolumns to promote mass transfer, distillation, foams,
sprays, cryogenics and nuclear reactor cooling, just to name a few.
The heat, mass and momentum transfer capabilities of the bubbly
regime are influenced by the coalescence of the bubbles in the flow.
Bubble coalescence plays a significant role in determining the
size, distribution and interfacial area of bubbles in many devices
and systems. In the cooling of a nuclear reactor, the coalescence
of bubbles can be detrimental because the process decreases the
surface area of the bubbles, decreasing the heat transfer
capabilities. In other applications, such as in separators, bubble
coalescence is advantageous since it produces larger bubbles with
higher rise velocities.
The flow structure around single bubbles has been studied for
years. Data have been collected and empirical relationships derived
for the drag force seen by the bubble and its terminal rise
velocity. The flow field around a bubble has been studied to a
lesser extent. The exact structure of the wake(s), its size, shape,
whether it's laminar or turbulent, is the topic of many
investigations.
The use of the modern computational techniques has been
employed in many codes to model the complex two phase flow problem,
by use of the governing mass, momentum and energy equations.
However, the understanding of the small scale processes, such as
2
bubble coalescence, is essential to produce meaningful results on a
macroscopic scale.
Scope
The purpose of this paper is to make a presentation of the
numerical model used in the simulation of the phenomenon of bubble
coalescence and to show its accuracy by comparison to existing
measured data. An existing one dimensional, two phase flow program
was altered to include the effects of bubble interaction, i.e.,
bubble coalescence. The existing code structure included the basic
governing equations (Chapter 2), but did not account for accurate
rise velocities of individual bubbles nor the interaction of the
bubbles themselves. The program is a mechanical code and does not
consider heat and/or mass transfer between the phases or internally
to the phase.
Additions to the existing program included interactions of
bubbles with each other as well as the continuous phase. Bubbles
are individually tracked throughout the flow field for the selected
problem time and area of interest. Accurate modeling of the drag
coefficients was included to predict terminal rise velocities of
bubbles. A wake model was developed to accurately simulate the
approach of a trailing bubble to a leading bubble. The wake model
was developed based on data in the bubbly flow regime. Slug flow,
annular, drop and vapor flow were not modelled.
3
CHAPTER 2
PROGRAM STRUCTURE
Governing Equations
The existing code is a one dimensional program which solves the
governing equations for a separated, two phase flow. The continuous
phase is modelled as an Eulerian system on a fixed grid. The
discrete, or particle phase, is based on a Lagrangian system, as
each particle(bubble) is individually tracked throughout the flow
field.
The two basic equations which govern the particle phase consist
of the momentum and mass continuity equations. A set of each of
these equations is included for each particle. The momentum
equation for a bubble is used in the following form
'bvb Kb Vb b + 'bÂ¥b fD<% uav> Vb
 f'atb
The term on the LHS side of the equation represents the inertia
effects of the bubble. The first term on the RHS of the equation
represents the force on the particle due to the pressure gradient of
the continuous phase. The pressure gradient, as seen by the
particle, is given by
2.2
() =
l3x;b
/3 U ^ ,,3Ul 1
Pav^9 ^3t + uax^
where the inertia and acceleration terms are based on the average
continuous phase properties near the bubble. The second term on the
RHS of the equation represents the gravity force on the bubble. The
third term represents the drag force on the bubble due to the
relative motion of the bubble and the continuous phase fluid. The
fourth term in the RHS of the equation represents the friction force
of the wall of the container on the particle. The last term is the
added mass force. This term arises when the bubble is accelerated
relative to the continuous phase and it sets up a two dimensional
flow around it which possesses kinetic energy. Extra energy is
required to move the particle. The extra energy requirement creates
an extra force seen by the particle.
The governing mass continuity equation for a bubble is
at(W = 0
2.3
which takes the form
n+lw n+1 n n
Pl Vu = Pl Vl
'b 'b ^b vb
The governing momentum equation for the continuous phase is
given by
2.4
+ u!^ = '+ apq f(uav ub} fu
w b
,3U ,,3U^
f(atb + uei>)
2.5
5
where the LHS side of the equation is the substantive acceleration
of the fluid which consists of the local and convective
contributions to the flow. The first term on the RHS is the
pressure gradient of the continuous phase, followed by the gravity
force. The third and fourth terms represent the effects on the
continuous phase from the drag and added mass forces on the
bubble(s). The final term represents wall friction for the
continuous phase.
The continuous phase mass continuity equation takes the form
at + ax u ,D
The void fraction, a, for the continuous phase is defined as
the total volume fraction of the continuous phase fluid minus the
volume fraction of the bubbles in a particular continuous phase grid
(cell), i.e. if the void fraction of the continuous phase is 1,
there are no bubbles present at that grid location. The void
fraction of the continuous phase, a, is based on the sum of the
void(s) caused by the bubbles and is given by:
a = 1 Eql
2.7
where the void caused by a bubble, c^, is defined as
V.
b = V
D vcell
2.8
and q(s) is the percentage of the bubble volume located in a
particular grid (cell). The percentage, n(s), is based on a fourth
order polynomial which spreads the bubble volume over a specified
6
number of cells in order to obtain a smooth interpolation of the
discrete phase into the continuous phase.
Certain fluid properties (density) for the discrete and
continuous phase are determined by the state equation
dP cd d
Program Scheme
The flow of the original code was not altered to include the
modelling of bubble interaction. A brief description of the program
flow is provided below, including the new features. Figure 2.1 is a
flow chart of the program, with altered and new subroutines
indicated.
The basic inputs to the code include: material properties for
both the continuous and discrete phases, initial particle
geometry/position, fluid parameters and grid specifications.
The initial continuous and discrete phase densities are first
determined in subroutine STATE using equation 2.9. The particle
volumes (equation 2.4) and void fractions (equations 2.7 and 2.8)
for the continuous and discrete phases are then computed in
subroutine PMASS. The drag components for each particle are
calculated in subroutine DRAG for later use in subroutine PMOM. The
DRAG subroutine calls several new subroutines (see Figure 2.1) to
determine wake geometry, wake velocities, bubble placement within
the wake, and the overall drag component.
7
INPUT
Fluid Propertto
Coometry
8
The particle velocity is then calculated using equation 2.1 in
subroutine PMOM. The continuous phase velocities are explicitly
estimated in CMOM using equation 2.5 (estimated only, since the new
time pressure, Pn+ is not determined yet, and the old time
pressure, Pn, is used). Next, the pressure is implicitly determined
in CMASS using equation 2.6, with the voids and velocities in the
equation expressed in terms of pressure. The continuous phase
velocities are then corrected for the pressure in subroutine CMOMFN.
At this point, all continuous phase and discrete phase
variables have been determined. The new bubble positions are then
checked in MERGE for a possible coalescence, merged if required, and
final bubble properties are determined.
9
CHAPTER 3
NEW CODE FEATURES
Bubble Hydrodynamics
Bubbles rising in a Newtonian fluid are generally considered to
belong to three broad regimes. The first regime is the spherical
regime, which usually applies to small bubbles (less than 1 mm in
diameter) that are actually spherical in shape or are only slightly
distorted. Ellipsoid bubbles, the second regime, cover medium sized
bubbles which'are usually oblate in shape, having a convex surface
over the entire bubble (viewed from the inside). The last regime,
spherical cap bubbles, applies to larger bubbles that, in general,
are spherical in the fore geometry and have a flat or skirted tail
in the aft. Figure 3.1 depicts the three bubble regimes and the
variety of bubble shapes which may fall in those regimes.
Three dimensionless groups (Clift, Grace and Weber 1978) are
used to establish which regime a bubble falls under: the Eotvos
number, the Reynolds number and the Morton number.
Eo = 2l
a
Re = ^ 3.2
v
M =
pza3
3.3
Spherical Bubbles
Ellispsoid Bubbles
Spherical Capped Bubbles
Figure 3.1. Bubble shapes
11
Figure 3.2, adapted from Clift, Grace and Weber (1978), shows
the relationship between these three dimensionless groups and bubble
regimes. The spherical regime includes bubbles with low Re, low Eo
and variable M. Bubbles are considered to be ellipsoid at
relatively high Re and intermediate Eo, and low values of M. For a
bubble to be considered a spherical cap, it must have both a
relatively high Re and Eo.
Spherical bubbles can be truly spherical in shape, or be
slightly distorted. Hetsroni (1982) considers bubbles to be
spherical when the aspect ratio (width to height) is a maximum of
1.1. Spherical bubbles are dominated by surface and viscous forces
(therefore low Re and Eo), while inertia and gravitational effects
are negligible. At very low Re, a small bubble may be considered as
a rigid sphere and bubble motions are predicted using Stokes law.
In an ultra pure system, internal circulation of the bubble reduces
the pressure and skin drag on the bubble, resulting in an overall
lower drag force of than of a true rigid sphere. In reality,
surface contamination and slight deformations tend to increase the
drag on the bubble. Surface contamination is generally swept to the
rear of the bubble, generating a surface contamination concentration
gradient. This surface concentration gradient causes a surface
gradient in the shear forces that opposes the surface motion,
increasing the drag. With a higher Re, inertial effects become more
important, leading to the onset of bubble deformation.
12
Bubble regimes adapted from Clift, Grace and Weber (1978)
Figure 3.2. Bubble regimes
13
Ellipsoid bubbles are usually considered to be medium sized
bubbles, having an equivalent diameter (diameter of a spherical
bubble with equal volume) ranging from 1 to 17 mm. The complex
interaction of viscous, surface, gravitational and inertial effects
causes a variety of "ellipsoid" shapes and motions to occur.
Ellipsoid bubbles may be oblate, prolate, or may lack for and aft
symmetry and wobble or oscillate while they rise. Ellipsoid bubbles
are characterized by the presence of secondary motion. Two types of
motion exist; primary rigid body motion, where the bubble can rock
from side to side, and have a zig zag or helical trajectory, and
oscillatory motion (shape dilations). The two motions may be
superimposed leading to extremely complex patterns.
Smaller bubbles in the ellipsoid regime tend to be slightly
oblate. The increasing oblateness tends to promote the development
of an attached wake, and eventually wake shedding. Surface
contaminants can also play an important role in wake generation and
wake shedding by damping out internal circulation. Wake formation
generally occurs at a Re of about 20 and wake shedding at Re of 200.
However, if significant deformation has occurred, the onset of wake
generation/wake shedding may occur at significantly lower values.
Alternatively, interfacial mobility in purified systems
significantly delays the formation of an attached wake and wake
shedding.
The onset of the secondary motion coincides with the onset of
wake shedding. Oscillations disrupt the internal circulation
14
patterns, leading to an increased drag, which in turn, decreases the
terminal velocity. The frequency of the oscillations and wake
shedding do not generally coincide. In fact, evidence suggests
(Otake, Tone, Nakao, and Mitsuhasi 1977) that when the frequency of
the oscillations and wake shedding do coincide, a resonance can
occur contributing to bubble breakup.
The spherical cap regime (d0 > 17mm) covers most bubbles of
engineering interest. These bubbles are important in fluidized bed,
liquid and metal processing and underwater explosions. Spherical
capped bubbles may have a truly spherical or ellipsoidal leading
edge and usually have a flat trailing edge. Dependent upon the
fluid and dynamic properties, the trailing edge may take on
different forms: for lower Re and higher M, the bubble may be
dimpled on the trailing edge, and for slightly higher Re, a skirt
may be formed on the trailing edge (See Figure 3.1). The trailing
skirt occurs when viscous forces on the outer rim are sufficient to
over come restraining effect of surface or interfacial tension
forces. The skirt thickness decreases with the length of the skirt
and the length can increase with an increasing Re. The skirts may
be smooth or wavy, but has little influence on the drag (and
terminal velocity) of the bubble. Surface contaminates have little
effect on the drag of spherical bubbles. Spherical capped bubbles
may have a closed primary wake at lower Re, or an open, turbulent
primary wake at higher Re.
15
Drag Coefficients and Transition Criteria
The previous section discussed different bubble regimes
considered for bubbles of different sizes in different liquids.
Bubble shapes are influenced by the size of the bubble as well as
the physical and dynamic properties of the bubble and its
surrounding fluid. One of the most extensive experimental studies
conducted in this area was performed by Peebles and Garber (1953).
In order to properly model the motion of a rising bubble,
accurate predictions must be made of the drag force seen by the
rising bubble. The numerical model for the drag coefficients and
associated transition criteria were based primarily on the work of
Peebles and Garber (1953). Bubbles were divided into the three
previously discussed regimes: spherical bubbles, ellipsoid bubbles
and spherical cap bubbles. As previously mentioned, cylindrical
bubbles, slug flow, annular and vapor flow regimes were not included
in the numerical model.
Spherical Regime
The spherical bubble regime is actually divided into two
subregions. Although all bubbles in this regime are considered to
be spherical, the underlining assumptions for the two subregions
are different. Figure 3.3, Kelley (1989) indicates that the
relationship between Re and the drag coefficient, C^, is decreasing
for increasing Re.
For subregion 1 (known as Region 1.1), the bubble is assumed to
16
Drag Coefficient (Kelley 1989) vs. Reynolds number, showing
bubble regimes.
Figure 3.3. Reynolds number vs. drag coefficient
obey Stokes Law for a solid sphere. Stokes solution is one of the
oldest known solutions for a creeping motion. A sphere is assumed
to be in a uniform flow, with the normal and tangential velocity
components to the bubble surface are equal zero (no slip
boundaries). Inertia forces are considered to be small in
comparison to viscous forces and gravity is assumed to be the only
external force acting on the bubble. Stokes solution solved the
NavierStokes equations and continuity for an incompressible fluid
where inertial terms were neglected, i.e.,
grad p = /iV2w 3.4
div u = 0 3.5
where: u = ui + vj + wk
The total drag over the bubble was obtained by integrating the
pressure distribution and shear stress over the surface of the
bubble yielding:
D = 6jrrbu5 3.6
Using the customary form of the drag coefficient (obtained by
dividing the drag force by the kinetic energy (dynamic head) and the
cross sectional area),
D
C 3.7
0.5^rb2ub2p
and combining it with equation 3.6, the drag coefficient for a
solid sphere was found to be:
CD = 24/Re 3.8
The Re in this case is defined as Re = 2rbpub/V.
18
The terminal velocity of a bubble in this region is found by
equating the buoyancy force and equation 3.6.
2Vg(' V
u =
9m
3.9
The upper limit for the Stokes region was experimentally found
by Peebles and Garber (1953) to be Re = 2.
The second subregion (Region 1.2) also assumes that the bubble
is of spherical shape, but is tending to become ellipsoid. It was
observed by Peebles that the drag coefficient for this subregion was
slightly lower than the predicted drag coefficient from Stokes
solution. The drag coefficient, similar to the Stokes drag, takes
into account both form and skin drag and was experimentally
determined by Peebles and Garber (1953) to be:
18.7
CD "
Re
0.68
3.10
The terminal velocity of this region was determined by equating
equation 3.10 with the assumed drag coefficient for the bubble's
terminal velocity. This drag coefficient was determined by assuming
that at terminal velocity the buoyancy force on the bubble, VApg,
was equal to the drag force (equation 3.7) resulting in
8 rb9
D 3 Ub2
3.11
Combining equations 3.10 and 3.11, the terminal velocity is
given by:
19
3.12
U = 0.33g0,76(Â£)52rb1,28
The upper limit for region 1, in terms of a transition Re, was
determined based on the drag coefficient relations for Region 2 to
obtain a smooth transition. Bubble behavior in Region 1 depends on
bubble size as well as liquid and gaseous properties. Bubble
behavior is fairly easy to predict, up until a point where the
relationship between the viscous, surface, gravitational and inertia
forces cause the bubble to cease behaving like a solid sphere, and
begin to distort.
The drag coefficient used in the numerical model for Region 1
was taken as the maximum of equations 3.8 and 3.10.
Ellipsoid Regime
Ellipsoid bubbles are also broken into two subregions. Figure
3.3 shows that in region 2, the increase with increasing Re. In
the spherical regime, the drag coefficient was not radically
different than the drag on a solid sphere. Figure 3.3 shows the the
drag coefficient departs radically from the solid sphere drag in the
ellipsoid regime. As the bubble distorts, secondary motions begin
to occur, and the drag on the bubble is greatly increased.
Peebles and Garber (1953) determined the drag coefficient in
this region by first making some assumptions about the terminal
velocity. The terminal velocity was assumed to be some function of
20
the bubble radius, liquid viscosity, density, surface tension and
gravity. Through dimensional analysis the following relationship
was obtained:
3.13
where a and d were determined from experimental data to
be a = 1.91 and d = 0.5. Solving for the terminal velocity
3.14
Equation 3.14 indicates that the terminal velocity in this
region is independent of the fluid viscosity. The drag coefficient,
found by substituting equation 3.14 into equation 3.11, is then
To define the transition between region 1 and region 2 the drag
coefficient relations were equated and solved for the governing Re:
From Peebles work, bubbles in Region 2.2 showed a marked
tendency toward a constant terminal velocity, independent of bubble
size. Figure 3.3 shows that the drag coefficient in this subregion
still increase with increasing Re, but not as sharply as in Region
2.1. The fluid viscosity is still assumed to have negligible effect
on the terminal velocity in this region.
CD = 0.0275 (2^) Re*
3.15
Re = 4.02fe)
*,0.214
3.16
21
Using the same procedure as before, Peebles and Garber (1953)
determined a relationship for the terminal velocity for this region:
Ub = 1.18(2S)025 3.17
Similarly, the associated drag coefficient is
CD = 0.82G1'*25Re 3.18
where G, =
1 pa3
The relationship for the drag coefficient determined by Peebles
and Garber (1953) was found to give a terminal velocity which was
too low. This was probably because the influence of the container
walls was not considered in the development of the empirical drag
relations (Harmathy 1960). Therefore, the relationship suggested by
Harmathy (1960) was used to determine the drag for this region:
CD = 0.488 G10,25Re 3.19
The drag coefficient for region 2 used in the numerical model
was taken as the maximum of equations 3.15 and 3.19.
Studying equations 3.15 and 3.19 one can see that the
combination of the Re and the dimensionless group Gj effectively
cancelled out the viscosity. Therefore, the motion of the bubble in
this region is independent of viscosity. Physically, this indicates
that the drag force is governed by distortion and the swerving
motion of the particle, and the change of the particle shape is
toward an increase in the effective cross section (Ishii and Chawla
1979).
22
Upon closer examination of the equations in region 2.1, if
equation 3.13 is rearranged, it is seen that it is equivalent to a
constant ratio between inertial and surface forces which acts to
deform the bubble.
Spherical Cap Regime
In the spherical cap region, the drag coefficient is assumed to
be constant with a value of 8/3. The transition to this region was
determined by equating the drag coefficient from region 2.2
(equation 3.19) and the 8/3 value, resulting in
Re = 5.47GJ0*25 3.20
The drag coefficients and associated transition criteria
summarized in Table 3.1.
Bubble Geometry
Relationships for bubble shape are considered for use in flow
visualization studies. For calculational purposes, all bubbles are
considered to be spherical, their volumes calculated from mass
continuity (equation 2.4), and equivalent radius determined by
rb=(g)0'333 3.21
Equation 3.21 is used to determine the size of bubbles in
region 1. Ellipsoid bubbles are assumed to maintain a constant
aspect ratio a/b, (See Figure 3.1), although in reality they are
distorting from a spherical shape to a spherical cap. From the work
23
Shaoe Reaime Draa Coefficients Transition Criteria
C Re
Spherical (maximum)
18.7 CD Re0.68 Re = 4.02 (^Â£)'0'214
CQ = 0.0275 (2^) Re*
Ellipsoid (minimum) CD = 0.488 Gj0'2^ Re = 5.47GJ025
Spherical Cap 8 3
Table 3.1. Summary of Drag Coefficients and Transition Criteria
24
of Kelley (1989), ellipsoid bubbles distort from an aspect ratio of
1 to an aspect ratio of 4.02 for spherical caps. Therefore, an
ellipsoid bubble is assumed to have an aspect ratio of 2, and the
equivalent diameter is
d = (4a3)0*333 3.22
t
Spherical cap bubbles, having an aspect ratio of 4.02 have and
equivalent diameter of
j 3a3 a3 *0.333
a M.02 4.023'
3.23
Bubble Coalescence
The coalescence of bubbles is important in determining the
size, distribution and interfacial area of bubbles in gasliquid
systems. Bubble coalescence influences the bubble induced flow and
mass transfer rate.
Bubble coalescence occurs in three stages:
1. A trailing bubble enters into the wake, or partially into
the wake of a leading bubble.
2. Since the relative velocity of the fluid in the wake is
greater than in the surrounding fluid, the trailing bubble
experiences less drag in the wake and approaches the
leading bubble until they collide.
3. After the collision, only a thin film separates the
bubble. A film thinning process takes place until
coalescence occurs.
Only the first two stages of bubble coalescence are considered
for the purpose of this paper. The wake structure behind rising
25
bubbles is clearly important in the coalescence process. The wakes
behind spherical, ellipsoid and spherical cap bubbles were briefly
discussed earlier. A single model was established for the wake
structure behind all bubbles for all three regimes. The model is
based on the wake behind a spherical cap bubble. Spherical capped
bubbles were selected due to their importance in many engineering
applications (and because of the abundance of data for this regime).
Flow visualization studies have clearly indicated the presence
of a wake behind isolated spherical capped bubbles. The structure
of the wake, is composed of two parts: a primary wake and a
secondary wake.
The primary wake is an attached, closed (at low Re) wake that
exhibits a recirculation pattern, or toroidal movement in the wake.
The primary wake is separated from the external flow by a single
streamline. In previous studies, the primary wake shape has been
assumed to have the shape equivalent to completing the sphere dr
i
spheroid of the spherical cap (giving a spherical or oval shape), or
as having a conical or exponentially decaying shape. The velocity
of the recirculating fluid contained in the primary wake may have an
upward velocity that is greater than the rise velocity of the bubble
(relative to a fixed frame). The primary wake structure is
important in the coalescence process of bubbles in close range
(within a few bubble radii).
The recent study by Bessler and Littman (1987) provided a good
description of the wake structure behind spherical cap bubbles.
26
Bessler and Littman (1987) investigated the wake structure behind
circular capped bubbles using aspirin traces. They discuss the
general wake structure adjacent to the primary wake as a free shear
layer which contains large scale vortices (See Figure 3.4) generated
near the bubble rim that remain essentially stationary in the lab
frame. These vortices dissipate energy, growing in size. In more
viscous fluids, the shear layers develop more rapidly and the size
of the primary wake decreases. Bessler and Littman (1987) found the
shape of the primary wake to be an airfoil shape with a cusped tail
due to the growth of the large scale vortices in the free shear
layer. They found the length of the wake to be about six radii of
the leading bubble for a glycerin solution and "slightly longer" for
water.
The secondary wake represents the general bulk movement of the
continuous fluid, due to the momentum defect caused by a rising
bubbles passage. It is the secondary wake motion which is
responsible for the coalescence of bubbles with large separation
distances.
Wake Structure
Primary Wake
Two primary wake models were investigated. Both models were
based on the work of deNevers and Wu (1971). Two wake shapes were
suggested in this study; exponential and conical (see Figure 3.5).
The exponential wake shape was of the form
27
Spherical Capped Bubble
Figure 3.4. Wake structure behind a spherical cap bubble
(Bessler and Littman 1987)
28
3.24
where: R = primary wake radius
Rrw = bubble equivalent radius
xD = distance between bubbles
k = decay constant, determined ,
experimentally to be 0.28 cm1
The conical primary wake assumed a wake length based on the
equivalent radius of the leading bubble. The length of the wake,
*
Lw, was determined using a dimensionless wake length, Lw which
*
remains constant for a particular fluid. For water, l_w = L^R^ =
*
10, as determined by deNevers and Wu (1971) (L = 7 and k = 0.37
w
for glycerin).
Two models were considered for the velocity structure in the
primary wake: 1) the velocity throughout the primary wake was
identical to that of the leading bubble and 2) the toroidal vortex
motion, described earlier, was assumed to give a center line
velocity that initially rose to a value greater than that of the
leading bubble (relative to a fixed frame) then decreasing to the
velocity of the leading bubble.
The toroidal motion model was based on measurements of the
primary wake velocity by Bhaga and Weber (1980) where the vortex
center was found to be located approximately at a distance 60% of
length of the primary wake. The velocity in the wake was assumed be
the same as the center line velocity, increasing from the tail of
the primary wake (which is moving at the same speed as the leading
bubble since the wake is attached) to a value of 1.5 the bubble
30
velocity, u^, at the vortex center. The velocity then decreases
back to (see Figure 3.6).
Secondary Wake
The secondary wake model was based on the asypmtotic wake
solution described by Batchelor (1967) where a long distance
downstream from a body, the wake obtained is independent of the
geometry of the body and is determined by the drag on the body. The
wake velocity along the center line decreases as the reciprocal of
the distance behind the body; and is distributed as a Gaussian curve
in the radial direction. The wake induced flow is taken from
Stuhmiller, Ferguson and Meister (1989) and is
uwCx,r) = uwc(x) e2 3.25
where: u (x,r) = wake velocity at r,x
u (x) = center line wake velocity
x = axial distance
E = radial distance
r = wake radius at x
Several relationships for the center line velocity of the
secondary wake were investigated. The first center line velocity
relationship considered (Stuhmiller, Ferguson and Meister 1989) was
V %(0.20 + 0.12#) + O.Ol#))'1 3.26a
b b
The relationship was based on the collection of center line
wake velocity data from several publications.
31
Figure 3.6. Primary wake velocity
32
Additional center line velocity relationships (listed below)
were given by Crabtree and Bridgewater (1971), Bhaga and Weber
(1980) and Batchelor (1967), respectively.
u. = 6.2
wc
Vb
AX
ubB
WC
( 24(^)'7 + (4)07 )>
3.26b
3.26c
wc
ubdB
AX
3.26d
The radius of the secondary wake was determined by using the
conservation of momentum flux. Far downstream of the leading
bubble, the total momentum is related to the total drag on the
bubble. It is assumed that a long distance downstream of the
bubble, pressure differences in the longitudinal and traverse
directions are small. Therefore, the loss of momentum due to
friction
J = p J ufu. u IdA 3.27
r J cw' b cw'
is equal to the total drag on the bubble
D = 0.5CDpAub2 3.28
The drag coefficient was determined by assuming the buoyancy
force and the drag force were equal (Equation 3.11).
Integrating equation 3.27, assuming that u 2 is small, yields
cw
0 'ubucw'rsw2 329
33
Equating 3.28 and 3.29 gives a relationship for determining the
secondary wake radius:
2 
CDub2
sw u
3.30
cw
Therefore, at distances far downstream, the secondary wake
velocity may be determined by using equations 3.25, 3.26, and 3.30.
At axial locations closer to the leading bubble, the fluid
motion is extremely complex and the assumption of negligible
pressure differences may not apply. Equation 3.30 cannot be applied
to determine the secondary wake radius. The axial location above
which the conservation of momentum flux is not applicable was
assumed to be at the point where the primary wake was of no
influence, i.e. at the tail of the primary wake.
A simple model for the secondary wake radius above this point
was used, based solely on the geometry of the leading bubble. The
radius of the secondary wake at Lw was first determined using
equation 3.30, assuming pressure differences were negligible. The
radius of the secondary wake above that point was then determined
assuming a linear relationship between the radius and L and the
bubble radius.
The velocity structure of the fluid located in the secondary
wake above Lw was considered separately for the primary and
secondary wake. No attempts were made to model the vortices in the
shear layer. The velocity along the streamline of the primary wake
is equal to the velocity of the leading bubble. The velocity at the
34
"edge" of the secondary wake approaches the free stream velocity,
therefore, the velocity in the secondary wake above l_w was assumed
to have a profile similar to equation 3.25, decreasing outwardly
from the velocity of the leading bubble. The complete model of the
wake structure, assuming a conical primary wake shape, is shown in
Figure 3.7.
Drag on a Bubble in a Wake
The total drag on a bubble caught in, or partially caught in,
the wake of another bubble is determined based on the various fluid
velocities "seen" by the bubble. The location of the following
bubble within the wake or wakes is determined first. The area of
intersection of the bubble frontal area and the wake are then
calculated. Finally, the area ratio of the area of intersection to
the bubble frontal area is calculated (the sum of all area ratios is
equal to one for each bubble). Area ratios were used in an attempt
to smooth or eliminate discontinuities in the total drag on the
bubble as it entered the region of the wake(s). If a following
bubble is effected by the wake of more than one leading bubble, only
the leading bubble with the highest (combined primary and secondary)
area ratio is considered. If a following bubble is completely within
the wake of one or more leading bubbles, the bubble closest to the
following bubble is considered.
35
Figure 3.7. Wake structure behind a spherical cap bubble
(code model)
36
The area ratios are used to determine the components of the
total drag force seen by the bubble. A drag component is calculated
for each fluid stream seen by the bubble. A drag coefficient for
each component is calculated using methods described earlier, and
based on the particular fluid stream velocity. The fluid velocity
seen by a bubble intersecting the primary wake is assumed to be that
of the center line velocity (Figure 3.6). The fluid velocity seen
by a bubble intersecting the secondary wake is a calculated average
velocity based on the frontal area of the bubble in the wake and its
location in the wake.
Drag components for a bubble partially within the primary wake
and secondary wake of a leading bubble, as well as partially in the
free stream, are calculated below.
The drag components are used in the momentum equation (2.1), in
subroutine PMOM to solve explicitly for the bubble velocity.
The merging of two or more bubbles is assumed to occur if two
bubbles touch. The film thinning process is not modelled. The
3.31a
3.31b
3.31c
where the subscripts f,p and s refer to
the free stream, primary wake, and secondary
wake, respectively.
Bubble Merging
37
equivalent radius, equation 3.20, was used to determine if a bubble
touched or overlapped another bubble.
The new, "merged" bubble required as estimate of its new mass,
volume and velocity. A simplistic model of the conservation of
linear momentum and mass was used to make an initial guess of the
new velocity:
Vbl + "2ub2 (ral + m2)ub new 3'32
The size of the merged bubble is determined by assuming the new
density is a weighted average of the individual bubbles and total
mass is constant
V = />av/( m1 + m2) 3.33
Although equation 3.32 assumes a constant linear motion, a new
horizontal (y) position for the merged bubble is calculated. The
new bubble position is calculated based on the relative volumes of
the two individual bubbles
ynew = (ylVl + y2V2)/(Vl + V2J 3,34
38
CHAPTER IV
RESULTS AND DISCUSSIONS
Comparison between calculated values and experimental data are
presented first to show the accuracy of the numerical models added
to the existing code. The drag coefficient calculations are tested
by comparing the calculated terminal velocity of individual bubbles
and experimental data. Second, wake structure models are tested by
comparing calculated bubble coalescence times with those from
experimental data. Finally, a flow visualization of the coalescence
of bubbles rising vertically in a pipe is presented.
Terminal Velocity
The terminal velocities of bubbles with varying radii were
calculated using the drag coefficient relationships and transition
criteria discussed in Chapter 3. The discrete phase for the
calculations was water and the particle (bubble) phase was air.
Figures 4.1 and 4.2 provide a comparison of calculated terminal
velocities with measurements of Peebles and Garber (1953) and with
that of Haberman and Morton (1953). Also shown in the figure are
the terminal velocities calculated by Peebles and Garber (1953).
The calculated values agree well with the calculated values by
Peebles and Garber (1953), as well as the measured data. Note that
Terminal Velocity, cm
Comparison of code calculated terminal velocities of bubbles
in water and measured velocities by Peebles and Garber (1953).
Also shown are terminal velocities calculated by Peebles and
Garber (1953) using equations in Chapter 3.
Figure 4.1. Terminal velocity of bubbles in water (Peebles and Garber 1953)
Terminal Velocity, cm/s
6j
4
3
2
1E+01=
6:
4
3
2
x
1E+OOH1rr
1E02 234
1E01 2 3 4 1E+00 2 3 4
Bubble Radius, cm
Comparison of code calculated terminal velocities of bubbles in
water and measured velocities by Haberman and Morton (1953).
Figure 4.2. Terminal velocity of bubbles in water (Haberman and Morton 1953)
the calculated data lies above the experimental data for bubbles
with equivalent radii greater the 0.3 cm. These values were
considered to be unrealistically low, as noted in Chapter 3, and a
different drag relationship was used for bubbles in this regime.
Additional comparisons to measured data are shown in Figure
4.2. The measured data was extracted from Haberman and Morton
(1953). Excellent agreement is shown for all bubble radii, except
where the transition from a spherical bubble to an ellipsoid bubble
occurs. A possible explanation is that the water used in the
Peebles measurements was contaminated, resulting in lower terminal
velocities.
Coalescence Times
The wake model was compared to data gathered by Crabtree and
Bridgewater (1971). Crabtree conducted experiments on the relative
motion of vertically aligned air bubble pairs, having volumes from
3 3
10 cm to 40 cm in a 67 percent wt solution of sucrose in water.
The bubble pairs were tested up to 70 cm apart.
The density and dynamic viscosity of the sucrose solution were
given, but the surface tension and speed of sound were not provided,
so estimates of those properties were made. A few test runs were
made to determine the sensitivity of the estimate properties. These
runs showed that varying the surface tension and speed of sound did
not effect the bubble velocities drastically.
42
Original Model
The first wake structure model attempted, used the following
1. Conical/exponential primary wake shape
(L = 10, k = 0.28)
2. Velocity in the primary wake = U.
3. Velocity in the secondary wake as in Eqn 3.23a
Figure 4.3 shows the coalescence of a bubble pair with an
initial separation distance of 25 cm. The leading bubble volume was
3 3
30 cm and the following bubble volume was 25 cm Results from
both the conical and exponential primary wake model are presented.
The leading bubble is unaffected by the following bubble. Its rise
velocity compares well with the experimental data. The following
bubble is initially affected by the secondary wake only, entering
the primary wake at around t = 0.3 seconds. The following bubble
data compare fairly well with the experimental data; however, near
coalescence the calculated velocity for the following bubble is
somewhat less that the velocity indicated in the experimental data.
This could suggest improper modelling of the primary wake.
Although not shown, additional runs were made assuming the drag
component on the following bubble was zero for the area of the
bubble considered to be within the wake. The results of this change
were small, although calculated coalescence times were shorter. The
conical model appears to given better results than the exponential
model.
The coalescence times of bubbles of the same size and for a
smaller leading bubble and larger following bubble are presented in
43
Bubble position vs. time of a trailing bubble (25 cc) approaching
a leading bubble (30 cc), in sugar water, using the original code
model. Measured data from Crabtree and Bridgewater (1971).
Figure 4.3. Bubble positions during coalescence original model
Figures 4.4 and 4.5, respectively. Again, the conical model gave
better results than the exponential. In both cases, the model
predicted the coalescence times closely to the experimental data for
bubbles with an initial separation distance of 40 cm or less. For
bubbles with an initial separation distance of greater than 40 cm
the calculated coalescence time were longer than the experimental
times.
Figure 4.5 gives a more favorable comparison with experimental
data. This was expected since the higher rise velocity of a larger
following bubble would play a larger role in the coalescence than
the wake of the smaller, slower leading bubble.
Final Model
Due to the relatively poor agreement of the data for both the
trailing bubble positions and the coalescence times, variations of
the model were made. First, a comparison of the available secondary
wake velocity relationships was made, selecting a consistent bubble
size and terminal velocity for comparison. The asypmtotic wake
velocity clearly gave too high a wake velocity (this is reported by
several authors). Equations 3.26b and 3.26c gave a similar velocity
profile, while equation 3.26d yielded a slower velocity. After many
comparison runs, Bhaga and Weber's (1980) relationship appeared to
give the most reasonable results. This relationship appeared to be
more appropriate than the others, since the velocity is based on
properties of the fluid and the bubble size, not just
45
Coalescence time vs. initial separation distance for bubbles
of equal volume (25 cc), in sugar water, using original code
model. Measured data from Crabtree and Bridgewater (1971).
Figure 4.4. Coalescence of same size bubbles original model
Coalescence time vs. initial separation distance for a larger
trailing bubble (40 cc) and a smaller leading bubble (10 cc),
in sugar water, using initial code model. Measured data from
Crabtree and Bridgewater (1971).
Figure 4.5. Coalescence of a larger trailing bubble original model
the bubble size. The conical primary wake shape was judged to give
better results and it was retained for the final model. A value of
*
Lw = 7 for the sugar water, as opposed to a value of 10 for pure
water as determined by deNevers and Wu (1971), was found to geve
better results. The primary wake velocity was assumed to follow the
recirculation profile given in Figure 3.6.
Figure 4.6 shows the coalescence of the same bubble pair as in
Figure 4.3, but using the final model. The calculated rise velocity
of the leading bubble compares well with the measured rise velocity
of the leading bubble, indicating the drag coefficient used for the
spherical cap was adequate. Additionally, the linear rise indicates
that the following bubble has no effect on the rise velocity of the
leading bubble.
The calculated positions of the following bubble agree quite
well with measured data throughout its approach toward the leading
bubble. The trailing bubble enters the primary wake of the leading
bubble at about t = 0.6 sec., after which the rise velocity of the
trailing bubble increases somewhat, but not to the extent of the
measured data.
Near coalescence, the calculated position of the trailing
bubble is slightly less than the measured data. Crabtree and
Bridgewater (1971) observed that as the following bubble came within
two or three radii of the leading bubble its aspect ratio increased,
elongating the bubble in the vertical direction. An increase in the
48
Bubble position vs. time of a trailing bubble (25 cc) approaching
a leading bubble (30 cc), in sugar water, using the final code
model. Measured data from Crabtree and Bridgewater (1971).
gure 4.6. Bubble positions during coalescence final model
aspect ratio would allow more of the bubble to enter the primary
wake, decreasing its overall drag, and increasing its velocity.
Figures 4.7 and 4.8 are for bubble coalescences in sugar water
over larger separation distances. For bubbles of equal volume
(Figure 4.7) the calculated coalescence times are much improved over
the previous model. The coalescence times match the measured data
quite well within 60 cm separation distances. The model predicts
slower coalescence time at larger distances, approximately 20%
slower than measured data at 70 cm. The coalescence of a larger
trailing bubble and smaller leading bubble (Figure 4.9) were also
improved using the new model. Coalescence times are predicted
reasonably, even up to 80 cm separation distances.
Figures 4.6 4.8 gave comparisons of calculated vs measured
data from bubble coalescences at fairly large distances with
reasonable success. At closer distances the importance of the
primary wake is increased. Few experimental data were found for
comparisons of coalescences at close distances. Experimental data
from the work of deNevers and Wu (1971) provided data in two liquids
for initial bubble separation distances of five radii. The two
liquids considered were glycerin and pure water.
The mechanism used to release the bubbles was a horizontal
sparger with a single circular hole. As bubbles were released into
the stagnant fluid, the motion of the bubble chain tended to cause
the general bulk movement of the fluid to have an upward motion.
This movement caused the terminal velocities of the bubbles to be
50
Coalescence time vs. initial separation distance for bubbles
of equal volume (25 cc), in sugar water, using final code
model. Measured data from Crabtree and Bridgewater (1971).
Figure 4.7. Coalescence of same size bubbles final model
Coalescence time vs. initial separation distance for a larger
trailing bubble (40 cc) and a smaller leading bubble (10 cc),
in sugar water, using final code model. Measured data from
Crabtree and Bridgewater (1971).
Figure 4.8. Coalescence of a larger trailing bubble final model
increased by as much a 30%. The upward motion, as it will be
discussed later, could cause discrepancies between calculated and
measured data.
The first coalescence event considered was for air bubbles in
glycerin (Figure 4.9). The calculated motion of the bubble is
somewhat different than the measured, although both coalesce in
approximately the same time. The velocity of the measured bubble
increases as the separation distance decreases. The calculated
velocity increases initially, but then appears to be constant.
Effects on the bubble motion by the release mechanism are unknown.
However, the primary wake model appears to have adequately predicted
the bubble motion in this case.
The available data on coalescence times for pure water, as
opposed to sugar water or other viscous fluids, were somewhat
limited. DeNevers and Wu (1971) obtained data from the coalescence
of two bubbles with equal radii in pure water. Attempts were made
to match the data, but with limited success. DeNevers and Wu (1971)
pictured the bubbles schematically to be spherical capped bubbles.
However, the equivalent radii of bubbles for the given bubble
volumes, fell in the ellipsoid range. In addition, the author
discussed a rocking motion of the bubble (causing the data to
scatter) which is characteristic of ellipsoid bubbles. The rocking
motion, as discussed next, causes difficulties in predicting the
coalescence.
53
Dimensionless Time, t* = Ut/r
Dimensionless Distance, x* = del x/r
Time vs. bubble position for a coalescence in glycerin. Measured
data from deNevers and Wu (1971).
Figure 4.9. Glycerin coalescence
Figure 4.10 shows the calculated and measured coalescence data.
According to the author, the trailing bubble was rocking back and
forth, moving in and out of the wake of the leading bubble. The
bubble was considered to be out of the wake of the leading bubble
ic ic ic
from about t = 15 to t = 20. The scatter of data near t = 0 also
suggests that the bubble may have been partially out of the wake.
If the general motion of the bubble at t =15 were continued,
*
it appears that the bubble would have coalesced near t = 14 or 15.
The calculated data indicated a faster coalescence time, occurring
*
at approximately t = 9. The assumed time of 14 or 15 gives a
coalescence time approximately 55% slower than predicted. Part of
the discrepancy may be due to the release mechanisms problems, but
more likely it is due to the oscillatory motion of the bubble
causing it to move in and out of the wake.
Flow Visualization
Figure 4.11 shows a flow visualization of bubble coalescence in
a 20 cm diameter pipe over a twelve second time period. At the
initial state, t = 0 seconds, three spherical bubbles exist in the
lower end of the pipe. Bubbles of various sizes and shapes are
released over the first five seconds. The stages of coalescence of
the bubbles are shown at times t = 4 seconds, t = 8 seconds, and t =
55
Dimensionless Time, t* = Ut/r
Dimensionless Distance, x* = del x/r
Time vs. bubble position for a coalescence in water. Measured
data from deNevers and Wu (1971).
Figure 4.10. Water coalescence
I
0 seconds
6 seconds
12 seconds
10 m
5 m
o
o
0
o
mi
Figure 4.11. Numerical simulation of bubble coalescence in a pipe
57
12 seconds. At t = 12 seconds the bubbles have coalesced to the
point where the onset of slug flow has begun, and the problem is
terminated.
Discussion
The empirical relationships derived by Peebles and Garber
(1953) and by Harmathy (1960) resulted in terminal velocities which
matched measured data very well. Accurate predictions of the drag
force (and terminal velocity) on a single bubble were essential to
the accurate prediction of bubble coalescence. The wake model
performed well in predicting the coalescence of bubbles from both
long and close separation distances.
The primary wake model, which at first was assumed to be a
major contributor to the coalescence effect, was modelled (and
remodelled) in great detail. However, in the review of the results
and the experience gained in the numerous test runs, it is not the
dominating process causing bubble to coalesce. The primary wake
affects only the end of the approach of a trailing bubble to a
leading bubble. It is the secondary wake which plays the dominant
role in bubble coalescence. The structure of the primary wake, the
internal circulation and the primary wake velocity are, however, an
interesting aspect of rising bubbles which help determine the flow
structure around the bubble. Many authors have produced excellent
photographs of the wakes behind spherical capped bubbles.
58
In any numerical model one must make a decision on the most
important effects of the actual physics, and the level of detail
which is necessary to adequately predict reality. Dependent upon
the required accuracy, it is possible that the primary wake model
may be sacrificed to decrease the complexity of large problems. The
modelling of each bubble in the Lagrangian method, as opposed to a
homogeneous twophase model, increases the complexity of the problem
and requires extra computational time, especially when large numbers
of bubbles are tracked. Therefore, any method of decreasing the
complexity of the model would clearly be advantageous.
The difference in scales between the continuous phase and
discrete phase can cause additional problems in modelling twophase
flow problems. The level on a more microscopic level of an
individual bubble required a different scale than for the continuous
phase (failure to recognize this fact or properly account for it can
give rise to executional errors in the program which may only be
solved by many late night test runs accompanied by caffeine
stimulants).
Additions to the code could be made to extend not only the
range of bubbles regimes (include slug flow, for example) but also
to handle more realistic scenarios. Breakup of bubbles could be
included in the model, as well as 2D effects of vessel walls. The
majority of the available data on bubble coalescence deals with
fairly viscous fluids, such as sugar water and glycerin and the
empirical relationships for the secondary wake velocity are laminar
59
in nature, decreasing by x"1. Problems for bubbles in turbulent
regimes have not been explored.
60
REFERENCES
Batchelor, G.K. 1967. An Introduction to Fluid Dynamics. London:
Cambridge University Press.
Bessler, W.F., and H. Littman. 1987. Experimental Studies of Wakes
Behind Circularly Capped Bubbles. Journal of Fluid Mechanics
185:137151.
Bhaga, D., and M.E. Weber. 1980. InLine Interaction of a Pair of
Bubbles in a Viscous Liquid. Chemical Engineering Science
35:24672474.
Clift, R., J.R. Grace, and M.E. Weber. 1978. Bubbles. Drops and
Particles. New York: Academic Press.
Crabtree, J.R., and J. Bridgewater. 1971. Bubble Coalescence in
Viscous Liquids. Chemical Engineering Science 26:83985.1.
deNevers, N., and J.L. Wu. 1971. Bubble Coalescence in Viscous
Fluids. AIChE Journal 17:182186.
Haberman, W.L., and R.K. Morton. 1953. David Taylor Model Basin
Washington. Report number 802.
Harmathy, T.Z. 1960. Velocity of Large Drops and Bubbles in Media
of Infinite or Restricted Extent. AIChE Journal 6:218288.
Hetsroni, G. 1982. Handbook of Multiphase Systems. New York:
McGrawHill.
Ishii, M., and T.C. Chawla. 1979. Local Drag Laws in Dispersed
TwoPhase Flow. Argonne National Laboratory NUREG/CR1230.
Kelley, J.M. 1989. Towards a Mechanistic Model for Multiphase Flow
The FourField Approach. Batelle, PNL (Draft Report).
Otake, T., S. Tone, K. Nakao, and Y. Mitsuhasi. 1977. Coalescence
and Breakup of Bubbles in Liquids. Chemical Engineering
Science 32:377383.
Peebles, F.N., and H.J. Garber. 1953. Studies in the Motion of Gas
Bubbles in Liquids. Chemical Engineering Progress 49:8897.
Stuhmiller, J.H., R.E. Ferguson, and C.A. Meister. 1989. Numerical
Simulation of Bubbly Flow. Electric Power Research Institute
NP6557.
61
APPENDIX
uu
PROGRAM OlSCtM
THIS PROGRAM SOLVES THE EQUATIONS FOR A MIXTURE OF A CONTINUOUS
AMO DISCRETE PHASE. TIC CONTINUOUS PHASE IS ON A FIXED KSM.
THE DISCRETE PHASE IS LACRANGIAN. THE CODE IS AS EXPLICIT AS
POSSIBLE WITHOUT A SONIC LIMIT.
INCLUDE (TT COMMON)
** INDEX STRUCTURE A PICTURE
aaa a a a a
1 2 VOLltCS 1ST a ISTP1 LSTP2
* a a
1 2 . S JUNCTIONS LST LSTP1 LSTP2
! BC 1REAL VOLUMES AMD JUNCTIONS! 1C I PHONY
a..............
ItCHIT DATA
CALL INPUT
ms PRELIMINARY INITIALIZATIONS
NCYC I
TIIC 1.1
IPRINT
C
c <
C Mia LOOP RETURN TO CONTINUE TIIC CYCLE
C
1M CONTINUE
C
IF (NCVC.EQ.B)GO TO 791
C MHI .......
C M ICRGE BUBBLES THAT TOUCH OR OVERLAP
C aaa
CALL ICRGE
C 
c aaa STATE PROPERTIES
c aaa
C
791 CALL STATE
IF (NCYC .EQ. 1) TWN
DO ISO I 1 1, IP
IF(IBUB(X).EQ.I)GO TO 15
RMOPAN(I) RHOPA(I)
VOLPN(I) VOLP(I)
151 CONTINUE
END IF
C
c
c aaa PARTICLE MASS EQUATIONS SOLVED FOR VOLltC
C (ALSO FINAL SHAPE DISTRIBUTIONS AMD VOIDS CALCULATED)
C aaa.....................................
C
CALL PKA5S
C
C IF (NCYC.EQ;0)GO TO 796
C aaa ................
C ICRGE BUBBLES THAT TOUCH OR OVERLAP MOT FIRST TIIC
c aaa
C CALL MERGE
CONTINUE
OUTPUT EDIT IF CALLED FOR
C
CALL OUTPUT
WRITE (6,9021) NCYC
DO 029 N 1,NP
XF(1BUB(H).EQQ)G0 TO 629
WRITE(4,9I23)N,XP(H),YP(N),VELP
029 CONTINUE
C
C aaa ............................
C LOGICAL BEGINNING OF THE TIIC STEP CYCLE
C aaa .................................
C
C
NCYC NCYC 1
TIIC TIME OT
XIN 1.1
IF (NCYC .LT. XWCYC) XIIt.0
C
C aaa SET UP CIO TIIC VALUES
C MW 
C
211
SB!
411
C
510
C
OO SIM 1, LSTP2
PRESN(K) PRES(K)
RHON(K) RHO(K)
VOIDN(K) VOXD(K)
DO 211 I 1, IP
XF(IBUB(X).EQ.O)GO TO 2M
VOIDPN(K,I) VOIDP(R,X)
PERCTN(K,I) PCRCT<*,I)
CONTINUE
CONTXMC
DO AM J 1, LSTP2
VELN(J) VEl(J)
VELSFN(J) VELSF(J)
CONTXMC
DO SOQ I 1, HP
IF(XBUB(1).Â£Q.0)G0 TO 501
RHOPAN(X) RHOPA(I)
VELPN(X) > VELP(I)
XPN(I) a XPU)
VOLPN(X) VOLP(X)
CONTXMC
63
C *
C GET INTERPOLATED CONTINUOUS HELD VARIABLES FOR PARTICLE EQUS.
C m*  
C
DO 010 1 1. NP
C
IFdBUB!l).EQ.O)CO TO 001
RHOAVE(I) 0.0
00 610 X I, LSTP1
RHOAVE(I) RHOAVEd) KRCT(K,X)*RHO(K)
600 CONTIWC
C
VELAVE(I) 0.0
ACCAVE(X) 0.0
DO 700 J 2, LSTPI
VCLAVEd) VELAVEd) PERCTJCJ,I)*VELSF!J)
ACCAVE(I) ACCAVE(I) PERCTJ!J,nACCÂ£l!J)
700 CONTINUE
C
800 CONTINUE
C
C ................................
C CET ALL DRAG FORCES
C IM  ........
c
CALL DRAG
C
c **...
C * ADVANCE PARTICLE NONENTW EQUATIONS FOR VELOCITY AND POSITION
C Ma  ...
C
CALL PNOM
C
C 
C ADVANCE CONTINUOUS MOMENTUM EQUATIONS FOR VELOCITY
C mmm
C
CALL CNOR
c
c aaa..................................
C * CONTINUOUS HASS EQUATIONS SOLVED FOR PRESSURE
C M  
c
CALL OUSS
C
c aaa ............. 
C m FINAL CONTINUOUS HONENTUH EQUATION ADVANCEMENT FOR VELOCITY
C MMN 
C
CALL CNOtEN
C
c **  . 
C CONTINUE TINE CYCLE LOOP
C aaa...................................
C
CO TO ltl
C
C aaa 
C mm FGRHATS
C 
c
9020 FORMAT!/3X'FOR CYCLE M5,' BUB X, Y POSITION A VELOCITIES:')
9023 FORIWTCSX,I3,10X,E12.S,3X,Â£1Z.S,3X,E12.S)
C
C
C
c
c
c
c
c
e
c
c
c
c
c
c
c
END
SUBROUTINE OUSS
THIS SUBROUTINE USES TtC CONTINUOUS PHASE NASS BALANCE WITH DC
VOIDS AND VELOCITIES EXPRESSED IN TERNS OF PRESSURE AND SOLVED
IMPLICITLY FOR PRESSURE
INCLUDE! TTCOI MOM)
GET OONORED PROPERTIES FOR NASS EQUATION
DO 101 J2, LSTPI
K J
IF (VEUJ) .67. 0.0) TIEN
RHOJ(J) RHONU1)
ELSE IF (VEUJ) .LT. 0.8) DEN
RHQJ(J) RHON(K)
ELSE
RHQJ!J) 0.S*!RHON!Kl)flHOM!K))
ENDIF
100 CONTINUE
BUILD PRESSURE EQUATION
DO 301 K 2, 1ST
J K
TEPHP VOIDJ! J*l)*RHQJ! J+1)"0T*ARÂ£A
TERMM VOIDJCJ )*RHOJ!J >OT*AREA
TERN XINVOLbRHON(R)>DAOP(K)
A(K) TERMM*V0P!J ) 0.00*TERN
C!K) TGWP"VOP! J*1) 0.00*TERN
B!K) VOL*VOIDN!K)"ORDP A!K) C(K)
1.30TERN
D!K) *TERMPVEL!J*l> TERMM*VEUJ )
DO 200 X &, HP
XF!IBUB(I).EQ.O)CO TO 200
D!K) D!K)
 RHON(K)*VELP!X)*DT*AREA*(VOXDPJ(Jld)VOXDPJ(JX))
208 CONTINUE
CHANCES DUE TO BOUNDARY CONDITIONS
IF(K .EQ. LST) D(X) DU) C!K)*!PRPRESM(LSTP1))
308 CONTINUE
TRIDIAGONAL PRESSURE SOLUTION FOR INCREMENTS DP(K)
(A(2) AND C(LST) NOT USED)
CALL TRIDIA
C
C UPDATE PRESSURE
DO 600 X 2, LST
PRESIX) 0P!K) PRESN(K)
400 CONTINUE
C
C BOUNDARY CONDITIONS
C LEFT IOUWARV
PRES!1) PRESI2)
64
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
RIGHT BOUNDARY
PRES(LSTPl) PR
DP(lSTPl) 0.0
RETURN
END
subroutiic CHON
THIS SUBROUTINE USES THE CONTINUOUS PHASE LINEAR MOMENTUM EQUATIONS
TO GET AM EXPLICIT ESTIMATE OF THE CONTINUOUS PHASE VELOCITY
INCLUDE (TTCOttOO
DO 100 J 2, LSTPI
X J
RHOC Q.S*(RHON(Kl) RHON(K))
VOIDC 0.5" VOXON(K))
ARV VOIDC*RHOC*VOL
CONVECTIVE ACCELERATION
IF (VELN(J) .GT. 0.0) T>CN
CONV VEIN( J)"(VELM( J)VELN( Jl) )/DX
ELSE IF (VELN(J) .IT. 0.0) THEN
CONV VELN(J)(VELN(J*1)VELH(J))/DX
ADVANCE CONTINUOUS VELOCITY
NTEST 1
IF (NTEST .EQ. X) THEN
SUN  ARVVELH(J) DT"ARVCONV
OTBARVGRAV
 DTV0IDC"V0L(PRE5NCX)PRESN(X1))/DX
 DTHXFACEP(J)
COEF ARV DT*FVALL(J>
VEL(J) SUM/COEF
VDP(J) DT/(RHOC"DX)
ELSE
SUM ARVaVELMU) DT"ARVCONV DTACCELP(J)
DTMARVMGRAV DTGRAVP(J)
 DTaVOL(PRESM(X)PRÂ£SN(Kl))/DX
 DTXFACEPW)
COEF ARV DTFVALL(J)
VEL(J) > SUN/COEF
VOPU) OT/(RHOCOX)
END IF
GRAVWKJ) DX*(ARV*G8AV GRAVP(J))/VOL
lit CONTINUE
BOUNOARY CONDITIONS
LEFT BOUNDARY'
VEL(2) VL
VEL(I) VEL(Z)
VELSF(l) VELSF(Z)
VDPC2) 1.0
RIGHT BOUNDARY
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
e
c
c
VEL(LSTP2) VELCLSTPl)
VELSF(LSTPZ) VELSF(LSTPl)
RETURN
END
SUBROUTINE CHOtTN
THIS SUBROUTINE MAXES AN ILLICIT CORRECTION TO THE CONTINUOUS
PHASE VELOCITIES (DUE TO DC ILLICIT PRESSURE)
XNCLUDECTT COMMON)
DO 200 J 2, LSTPI
X J
VELCJ) VEL(J) VDP(J)"(DP(X)'*DP(K1))
VELSF(J) VOIDJ(J)"VEL(J)
DO 100 I I, W>
IFdBUB(I) .EQ.O)GO TO ll
VELSF(J) VELSF(J) VOIDPJCJ,I)VELP(I)
100 CONTINUE
ACCEL(J) (VELSF(J)VELSFN(J) )/DT
* VELSFH(J)MVELSFN(J1)
210 CONTIWC
BOUNDARY CONDITIONS
LEFT BOUNTRY
VEU2) VL
VEL(I) VEL(2)
VELSF(l) VELSFC2)
RIGHT BOUNDARY
VEL(LSTP2) VEL(LSTPl)
VELSFUSTP2) VELSF(LSTPI)
RETUM
Eta
SUBROUTIIC DRAG
THIS SUBROUTINE CALCULATES ALL DRAG TERMS. NALL DRAG FOR TtC
CONTINUOUS PHASE* INTERFACE DRAG BETMEEH T>C PARTICLE AND CONTINUOUS
PHASE* AND HALL DRA6 FOR TtC PARTICLE (THIS IS TO HOLD PARTICLE TO
THE HALL BEFORE RELEASED INTO THE FREE STREAM).
XNCLUDC(TTCONMM)
INITIALIZE AREA RATIO AND BUBBLE INDEX
DO 99 N 1,NP
IKOÂ£XR(N)M
99 CONTINUE
C
C CONTINUOUS PHASE HALL DRAG
C
DO 100 J 2* LSTPI
FHALKJ) l.t
101 CONTINUE
C
C PARTICLE HALL DRAG
C
65
o
on oooon
n o
n
cn
o>
*
vi
no nnon
DO 20a I 1, w
IF(IBUB(I).EQ.0100 TO 209
FPMLLdJ 1.9
IF (TI .LT. RELTIH(I)) FPMLl(l) 1.9E20
200 CONTINUE
5992
C
C
c
c
c
Cl24.
C2Z4.
GO TO 5992
EKDIF
Cl24.0(1.0/REP(I))
C2i6.7REP"(0.66>
CSUBD2 DHAX1(C1,C2)
ELSEIF
Cl0.027561*REP( 1 >4
C20. 466*61 (0.25) REP< I)
CSUB02 DNIM1(C1,C2)
ELSE
CSUBD28./3.
RETURN
EHO
SUBROUTINE INPUT
THIS SUBROUTINE IWVTS AND PROCESSES ALL THE PROBLEM DATA,
INITIAL CONDITIONS AMD BOUNDARY CONDITIONS
INCLUDE (TTCOfBOf)
C
c
c
c
c
c
c
c
MTL6 0
BASIC NUMERICAL DATA
OT 9.001
0T 0.01
OX o.os
DX 0.010
NAXCYC 10
NAXCYC 1000
NCELLS 196
NCELLS 199
NCELLS 496
PIPELEMOXFLOAT(NCELLS)
NP 25
*RINT 1
IGSAPHM
NOUT 1
NOEBUG 6
C INITIALIZE BUBBLE FLAG
DO 57 I 1,1*
I6UB(I>1
57 CONTINUE
C
C PERMENANT INDEXES (CEDES FOR LOOP STRUCTURE
LST NCELLS 1
LSTP1 LST 1
LSTP2 1ST 2
LSTM1 1ST l
C
c
CVW
c
c
EQUATION PARAICTERS
PIE 5.14159
NAME PIPE DIAMETER LARGE ENOUGH MOT TO EFFECT BUBBLES W/ D*6 CM
AREA PIE1.0**2/4.
P1PD1A 6.1
P1PDXA *0.2
AREA PIEPIPDIA2/4.6
VOL AREAOX
C CALC TIC MAX EOUIV RADIUS FOR ABOVE PIPD1A W/ SPEWICAL CAPPED
C BUBBLES WITH ASPECT RATIO 4.02, CAP VOL (1/6)6(3A*2*B2)
C RE^ O.IKoOC PIPDIA/(C3.D00/4.02IKOO)(1.D*00/4.02IHOOm5)))
Z t (3.0*00/4.020*06) (l.D*00/4.02D*003) )(I.D*00/3.D*00)
REMAX O.SD*00 P1PDXA l
WRXTE(6,392)REHAX ______
392 FORNATC REMAX (,E12.S,* METERS*)
C
GRAV 9.6
CSUBD 0.4
CADDMS 0.5
CXNERA 1.0
XLEN 1.0
IWCYC 20
C INITIALIZE SHAPE FACTORS
DO 94 M 1,MP
XSHAPE(H) 0
94 CONTINUE ____ _________
CVU INPUT BUBBLE DIAfCTERS (METERS)
C DIA1RENAX0.99990*06
DIA10.005
DIA20.0100
DIAS0.0050
DIA40.0106
DIAS0.0106
DXA60.010I
DIA70.0160
DIA60.0106
DIA90.0600
OIAIOO.OIOO
DIA110.0140
D1A12*0.0400
DXA15*0.0106
DIA140.0100
DIA1S0.U00
DIA160.0606
DIA170.0606
DIA160.0500
DIA190.1006
DXA20O.00S0
DXA210.0106
DIA226.00S
0IAZ5M.0900
DXA24O.0I06
D1A2SI.10#
WRITE (6,) *BUB DXA 1 ,D1A1
WRITE (6,) 'BUB DIA 2  ,DIA2
WRITE (6,) 'BUB DIA 3 ',DZA3
WRITE (6, ) 'BUB DIA 4 .DIA4
WRITE (6,) 'BUB DIA 5 ,DIAS
WRITE (6.) BUB DIA 6 *,DIA6
WRITE (6,) 'BUB DIA 7 "',DIA7
WRITE (6,*) BUB DIA 6 ,DlA6
WRITE (6,) 'BUB DIA 9 *,0IA9
WRITE (6,> BUB 0IA10 '.DIA10
WRITE (6,) BUB DIA11 >',DIA11
WRITE (6,) *BUB DIAI2 ,DlA12
WRITE (6,*) 'BUB DIA13 *,DIA13
WRITE (6,*) BUB 0IA14 ,DIA14
67
c
CVH
WRITE 16,) 'BUB 0IA1S .9IA1S
MITE (6,") 'BUB DIA16 ' .01A16
MITE (6,"> 'BUB 0IA17 *,0IA17
MITE <6,) BUB OUia '0IA18
MITE (6,*) 'BUB DIA19 *,OIA19
MITE (6,*> 'BUB DIA2B **,01620
MITE (6,**) 'BUB DIA21 >',DIA2I
MITE (6,*) 'BUB 0IA22 "'.DIA22
MITE (6,*) BUB DXAZ3 ',0IA23
MITE <6,) 'BUB DXA26 *',1)1626
MITE (6,) 'BUB 0IA2S '.OlAZS^^
MITE (6,0 'PIPE LENGTH PIPELEN
MITE (6,*l 'PIPE OIA PIPDIA
MITE (6,0 'AREA ' AREA
MITE (6,) 'VOl VOL
MITE (6,) 'CRAV GRAV
MITE (6,) 'CSUBD CSUBD
MITE <6,> 'CAOONS CADMS
MITE (6,) 'CINERA C1MRA
MITE (6,M 'XLEH XL EM
MITE (6,) IWICYC IMCTC
MITE (6,) 'RAXCYC HAXCYC
MITE (6,> 'DX OX
MITE (6,) 'OT *, OT
MITE (6,1 'HCELLS HCELLS
MATERIAL PROPERTIES
C VIS1.005E3
C SURTEN*7.3S0E2
C FOR GLYCERIN
C VIS1.49
C SURTEM0.063
C FOR SUGAR WATER
VIS* 0.2626
SURTEH* 0.073
CVH ABOVE PROPS: VIS(MS/M**2>,SURTE>N/M> AT 20 DEC C
STATE PARAMETERS
PREF l.OES
WATER
RREF I.0E3
GLYCERIN
RREF 1.2560
SUGAR WATER
RREF 1.3330
SOUND l.OES
DROP 1.0/(SOUND*SOUND)
RPREF 1.0
SOLMDP 3.0E2
DRPDP 1.0/(SOUNDP"S0UNDP)
MITE (6,) 'SOUND SOUND
MITE (6,*) 'SOUNDP ', SOtMDP
CONDITIONS AND CELL LOCATIONS
MITE (6,*> 'CELL CENTER POSITIONS'
DO 100 X 1 LSTP2
PRES(X) 1.00E5
X(K) DX*(X1) l.5*DX
MITE (6,> X, X(K)
100 CONTINUE
MITE (6,*) 'JUNCTION POSITIONS*
DO 200 J 1, LSTP2
VEL(J) 0.0
VELSF(J) 0.1
XJ(J) DX"(Jl) l.ft*DX
ACCEL(J) 0.0
XJ(J>
c MITE (6,) J,
200 CONTINUE
c
C XP(1) 0.05
c YPC1) 0.10
XP(1) 0.30
YP(1) 0.05
XP(2) o.st
YP<2) 0.11
XP(3) a 0.15
YP(3) a 0.1S
XP(6) a o.os
YP(6) a 0.06
XP(5) a 0.05
YP(5) a o.u
XP(6) a 0.05
YPC6) a 0.16
XP(7> a O.OS
YP(7) a 0.06
XP(ft) a O.OS
YP(ft) a 0.09
XP(9) a O.OS
YP<9) 0.16
XP(10) 0.05
YP(10) 0.07
XP(ll) 0.05
YP(ll) 0.12
XP(12) 0.05
YP(12) 0.1ft
XPC13) 0.05
YPC13) 0.05
XP(16) 0.05
YP(16) 0.10
XPI1&) 0.05
YPUS) 0.15
XP(16) 0.05
YP(16) 0.06
XPC17) 0.05
YP(17) 0.09
XP(1B) 0.05
YP(lft) 0.16
XPU9) 0.05
YP(19) 0.03
XP<20) 0.05
YP(20) 0.16
XP(21) 0.29
YP(21) 0.05
XP<22) 0.05
YP(22) 0.05
XP(25)
YP(23)
o.os
0.09
XP(26> 0.05
YP<26) 0.11
XP(2S) 0.05
7S5 YP(25) 0.12 MITE(6,783) FORIUTC BUBBLE X AND Y POSITIONS',/,5X, BUB 1',* XP
1 ' YP ')
576 DO 733 I I,MR MITE(6,576)1.XP(X),YP(1> F0RHAT(6X,I3,2XFA.6,2X,FA.6)
733 CONTINUE
DO 400 I 1,
c IF(1BUA(I).EQ.O)GO TO 600 XP(1) 1.5 (FOR OXO.Ol)
CVW SET INITIAL BUBBLE POSITION TO 6TH JUNCTION
C XP(I) 0.02
VELP(I) 0.0
CVW INSERT DIAICTER (CR) OF INITIAL BUBBLE
VOLPtl) PIE"(DIA1)aaS/6.0
VOLPt2) PIE"
VOLP(3) PIE(DIAS)3/6.0
V0LPC4) PIE*(DIAt)**3/60
VOLP(S) PIE(0IAS)3/6.0
V0LP(6) PIEa
V0LP<7) PIE*(DIA7)**S/6.0
VOLP(B) PIE"(0IAA)*a3/6.0
V0LP(9) PIEa(DIA9)aa3/6.0
VOLP(IO) PIEa(DIA10)aaS/6.0
VOLPt11> PIE"(DIA1I>aa3/6.0
VOLP(IZ)* PIEa(DIA12)aaS/6.0
VOLPt13> PIEa(DIA13)flaS/6.0
V0LPU4) PIEa
VOLP(IS)" PIEa(DIA15)aa3/6.0
VOLPt16)* PIEa<0IA16)a"l/6.0
VOLP(17)" PIEa(DIA17)aa3/6.t
VOLP(IA) PIEa(DIAlA)aal/6.l
VOLPt19) PIEa(DIA19)""3/6.0
VOLPt20) PIEa(DIA20)aai/6.0
VOLPt 2D PIEa(DIA2I)aa3/6.0
VOLPt22) PIEatDIA22)aa3/6.l
VOLPt23)* PIE*(DIA23)aa5/6.0
VOLPt24> PIEa(DIA24)aaS/6.9
VOLPt25) PIE*(0IA25)aa3/6.0
CVW
RELTIM1)
RELT1M2) 
RELTIM(S) 
RELTIH(A) 
RELTIH(S) 
RELTIH(6) 
RELT1H(7)
RELTXMS)
RELTIM9)
RELTIM10)
RELTIHtll)
RÂ£LTIM(12)"
PÂ£LTIN(13)
RELT1HC16)*
RELTIM(IS)
RÂ£LTIH<16>
RELTIMt17)
RELTXHtlB)*
RELTIMt19)
.on
.001
.101
0.3
0.36
0.5
2.1
2.0
1.9
2.6
3.3
3.6
3.6
3.7
3.6
6.5
6.A
6.0
5.1
C
C
c
c
c
c
c
c
c
c
c
RELTIMt20)" 6.9
RÂ£LTIM(Z1) 0.001
RELTIMt22)" 5.3
RÂ£LTlM(23)a 6.9
RELTIMt24)" S.l
RELTIMtZS) 5.1
RADUSP(I) (0.75aVOLP(X)/PIE)aaCl.0/3.0)
SCALEP(I) XLEHaOX
DO 300 K > 1. LSTP1
CALL%RACT
300 CONTINUE
NPLOCJ(l) IOINTC CXPam.5DX)/DX ) 1
MPLOCV(I) IDZNTt (XPtl)* DXJ/DX ) 1
600 CONTINUE
MITE <6,0 'RELTIH1 *, RELTXN(l)
MITE <6,) 'RELTIH2 *, RELTIMt2)
MITE (6,") 'RELTIK3
MITE (6.0 RELT1M4 V
MITE (6,0 'RELTXMS V
MITE (6,*) 'RELTIH4 '
MITE (6,) 'RELTIH7 V
MITE (6,0 'RELTINA V
MITE (6,0 'RELTXM9
MITE <6,0 'RELTIH1B"',
MITE (6,) 'RELTIH11"'
MITE (6,") PELTIN 12*V
MITE (6,0 *RELTXN13>1
MITE (6,a> RELTIN16*
MITE (6,0 'RELTIN15"'.
MITE (6,") 'RELTIHX6"',
MITE (6,a> 'RELTIM17'
MITE (6,o) RELTIHIA**.
MITE (6,0 'REITIM19', RELTXN(19)
MITE (6,0 RELTXN2M*, RELTIM(20)
MITE C6,"> RELTIN2I*V
MITE (6,a) 'RELTXH22"*.
MITE (6,0 RELTIMZS',
MITE (6,a) RELTIN26*.
MITE (6,a) 'ICLTXH25*1.
BOUNDARY OATA
VL 0.0
PR 1.0E5
RETURN
EMO
SUBROUTIIC TRIDIA
THIS SUAROUTI* SOLVES THE TRIDIAGONAL PRESSURE EQUATION THAT
RESULTS FROM DC CONTINUOUS PHASE HASS BALANCE. TIC UMNOMN IS
THE PRESSURE XHCREICMT (ICWOLD)
IHCLUDE(TTCONHON)
RELTIH(3)
, RELTXH(6)
RELTIN(S)
RELTIH(6)
RELTZH(7)
RELTXH(A)
RELTIN(9)
RELTimiO)
RElTIH(ll)
RELTINU2)
REITXH(IS)
RELTIN(I6)
RELTIMt15)
RELTIMt16)
RELTIH(17)
RELTXH(IB)
, RELTIM16)
RELTW17)
RELTIMt1A)
RELTIMt19)
REITXM20)
TRIDIAGONAL PRESSURE SOLUTION (A(2) AND C(LST) MOT USED)
KPPC2I C(2)/B(2)
X0Q(2) 0(2)/B(2)
00 1100 R 3, LSTH1
69
(rnrdoiQA (DroxoA cnraiOA
nm ( t p )raaxoA( oo*ao i is* (x' r >rdOXOA >di
O
I o
SS I
no onnnn n
sg5sg5
hOI/ihOo
i5e:fis:
s3 i
Â§
l
sj
5
I
S iA
a o
l*N 
non
a
S
It 1
M9
Xo
jO
X
0
1
9
s
V
H
n
ns
M
I*
u
s
s
WRITE (6,9*20) X, XP(X), VELP(I), RAOUSP(X), VOIP(I),
RHOPA(I), SCALEP(I)
100 K X, LSTP2,*
WRITE (6, 9020) K, VEICKJ, PRES(K), VOI0(K)
 ,PRES(K)PRESH0(K) PRESHOCK), VELSF(K)
n n no non n o n n n n onnnn n
M0N3a/((IX)MX(a>V
HOH9Q/CXO (X)ddX
(IX)ddXH(X)V (X)8 > HON3d
300 CONTINUE
C
MO CONTINUE
C
RETURN
END
SUBROUTIW PNON
C
C THIS SUBROUTINE SOLVES TW PARTICLE LINEAR mwENTUH EQUATIONS FOR
C THE VELOCITIES AN0 THEN UPDATES THE PARTICLE POSITIONS. IT ALSO
C GETS THE INTERFACE TERNS IN THE PARTICLE MNCNTUH EQUATIONS FOR USE
C IN T* CONTINUOUS PHASE HONENTUH EQUATIONS. IT ALSO CALCULATES
C THEPARTICLE DENSITY TERNS CEDED IN THE CONTINUOUS HASS EQUATIONS.
C
IHCLUDE(TTCOfKM)
C
C
DO 100 J 1, LSTP2
XFACEP(J) 0.0
ACCEIP(J) 0.0
GRAVP(J) 0.0
100 CONTINUE
C
DO ZOO k 1, LSTP2
QAEPUJ 0.0
200 CONTINUE
C
DO 700 I 1,
IF
C
C NOICNTUN EQUATION
C
C ADDED NASS COEFFICIENT
ADDHAS CAODMS"RHOAVE
C ADOKAS (RHOPAN(I)M.$hRHONII))VOLPN(I)
C
C PRESSURE GRADIENT SEEN BY PARTICLE
DPDXAV RH0AVE(I)"(G8AV'CItCRA"ACCAVE(I))
C
c IHIIMMIIIHIWMMIIIIIHMIIHIHIHWHIMIHMIlWmHIWHI
C CALCULATE TIC CORRECT VELOCITY FOR THE DRAG TERM
NBUB IKDEXR(X)
IF(HBUI.EQ.O)THEM
FIZ
FI3(X) 0.1
VP 0.0
VS 0.0
ELSE
VP VELPRKX'HBUB)
VS VELSEC(IN8UB)
ENDIF
C ADVANCE PARTICLE VELOCITY
SUM  RHOPANCI>"V01PH(I)"VELPN(I)
* DTftMOPAM(X)"VOlPN(I)"C&AV
 DT*VO10N(I)*OPDXAV
DT" (Fill)"VELAVE(I) FI2
ADDHAS"VELPN(1) 0T"A00MAS"ACCAVE(I)
COEF RHOPAJMI )*V0LPN( I) ADDHAS
OT (FT(I)FX2(U*FnCim DT"FPNALL(I>
VELP(l) SUH/COEF
903 FORHATf IN PNOM, FOR BUBBLE*13>/* FI1,21,3E1Z.S,/,
I * VELAVE,VP,VS',3E12.S)
C
C H* ADVANCE POSITION EQUATION
C
XP)
letoeju) IDINT( (XP/DX ) 1
NPLOCV(X) IDIHTC (XPUJ* DXJ/OX ) 1
C
C GET INTEGRATED EFFECTS OF PARTICLES FOR CONTINUOUS EQUATIONS
C
DO SCO J Z, LSTP1
K J
ACC ( VElPt UVELPNt I) 3/DT
C XFACEP(J) XFACEP(J)
C PERCTJ(J,I)( FI
C AODNAS"(ACCACCAVEU)> )
XFACEP(J) XFACEP(J)
perctj>
FI2(X)
ADDHAS"(ACCACCAVE(D) )
C XFACEP(J) XFACEP(J) PERCTJ(J.DFPVALLU)VLP(I)
TERM PERCTJ(JPI)RHOPAN(X)VOLPN(I)
ACCELP(J) ACCELP(J) TERM"ACC
GRAVP(J) GRAVP(J) TESM"GRAV
500 CONTINUE
C
DO 6## K 1, LSTP1
J I
QADP(K) DAOP(K) VOIDP(KfZ)"ORPOP/RNOPAN(I)
CALL FRACT(XJ(J),XP(I),XJ(J+1),SCALEP(I),PERCT(A,D)
60# CONTINUE
C
70# CONTXttiE
CW INDEX DRA6 COUNTER FOR ADDITIONAL DRAG CALCULATION
I0RA62MDRAG2+1
RETURN
END
SUBROUTINE STATE
THIS SUBROUTINE SETS ALL STATE PROPERTIES WITH PRESSURE AND
TEWERATISE AS TIC INDEPENDENT VARIABLES
INCIUOECTTCONIOH)
DO 100 X I. LSTP2
RHO(K) RREF 0R0P(PRÂ£S(K)PREF)
RHOP(K) RPREF DRPDP(PRÂ£S(K)PREF)
100 CONTINUE
00 300 I l,
XF(IBUB(I).EQ.#)GO TO 300
RHOPA(I) 0.0
DO 200 K 1, ISTP1
RHOPA(I) RHOPA(I) PERCT(X,I)"RHOP(R)
200 CONTINUE
300 CONTINUE
71
RETURN
CKO
SUIROUTIX WAKE(K,RAT101,RATI02,INDEX)
C THIS SUBROUTINE CALCULATES NEW 0RA6 COEFFICIENT (CSUBO) FOR BUBBLES
C THAT ARE INFLUENCED BY TX WARE OF ANOTHER BUBBLE.
C
INCLUDE!TTCOHHON)
C SET AREA RATIO TD ZERO FOR EACH BUBBLE
DO BOS L"1,NP
IF!L.EQ.K.OR.IBUB!L).EO.B) CO TO 603
ARATIO!L)*0.0
ARATI02(L)0.0
603 CONTIHUE
C
C WRITE(6,891)K
891 FORMAT!/,* BUBBLE NUMBER \I3)
IFLC 0
DECAY 28.0
C FOR THE CURRENT BUBBLE (K), CX TO SEE IF IT IN TX WAKE OF ANOTHER
C BUBBLE tJ)
DO 100 Jl.NF
IF! J.EQ.OGO TO 100
IF!IBUB!J).Efl.0)SO TO 100
IWW 0
IPRIM 0
C IF COMPARISON BUBBLE HAS PREVIOUSLY XRGED OR IT K, SKIP OUT
IF!J.EQ.K.OR.IBUB!J).EQ.0) CO TO 100
C CM TO SEE IF BUBBLE COULD BE IN TX WAKE OF ANOTHER BUBBLE, BASED ON
C SOLEY ON ITS VERTICAL (XP) POSITION. IF TIC BUBBLE IS ABOVE ANOTHER
C BUBBLE, DON'T CONSIDER IT
IF!XP!K).GE.XP!J))TXN
CONTIHUE
CO TO 100
EKDIF
CHIIMIIMtMHIHMNIRIIIMIIIII
C IF TX VELOCITY OF TX BUBBLE CAUSING WAKE IS SMALL, TXN TX CURRENT
C BUBBLE IS NOT EFFECTED BY TX WAKE, SO SKIP OUT
IF!VELPN!J).LT.1.0ES)SO TO 100
CHIIHHMUIIIIIHmillllllll
C CALC TX LENGTH OF SECONDARY WAKE FOR BUBBLE J
C IF TX DISTANCE BTW BUB K AND BUi J > SECOND WAKE LENGTH SKIP OUT
C WAKE LENGTH BASED ON STMMILLEft
PERCENT 0.01
SOLU ( !0.)2*PCRCEXT)*0S0RT( (0.12PERCENT>2 
1 4.*0.tlPÂ£RCENT"!PERCENTH.2~l.0) ) )/ !2.0.01PEftCEMT)
. WAKELS SOLU RAOUSP(J)
BUBDIS XP(J) XP(K)
IF(BUB0IS.6T.WAKELS) TXN
C WRITE(6,70S)K.J, WAKELS,BU101S
70S FORMAT!* BUBBLE*,IS.* IS HOT EFFECTED BY BUBBLE,13,' WAKE BASED
1 ON TXIR VERTICAL DISPLACEXNT i*, /,
2 SECONDARY WAKE LENGTH"' ,E12.S BUBDIS *,E12.S)
GO TO 100
EKDIF
CWMMWWW
C NOW CHECK IF ITS IN TX WAKE DX TO HORIZONTAL POSITION
C CW TO SEE IF TX BUBBLE TOUCXS OR OVERLAPS TX WAKE OF BUBBLE J
C AT XP !CONTINUE IF TXY DON'T TOUCH)
C
C CHECK TO SEE IF BUBBLE OVERLAPS SECONDARY WAKE FIRST:
CALL SECWAKE!J,K,VELSW,WAKRDS)
YY11 YP(K) RADUSP(K)
IF!YYI1.LT.O.O)YY11O.I
YY5S YP(K) RADUSP(K)
IF!YY3S.6T.PIPDIA)YYS3PIPDIA
YY5S YP(J) WAKRDS
IF(YY5S.LT.0.Q)YY550.B
YY66 YP!J) WAKRDS
IF(YY66.CT.P1P0IA)YY66PIP0XA
IF! YY11. ST. YY66. OR. YYSS. LT. YYS5 ) TXN
APROJB > PIE RADUSP(K) mi 2
GO TO 100
EKDIF
C
C IF BUB DOES OVERLAP SEC WAKE (ABOVE), CXCX TO SEE IF BUB IS lELON
C TX ESTIMATED WIDTH OF TX PRIMARY WAKE
XWS 7.1
XW XVS RADUSP(J)
IF!BUBDIS.6T.XW)G0 TO 7777
C IF BUB CAN OVERLAP PRIM WAKE AND SEC WAKE, CALC PRIM WAKE RADIUS
ALPHA DATAN!RADUSP(J)/XW>
TR1LEN XW BUBDIS
WAKRDP TRILEN DTAN!ALPHA)
CWUfMMlWWHUMMimMIIIW
C CM! PLACEMENT OF BUBBLE TO TX WAKE AT XP(K), CM TO SEE IF BUBBLE
C OVERLAPS TX WAKE ON TX RIGHT OR TX LEFT OF TX WAKE OR IF TX
C BUBBLE IS ENTIRELY WITHIN TX HAKE
C BUBBLE IS LEFT OF WAKE
C SKIP OUT IF ENTIRELY IN WAKE OR PROJ WAKE AREA W/IN BUBBLE APROJB
APROJB PIE RADUSP(K) mi 2
C WRITE!6,243)K,J.WAKRDP,WAKRDS
243 FORMAT!* FOR BUB*,IS,' 1XLUEHCE BY WAKE* ,13, *RADP/RABS* ,2E12.5)
YYZ YP!J) WAKRDP
IF!YY2.LT.O.O)YY2O.I
YY4 YP!J) WAKRDP
IF!YY4.6T.PIPDIA)YY4PIP0IA
YY1 YP!K) RADUSPU)
IF!YY1.LT.O.O)YY1"O.Q
YY3 YP(K) RADUSP!K)
IF!YY3.CT.PIPDIA)YY5*PIPDIA
YY5 YP(J) WAKRDS
IF!YYS.LT.O.O)YYSat.l
YY6 YP!J) WAKRDS
!F(m.9T.PtPDU)VV**PXPDXA
AWAKES PIE WAKRDSM2
AWAXE1 PIE WAKRDPZ
C
C DEFIX BUBBLE POSITION IN WAKE (IWW)
C 1 BUBBLE ENTIRELY IN SECONDARY WAKE AND COIPIETELY OVERLAPS PRIMAR
C Y WAKE
C 2 BUBBLE ENTIRELY IN PRIMARY WAKE
C 3 BUBBLE ENTIRELY IN SECONDARY WAKE BUT DOES NOT OVERLAP PRIMARY
C 4 BUBBLE INTERSECTS SECONDARY WAKE, BUT NOT PRIMARY
C 5 BUBBLE INTERSECTS BOTH PRIMARY AND SECONDARY WAKE
C 6 BUBBLE INTERSECTS SECONDARY AND OVERLAPS PRIMARY WAKE
C 7 BUBBLE ENTIRELY IN SECONDARY AND INTERSECTS PRIMARY
C S SECONDARY WAKE "ENCLOSED** IN BUBBLE RADIUS
C
IF!YYS.LE.YY1.AND.YY1.LT.YY2.AND.YY3.GT.YY4.AND.YY3.LE.YY6)IWW1
72
IFCYYZ.LE.YY1.AND.YY3.LE.YY4)IWN2
IF(YY5.LE, YY1. AND.YY3.LE. YYZ) IW*S
IF(YY4.lE.YYl.AMD.YY3.LE.YY6)iyU*S
IFCYV5.LT.YYS.AMD.YYl.lT.YYS.AND.YY3.lE.YYZ)IWW4
XFCYY4.LÂ£.m.AND.m.LT.YY6.AHD.YY3.6T.YY6)IWN*4
tFCYY3.GT.YYZ.AMD.YY3.LT.YV4.AND.YYl.lT.YYS)Ilrtf5
IFCVY1.GT. YY2.AMO.YY1.LT.YY4.AND.YY3.6T.YY6JIVW*5
IF C VY1.GT.YY5.AMO.YY1.LE.YY2.AMO.YY3.GT.YY6)IWM*6
IFCYY3.GE.YY4.AND.YY3.LT.YY6.AMD.YY1,LT.YY5)IWU6
IF(YY1.CE.YY5.AM0.YY1.LT.YY2.AMD.YY3.CT.YY2.AMD.YY3.LT.YYA)IW7
IF(YYl.CT.YY2.AN0.YYl.LT.YYA.AMD.YY3.GT.YYA.AM0.YY3.LE.YY6)lMM7
IF
IF(NCYC.EQ.S17.0R.KCYC.EQ.518)THEN
C WITE(6,152)NCYC
152 fORNATC' FOR CYCLE**,13)
EN02F
C SET ARU RATIOS RASED ON BUBBLE POSITION IN UC UAKE(S)
IFdWW.EQ.DTHEN
ARATIOU) AVAXE1/APG0J1
ARATI02CJ)
ELSE1F(1W.EQ.2)TN
ARATIOCJ) l.l
ARATI02CJ) 0.1
ELSE2F(IM.EQ.3)T(CN
ARATIOU) 0.0
ARATIOZU) 1.0
ELSEIFCIMW.EQ.AJTICN
ARATIOU) > 0.0
CALL XNTEft(K,J,YY5,YY6.AIKT2)
ARATIOZU) A1MT2/APR0JB
ELSE1F(XM.Â£Q.5>TICN
CALL XMTER(K,J,YY2,YYA,AINT1)
CALL INTER(K,J,YY5,YY6,AINT2>
ARATIO(J) AINT1/APROJB
ARATIOZU) QABSCAINTZ AIMTD/APWUB
ELSEIFCIW.EQ.UTICN
CALL INTER
ARATIOU) > AHAKE1/APRQJB
ARATX02CJ) DA1$(AINT2~AMAXE1)/APR0JB
ELSElF(IUM.Efl.7)THEH
CALL INTER(K,JYY2,YY4,AIHT1)
ARATIO(J) AXNT1/APROJB
ARATI02U) DABS(APROJi A1XTD/APGQJB
ELSEXFCIW.EQ.8)TICH
ARATIOCJ) AHAKE1 / APRQJB
ARATI02CJ) DABSCAUAKC2 AMAKE1)/APROJB
ELSE
VRITEC6,429)IIM
420 FORMAT!* TTTT????? IM DOES NOT 18, IW*, 13,*?????????????*)
ENDIF
IFL6 1
GO TO 108
C INTERSECTION WITH SECONDARY HAKE ONLY
7777 CONTINUE
C CHECK TO SEE IF BUBBLE ENTIRELY HITHXN HAKE OR JUST INTERSECTS
YP1 YP(K) RAOUSP(K)
ifcypi.lt.o.o)YPio.e
YP3 YPUO RAOUSP(K)
IF( YP3. GT .PIPDIA)VP3*IPDU
YP5 YPCJ) UAXRDS
IF(YPS.IT.0.9)YP5*I.I
YP6 YPCJ) MAKROS
IF(YP6.GT.PIPDIA)YP4PIPDIA
C
IF(YP(K).Â£g.YP(J))TtCN
IF ( YP1. GE. YPS. AND. YP3. LE. VP8 )T)CN
ARATIO(J) l.l
ARATZ02CJ) 1.8
ELSE
ARATJO(J) l.l
ARATIOZU) (PIENAKRDSm2)/(PXERAOUSP(K)mi2)
ENDIF
ELSElFCYPl.SE.YP5.AND.YP5.LE.YP8)THOI
ARATIO(J) *1.1
ARATI02CJ) 1.0
ELSElFCYPl.LT. YP5. AND. YP3.GT.YP8)THQI
ARATIOC J)M.I
ARATI02CJ)* PIE*ttAKRDS*"2/(PlE*ftA0U5P(K)*"2)
ELSE
C CALL XHTER(K,J,MKRBS,AINT2)
CALL XNTERCK,J,YP5,YP6,AIMT2)
ARATIOCJ) 0.0
AAR* PIE CRADUSP(K))2
ARATI02CJ) AINTZ/AAR
ENDIF
IFLG 1
C MWIWIBMMWMMMMMMBWIiWHWIWmWWWmWWmiWHWHM
110 CONTINUE
C
C FIND LARGEST ARATIO AND CORRESPONDING BUBBLE, THIS IS TK HAKE THAT
C HIU IWLUENCE BUBBLE K (IFLC1 AHS THAT TIC BUBBLE IS IICLUENCED
C BY THE MAKE OF AMQTtCR BUBBLE)
IFdFLG.EO.DTICN
BIG 0.9
INDEX I
C MRITEC6,S01 )K
SOI FORNATC/, FOR BUBBLE*,13,/)
DO 300 H 1,1*
IF(IBUB(H).EQ.O)GO TO 301
IF(N.EO.K)GO TO 300
8D(N) XPCN) XPCK)
SSU8NN) ARATIOCH) ARATI02CH)
IFCBDOO.LE.O.OIGO TO 301
IF(SSUW(N).6T.I6)TTCN
BIG SStfWCH)
INDEX N
ENDIF
301 CONTINUE
IFCB1G.LT.l.DGO TO 351
XL IT 1.0E5
DO 3ZS N 1,M>
IF(IBUB(N).EO.O)GO TO 3ZS
IFCN.EQ.K)GO TO 32S
DIS XPCN) XPCK)
ADD ARATIOCH) ARATX02CN)
C URITE(6,395)MADD,DISXLIT
395 FORNATCBN*',13,* ADD*',E15.9,/,' DIS*.E15.9,* XLIT**.Â£12.5)
XFCDXS.LT.O.DGO TO 3Z5
IFCADD.LT.0.99999)60 TO S2S
73
C HRITE(6,396)N.ADD.DIS,XIIT
3H FORMAT! AMML.* ADD* ,Â£14.6,/,* DIS"*,E12.S,* XLIT* .E12.S)
IF(DIS.Lf.XLIT)ThCN
XLIT DXS
INDEX N
C HRITE(6,830)XL1T,DIS,IMDEX,N
B39 FORMAT! IN IF XLIT* ,E12.S, DIS* ,E12.5,* IWEX/)** ,213)
ENDIF
325 CONTINUE
350 RATI01 ARAT10(INDEX)
RATIOS ARATI02!INOEX)
C *
ELSE
RATI01 9.9
RATIOS I.*
INDEX 9
C WRITE!6,713)K
ENDIF
717 FORMAT!/,* BUBBLE*,13, IS INFLUENCED BY BUBBLE*,13,/,
1 * WITH RATI01*'.E1Z.5,1 RATI02 \E12.S,* RADUSP*,Â£12.5,/)
713 FORMAT!/,' BUUIE',13,* IS NOT I1FLUENCE0 BY ANY HAKE',/)
SOO CONTINUE
600 CONTINUE
RETURN
END
SUBROUTZfC INTER!K,J,PS,P6,AINT)
C SUBROUTINE TO CALCULATE AREA OF INTERSECTION OF BUBBLE AW HAKE
INCLUDE!TTCOIMN)
AL 1.0
AR  9.0
THETAU l.l
THETAl t.l
TtCTAR 9.9
PIEH PIE/2.9
C SET CORRECT BUBBLE OR HAKE TO RI6KT/LEFT
IF!YP!K).GT.YP!J)) TtCN
RR RADUSP!K)
C RL  HAKRAS
RL Pi YP(J)
NT RL
ELSE
C RR  WAX SAD
RR YP(J) PS
HT RR
RL RADUSP(K)
ENOIF
C IF TTC BUB VELOCITY IS ZERO, TfERE IS NO HAKE, SO AINT9 9 RETURN
C XF(HAXRAD.LT.1.0Â£5)TtCN
IF!HT.LT.1.0C5)T>CN
AINT 9.9
RETURN
ENOIF
C CALCULATE SEPARATION DISTANCE
DISOABS!YP(K)YP!J))
C FIND LAR6EST SIDE OF TRIANGLE AW THEN FIND ANCLES
XMAX DIS
IFtRL.6T.XMAX) TWN
XMAX RL
NO* 2
ELSEIFtRR.6T.XHAX) TXM
XMAX RR
NCtft 3
ELSE
CONTINUE
ENDIF
IFtNCM.EO.imCN
THETADDACOS! !DIS*a2Rl2RRM2) / ( (2.9)Wl"Rft ) )
TtCTARDASINI RLOSIN(TICTAU)/DIS )
THETAL*DA$INtRR"DSIN!1TCTAU)/DXS)
ELSEIFt NCHC. EO. 2)T)EN
TtCTARDACOS! (RL2RRm2DIS2) / t !2.l)RR"0XS) )
TfCTALDASINt RRDSXNt TNETARJ/RL )
TCTAI>*PIET)CTALT)eTAfi
ELSEXF!NCH(.EQ.3)TfCN
T>CTALaDACOSf !RRH2RLm2DXSH2) / t2.0)*RL*DIS) )
T)Â£TAR"DASIN(RL"DSINtTtCTAL)/RR)
TJCTAUPIETHETALTHETA*
ELSE
CONTINUE
ENDIF
C SET CASE
ICASE 1
XF!T)CTAR.6T.PIEH)ICASE 2
IF!DCTAL.6T.PIEH)XCASE 3
C CALCULATE INTERSECTION AREAS
IF(ICASE.E0.1)TtCN
AR 0.5"RRM2t2.TTCTAR DSXMtZ.OaTHETAR))
AL O.SRLM2a(2.0MTHETAL BSXN!2.0aTTAl>>
ELSEIF! I CASE. EQ. 2)T)CM
AR tPIEttRR**2)(9.5*RR**2<*(2.9*tPIET)CTAR)
1 OSIN!2.!(PIETHETAR))))
AL I.SHRL*a2(2.0TtCTAL DSIN!Z.0aT>CTAl)>
ELSEIFt ICASE.EQ.3)T)CN
AR 9.S*RR**2a(2.0aTtCTAR DSIN!2.9aTHETAR)>
AL tPIEaRLa2)!9.5*RLaa2B(2.9a!PIET>CTAL)
1 DSINC 2. Oat PIET)CTAl))) )
ELSE
C HRITE!6,761)
7A1 FORMAT!* YOU MESSED UP ICASE DOESNT EQUAL 1,2 OR 3*)
ENDIF
IF!1TCTAU.EQ.I.9.0R.T>ETAR.EO.I.9.0R.THETAL.EO. 0.0 )HRITE( 6,516)
516 FORMAT! aBMMWHMBMMwi A THETA EQUALS
C CALCULATE AREA RATIO
AINT AR AL
RETURN
END
SUBROUT1W fCRSE
C SUBROUTINE TO CHECK IF BUBBLES HERCE, OCCX FOR SLUG FION, AW CALC
C WH (CRGED BUBBLE PARAMETERS
INCLUDE!TTCOMHON)
IEND 9
DO 101 I"1,W
IF!IBUBCD.EQ.DGO TO 109
C OCCX TO HAKE SURE BUBBLE IS RELEASED
IF(RELTXH!I).6T.TXrC)GO TO 100
C OCCX TO SEE IF BUBBLE OUT OF TOP OF PIPE
IF! (XPtI)RADUSPtD).GE.XJ(lSTFl) )THEN
WRITE!2,234)1,TIME
74
234 FORMAT(* BUBBLE*,13, OUT OF TOP OF PIPE AT TIIC (S)*,FB.4)
IBUS!I) 0
60 TO 101
EKZJZF
e oca for siue flow based on max eouiv oiaicter stop
IF(RADUSPII) .GE.REMAXITICN
I END 1
WRITE!6,409)HCYC,I,YP!I),XP!t),RAOUSP!I).ISHAPE!I) 
WRITE!2,409)NCYC,I,YP!I).XP!I),BADU$P(I),ISHAPECI)
409 FORMAT! IIHIIIHIHHimimillHIIIIIIIHIIIHIIIHItlir ,
1 /,' AT NCYC*,15,* BUBBLE* ,13, REACHED SLUG FLOW SIZE* *
2 /,* TP*>E1?.S* XP* ,E12.S* R*,Â£12.S, SHAPE',13,
3 /, tillIlf49llllt*i0tllll14499404MilIttlIII44190IIIIt' )
60 TO 543
ENUIF
C OCa TO SEE IF BUBBLE OVERLAPS SIDE OF PIPE
IF (I SHAPE (I). Â£ Q .1) TtCH
RIGHTYP(I) RADUSP(X)
XLEFTYP(I) RAOUSP(I)
ELSEIFCISHAPE!I).Â£Q.2)TtCN
AZZ 0.629960*01 2.D*tO RADUSP(I)
A2Z* A22/2.D*IB
SI6HT YP(I) A2Z
XLEFT YPCI) A2Z
ELSEXF (ISHAPE! I). EG. 35TWN
AZZ* (!2.D*00*RADUSP!I))"*S/!!S.D*40/4,020*04>
1 (1.0*00/4.02D*0t*3))) m (1.0*40/3.0*00)
422* 422/2.0*00
RIGHT YP(X) 422
XLEFT YP(I) A22
ELSE
tMTE!6,59l)NCYC,l
590 FORMAT!* ?T??TT?? ISHAPE 1C 1,2/3 AT HCYC*,13,* BUt',I3,*mT*)
CONTINUE
EWIF
IF(RIGHT.6T.PIPDXA)YP(X)YP(X)!RICHT'PIPDIA)
IF(XLEFT.lT.i.4)YPCJ)YP!I)*0ABS!XL6FT)
154 CONTINUE
C
DO 114
IF(IBUB(J).EO.O)60 TO III
IF( J.EQ.DGO TO 114
C OCa TO HAKE SURE COWARXSON BUBBLE IS RELEASED
IFCRELTIH! J).CT.TXIC)GO TO lit
C
asm m RADUSP(J) RADUSP(I)
TERM1 DABS! XP(J)XP(IJ )
TERM2 DABS! YP!J)YP!I> )
RCALC DSBRT! TERH1*2 TERM22 )
IF!RCALC.CT.HSUH)GO TO 111
C CALCULATE ICW VELOCITY, RADIUS ISXZE) AND POSITION OF ICW BUBBLE
741 XHASS1 RHOPANIX) VOLPN(X)
XNASS2 RHOPAN(J) VOLPN(J)
VEUCH (XMASS1 VELP(I) XMASS2 VELP!J))/!XHASS1 XMASS2)
VOOCV (XHASSr XMASS2)/!(RH0PAN!I)*RH0PAN!J))/2.D*0IJ
RADWW ((V0LPCIfa3.l)/(4,l*PXE))mi.33BS
C OCa FOR SLUG FLOW
IF (RA9CV.6T. REMAX )T>CM
IEKD 1
VRlTEl4,4ll)tCYC,J,I,RAIBCW
WRITE! 2,441 )NCVC,J, I,RAUCH
440 FORMAT!* IM4IM4ll4t4Mtlt944MMM4MMt4IMMM4IMtM*,
1 /,* AT NCVCM5,' BUBBLE,13,* MERGED INTO BUBBLE*,13,
2 /,' ICV RADIUS *',E12.S,' ONSET OF SLUG FLOW*,
3 /,* litOilOOttllOOOOtOOOOOOOtttiftOltllltOltlltlOOIOOft')
GO TO 54S
EWIF
C
YTC** YP(n*VCLPN(D* YP!JlVOLPN!J))/ !V0LPN!I)*V01PN(JJJ
RHXCW (RHOPAN(I) RMOPAN(J)) /2.0
C REASSIGN TO VALUES TO TtC ICRGED BUBBLE AND TAKE QIC OUT OF SYSTEM
C REASSIGN TO ICW VALUES AND T>CY WILL IE ASSIGNED OLD VALUES IN
C TIC MAIN PROGRAM
IBUB(J) 0
VELP(X) VELtCV
YOLPNdl* VOUCH
RADUSPU) RAIWCV
YP(IJ YKW
RHOPANIX)* RHOCW
WRXTE(6,5IO)J,I,VOLPN(I),RHOPAN(I),RADUSP(I).VELP(1),YP(X),
1 NCYC
FORMAT!* BUBBLE*,13,* ICRGED INTO BUBBLE M3,* ICW STUFF:',/,
VOL*',E12.5,' RHO *,E12.5,' RAO**,E12.S,* VELP*',
E12.S,' YP',E12.5,/,' ICRGE AT MCYC',I5)
WRXTE(2,499)J,Z,V0LP(I),RH0PA(1),RADUSP(I)VELP(Z),VP(X)
NCYC
BUBBLE',13,' ICRGED INTO BUBBLE ',13,* ICW STUFF:*,/,
VOL"',E12.5,' RHO* *,E12.5,' RAD*,E12.S,/, VELP*',
E12.S#' VP**,E1Z.S* ICRGE AT NCYC*.15)
114 CONTINUE
104 CONTINUE
14 CONTIMC
IGRAPHIGRAPH*1
XFINCYC.E0.21GO TO 543
IFCIGRAPM.EQ.O.OR.XG8APH.E0.5B)TTCN
IGRAPH*#
MlTE!2,344>TlfC
FORMAT!' FOR T *',F6.4,( SECONDS:'//,
SW.* HP *,* YP *,* VELOCITY *,
RADIUS *1
DO 44S N 1,*
IF!IBUB!H).EQ.I)GO TO 445
WRITE!2,345)H,ISH4PC(H),XP(H),VP!H),VELP(M),RADUSP(N)
. FORMAT!IX,13,1X,I32X,FB.4,SX,FB4,SX,E12.S,5X,Â£12.5)
CONTINUE
IF(XEMD.EQ.1)ST0P
EWIF
RETUM
END
SUBR0UT1W SECNAKC(I,K,VELS,WAKRDS)
C SUBROUTIIC TO CALCULATE SECONDARY WAKE RADIUS AND VELOCITY
INCLUDE!TTCONHON)
DECAY 26.0
C BUBBLE SEPARATION DISTANCE
OEITAX XP(I) XP(K)
C PRIMARY WAKE LENGTH !XW) AND WINARY WAKE RADIUS (RAOPRI)
C RAOPRI RADUSP(Z)BDEXP(DCCAYBDELTAXT
XWS 7.0
PRIKX XWS RADUSP(I)
C DONT CALC PRIMARY WAKE RADIUS IF BUBBLE IS BELOW PRIMARY HAKE
541
1
2
499 FORMAT!'
1
2
543
344
1
2
345
449
75
IF(PRIHX.LT.DLTAX)GO TO SI*
ALPHA DATAN(RADUSP(I)/PRINX>
TRXLEM PRIHX DELTAX
RAOPRX TRIlEM DTAN(ALPHA)
C PfilHX (1.0/DECAVI a DL06(DSQRTC0.01aRADUSPmaa2>/RADUSP(X
510 CONTINUE
RHOREF RHOPANCI)
IF(VELPN( I I.LT. 1.0E*5)T>CN
CD 24.0
60 TO 10
ENOIF
CO (1>0 RHOREF/PHQAVE(I)) DABS(CRAV) a 6.0 a RADUSP(I) /
1 (VELPHCI)2 3.0)
IFCCD.GT.24.0)CD*24.I
10 CONTINUE
C
C CALC SECONDARY HAKE VELOCITY BASED ON STUNILLER
C UIXJ1 1.0 / ( 0.2 (t.l2B0ELTAX/RADUSP(D) 1.01
C 1 (DELTAX/RADUSPCI))aa2)
C UD1* UDUla VELPNCI)
C UDlPRlHXl./( 0.2 (l.l2aPRXHX/RADUSP(X)) 0.11 a
C 1 CPRINX/RAflUSP(I)>aa2>
C UDPRINX1* UDIPRIHX a VELPN(I)
C CALC SECONDARY HAKE VELOCITY BASED ON CRABTREE
UDU1 (6.2 RADUSP(X / DEL TAX
UD1* UDUla VELPN(X)
UDIPRIHX* (4.2 a RADUSP(I)> / PRINK
UDPRIHX1* UDIPRIHX a VELPN(X)
C CALC SECONDARY HAKE VELOCITY BASED ON BHA6A
C (DON'T DIVIDE BY 2ERO ON FIRST TI*STEP)
IFCVELPN(I).LE. 1.0Â£4>THEM
UD2*0.I
UDPRINX2*0.I
CO TO 20
ENDXF
BVAL"(DAlS(GRAV)a(2.HNADUSP(X))"a2 RHOAVE(X))/ (VEIPN(I)aVIS)
EVAL 1.1/0.7
C URITE( 6,769 )BVAL,EVAL,DELTAX PRINX,RADUSPCX)
C 769 FORHATC* IN SCCVAKE: B\E12.5,' EVAL*'E12.S,/,SX,' DCLX*',E12.5,
C 1 PRINX',E12.S' RADP*,E12.5>
C EQUATION IS BASED ON DISTANCE FRON REAR STAGNATION POINT, PRIKC, TO
C T)C NOSE OF TIC FOLLOWING BUBBLE
DIFFER (DELTAX RADUSPCX)) PRIKC
IF(DIFFER.LT.O.O)DXFFER 0.0
C UDU BVAL/C24.a(CBVAL/24.)aa*.7*CDXFFER/C2.aRADUSPCI>)> )aaCVAL>
UDU2 BVAL/(24.a
1 ((BVAl/24.0)aa.7 (DIFFER/C2.RADUSPCX)) >aa.7>EVALI
UD2* UDU2a VELPNCI)
UDPRINX 2* VELPNCI)
20 CONTINUE
C UD (UDlUD2)/2.0
C UDPRINX (UDRRXHX1 UDPRINX2)/2.0
UD UD2
UDPRINX USPRINX2
VELS UD
VELC UD
IFCVELC.GT. VELPNCI) )VELC* VELPNCI)
IFCUD.LT.l.0E5)THEN
HAKRDS *0.1
VELC B.B
cc
CO TO 659
ENDXF
C CALC SECONDARY HAKE RADIUS
IF< PRINX. 6T. DEITAX)T)CN
RR DSORTC O.S a CD a VELPNCI) / UDPRINX) RADUSP(I)
R2 DABSCRR RADUSP(D)
C IFCRR.LT. RADUSPCI>>URITE(6,659>RR,RADUSPCI)
659 FORHATC 111111 SEC NAK RAD IS LESS THAN BUB RAD* ,2E12.5)
THETA DATAN CPR1KC/R2)
HAKRDS CDELTAX / DTANCTICTA)) RADUSPCX)
C XF(MFL6.EQ.0)URITE(6,496)KX
496 FORHATC* aaa* BUBBLE*,IS,' ENTERED PRXN HAKE OF BUB*.13)
WTL6 1
ELSE
HURDS DSORTC0.5 CD a VELPN(X) / UD ) a RADUSPCX)
ENDIF <
659 CONTINUE
C CALCULATE TTC AVERAfiE VELOCITY SEEN BY TTC BUBBLE
C
IFCYPCK).EO.YP(X))THEN
C
IF(PRIHX.CT.DELTAX)TTCN
RAD RADUSPCX) RADFRI
IF ( HU RDS. L E. RADPRX. OR. RAD. LE 0.0 ) T)CM
VELS 0.0
60 TO 571
ENDIF
IFCRADUSPCK).CT.HAKRDS)RAD HURDS
VELR VELC DEXPC CRAD/VURDSJa2)
VELS > (tÂ£LRVELC)/2.l
ELSE
RAD RADUSPCK)
IFtRAD.BE.HAKRDS)RAD HURDS
VELR VELC a DEXPC CRAD/WAKRDSa2 )
VELS CVELR VEL0/2.0
ENDXF
571 CONTINUE
ELSE
CALL SECPOSCYPCK),YPtX),HAKRDS,VELC,DELTAX,PRINX,RADUSP(K),VELS,
1PIPDIA)
ENDIF
RETURN
END
SUM0UT1IC PRMVf XT XL VWI RT iRL VW)
C SUBROUTIIC TO CALC PRIMARY HUE VELOCITY AT BUBBLES NOSE
C PRIN HAKE LENGTH FACTOR
XVS 7.B
C NOSE POSTION OF TRAILING BUBBLE
XTN XT RT
C HUE LENGTH
UULEN XNS a RL
C HUE TAIL POSITION
HUPOS XL UULEN
6 VORTEX POSITION
XVC (0.6 UULEN) UAXPOS
C
IFCXTN.LE.XVOTHEN
XX (XVC XTM)/UULEN
IF(XX.GT.0.61H0)XX 0.6D0
76
VW VKX> (DSflRTC (l.OEKt (XX/0.6D0)2) 0.SD*02> 1.0D*0>
ELSE
XX (XTN XVO/WAKLEN
IF(XX.6T.0.41H00)XX0.*0*19
VW VKX* (DSQRTC (l.QDM (XX/0.*D*0>*2) 0.5D*02> 1.0D+0)
EKDXF
RETURN
END
SUBROUT DC SECPOS(YIYJ,URVELC,DELTAX,PRINX,RI,VELB,PIPDIA)
C SUBROUTINE TO CALC VELOCITY OF SEC WAKE FOR BUBBLE IN VEL PROFILE
IMPLICIT DOUBLE PRECISION (AH,0Z)
IF (DC L T AX. L T. PR I MX ) TICN
VELB (VELC VELC*DEXP(I.#0*00)) / 2.0
EKDIF
C
D DABS (VI 
POSL YI RX
POSR YI RI
WAKL YJ UR
WAKR YJ NR
URL UR
YJ)
II
IR
C HAKE SURE SEC HAKE DOES HOT OVERLAP PIPE BOUNDARY
XF(WAKL.LT.I.O)THEN
WAXL t.a
XL 1
EHDIF
IF(UAKR.6T.PIPDIA)TtCH
WAKR PIPDIA
IR *1
EHDIF
C
C IF BU1BLE CENTER IN SECUAK THEN VELOCITY IS WHAT BUB CENTER SEES
XF(YX.QT.UAKL.AM)YI.LT.WAKR)THEN
IFCYI.LT.YJ.AND.IL.EQ.DURL YJ
IF( YI.6T.YJ. AND *IR.EQ. DURR PIPDIA YJ
UR DHIN1(URL,WRR)
VELB VELC 0EXP(
ELSE
IF ( POSL. LE YJ. OR. POSR. GE YJ ) TtEN
VELB (mC VELC"DEXP(1.0DM) )/2.l
ELSE
R DABSfORI)
IF(YX.LT.YJ.AHD.IL.EQ.1)URL YJ
IFCYI .6T.YJ.AND. XR.EO.DURR PIPDIA YJ
UR DHINl(VRL.WRR)
VELB VELC DEXP((R/VR)a2)
ENDIF
EHDIF
RETURN
DO
77

PAGE 1
THE NUMERICAL SIMULATION OF BUBBLE COALESCENCE by Valerie Ellen Walker B.S., Colorado State University, 1984 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of Mechanical Engineering 1991
PAGE 2
This thesis for the Master of Science degree by Valerie Ellen Walker has been approved for the .Department of Mechanical Engineering by JOhTi A. J. Kenneth E. Ortega I r;? j Date
PAGE 3
Walker, Valerie Ellen (M.S., Mechanical Engineering) The Numerical Simulation of Bubble Coalescence Thesis directed by Professor John A. Trapp A computer model was developed to accurately simulate the process of bubble coalescence in a bubbly flow regime. The numerical model included proper drag relationships for the various bubble shapes and regimes occurring in bubbly flow. A model of the wake structure behind spherically capped bubbles was developed which included models for both primary and secondary wake structures. The numerical model was incorporated into an existing one dimensional, mechanical, two phase flow program. The revised program having the capability of tracking individual bubbles throughout the flow area, with the ability to simulate bubbly flow through the onset of slug flow. Calculated values of bubble terminal velocities in water compared well with measured data. Simulations of bubble coalescence also compared well with measured bubble positions during coalescence and coalescence times. The comparison were made for sugar water and glycerin coalescences, yielding favorable results from close distances (a few bubble radii) to separation distances up to 80 em. The form and content of the abstract are approved. I recommend its publication. Signed c;M /1. John A. Trapp
PAGE 4
CONTENTS Figures . . . . . . vi Acknowledgements . . vii Nomenclature viii CHAPTER 1. INTRODUCTION . . 1 Background 1 Scope . 3 2. PROGRAM STRUCTURE 4 Govern1ng Equations . . 4 Program Scheme 7 3 NEW CODE FEATURES . . 10 Bubble Hydrodynamics . . . 10 Drag Coefficients and Transition Criteria 16 Spherical Regime 16 Ellipsoid Regime 20 Spherical Cap Regime 23 Bubble Geometry . 23 Bubble Coalescence 25 Wake Structure . . 27 Primary Wake 27 Secondary Wake . . 31 Drag on a Bubble in a Wake . . 35 Bubble Merging . . 37
PAGE 5
4. RESULTS AND DISCUSSIONS Terminal Velocity Coalescence Times Ori gina 1 Mode 1 Final Model .. Flow Visualization Discussion REFERENCES APPENDIX A. Code Listing ... v 39 39 42 43 45 55 58 61 62 63
PAGE 6
Figure 2.1 Program flow chart 3.1. Bubble shapes 3.2. Bubble regimes FIGURES 3.3. Reynolds number vs. drag coefficient 3.4. Wake structure behind a spherical cap bubble and Littman) . . 3.5. Primary wake shapes 3.6. Primary wake velocity 3.7. Wake structure behind a spherical cap bubble 8 11 13 17 (Bessler 28 29 32 (code model) . . . . . . . ...... 36 4.1. Terminal velocity of bubbles in water (Peebles and Garber 1953) ............... 40 4.2 Terminal velocity of bubbles in water (Haberman and Morton 1953) ..... . 41 4.3. Bubble positions during coalescence original model 44 4.4. Coalescence of same size bubbles original model 46 4.5 . Coalescence of a larger trailing bubble original model 47 4.6. Bubble positions during coalescence final model 49 4.7. Coalescence of same size bubbles final model . 51 4.8. Coalescence of a larger trailing bubble final model 52 4.9. Glycerin coalescence . 54 4.10. Water coalescence 56 4.11. Numerical simulation of bubble coalescence in a pipe 57
PAGE 7
ACKNOWLEDGMENTS I wish to thank John Trapp for the thesis topic and for the use of the existing computer. program. I would also like to thank Gary Young for his support and confidence throughout the project.
PAGE 8
Variables p density V Volume u velocity U terminal velocity t time NOMENCLATURE x horizontal distance or position y vertical distance or position g gravity fw wall friction coefficients f0 drag component f added mass component a void fraction c speed of sound r radius (equivalent of a sphere) D total drag C drag coefficient d0 equivalent diameter P pressure B gd2p/Up. A area a surface tension p. dynamic viscosity o(s) percent of bubble in grid cell LHS left hand side RHS right hand side a width of bubble b height of bubble L' distance from end of primary wake to following bubble u velocity with u,v,w as the i,j,k Subscripts b discrete (bubble) pw primary wake we wake center line sw secondary wake f free stream p primary s secondary Superscripts n previous time step n+l new time step phase
PAGE 9
CHAPTER 1 INTRODUCTION Background Two phase systems are essential in many engineering applications as well as in nature. Rainfall, fog, and fermentation are all considered as two phase systems. Industrial applications utilizing two phase flows are numerous. Mechanical applications include combustion, fluidized beds, nuclear reactor cooling, foams, jets and sprays. Distillation, absorption and floatation are also considered to be two phase systems. Many two phase flows, such steam and water, may exist in various flow regimes. Bubbly flow is characterized by discrete bubbles in a continuous fluid. Slug flow contains very large, discrete bubbles which nearly fill the cross sectional flow area. In annular flow a thin liquid film flows alting the walls of the pipe or vessel and steam flows along the core of the pipe. Drop flow is characterized by drops of liquid which flow along the pipe and surrounded by a continuous steam phase. In addition to these flow regimes, the flow regimes may be superimposed, resulting in bubblyslug or annulardrop flow. Of particular importance is the bubbly flow regime. In this regime void fractions can range from that of a single discrete
PAGE 10
bubble to foam or frothy flow, up to the slug flow regime. Engineering applications in utilizing bubble flow include bubblecolumns to promote mass transfer, distillation, foams, sprays, cryogenics and nuclear reactor cooling, just to name a few. The heat, mass and momentum.transfer capabilities of the bubbly regime are influenced by the coalescence of the bubbles in the flow. Bubble coalescence plays a significant role in determining the size, distribution and interfacial area of bubbles in many devices and systems. In the cooling of a nuclear reactor, the coalescence of bubbles can be detrimental because the process decreases the surface area of the bubbles, decreasing the heat transfer capabilities. In other applications, such as in separators, bubble coalescence is advantageous since it produces larger bubbles with higher rise velocities. The flow structure around single bubbles has been studied for years. Data have been collected and empirical relationships derived for the drag force seen by the bubble and its terminal rise velocity. The flow field around a bubble has been studied to a lesser extent. The exact structure of the wake(s), its size, shape, whether it's laminar or turbulent, is the topic of many .: fnvestigations. The use of the modern computational techniques has been employed in many codes to model the complex two phase flow problem, by use of the governing mass, momentum and energy equations. However, the understanding of the small scale processes, such as 2
PAGE 11
bubble coalescence, is essential to produce meaningful results on a macroscopic scale. The purpose of this paper is to make a presentation of the numerical model used in the simulation of the phenomenon of bubble coalescence and to show its accuracy by comparison to existing measured data. An existing one dimensional, two phase flow program was altered to include the effects of bubble interaction, i.e., bubble coalescence. The existing code structure included the basic governing equations (Chapter 2), but did not account for accurate rise velocities of individual bubbles nor the interaction of the bubbles themselves. The program is a mechanical code and does not consider heat and/or mass transfer between the phases or internally to the phase. Additions to the existing program included interactions of bubbles with each other as well as the continuous phase. Bubbles are individually tracked throughout the flow field for the selected problem time and area of interest. Accurate modeling of the drag coefficients was included to predict terminal rise velocities of bubbles . A wake model was developed to accurately simulate the approach of a trailing bubble to a leading bubble. The wake model was developed based on data in the bubbly flow regime. Slug flow, annular, drop and vapor flow were not modelled. 3
PAGE 12
CHAPTER 2 PROGRAM STRUCTURE Governing Equations The existing code is a one dimensional program which solves the governing equations for a separated, two phase flow. The continuous phase is modelled as an Eulerian system on a fixed grid. The discrete, or particle phase, is based on a Lagrangian system, as each particle(bubble) is individually tracked throughout the flow field. The two basic equations which govern the particle phase consist of the momentum and mass continuity equations. A set of each of these equations is included for each particle. The momentum equation for a bubble is used in the following form du aP pbVb dEb = Vb(ax)b + pbVbg fo(ub uav) fwub f(grb + u 2.1 The term on the LHS side of the equation represents the inertia effects of the bubble. The first term on the RHS of the equation represents the force on the particle due to the pressure gradient of the continuous phase. The pressure gradient, as seen by the particle, is given by
PAGE 13
2.2 where the inertia and acceleration terms are based on the average continuous phase properties near the bubble. The second term on the RHS of the equation represents the gravity force on the bubble. The third term represents the drag force on the bubble due to the relative motion of the bubble and the continuous phase fluid. The fourth term in the RHS of the equation represents the friction force 6f the wall of the container on the particle. The last term is the added mass force. This term arises when the bubble is accelerated relative to the continuous phase and it sets up a two dimensional flow around it which possesses kinetic energy. Extra energy is required to move the particle. The extra energy requirement creates an extra force seen by the particle. The governing mass continuity equation for a bubble is which takes the form 2.4 The governing momentum equation for the continuous phase is given by 2.5 5
PAGE 14
where the LHS side of the equation is the substantive acceleration of the fluid which consists of the local .and convective contributions to the flow. The first term on the RHS is the pressure gradient of the continuous phase, followed by the gravity force. The third and fourth terms represent the effects on the continuous phase from the drag and added mass forces on the bubble{s). The final term represents wall friction for the continuous phase. The continuous phase mass continuity equation takes the form 2.6 The void fraction, a, for the continuous phase is defined as the total volume fraction of the continuous phase fluid minus the volume fraction of the bubbles in a particular continuous phase grid {cell), i.e. if the void fraction of the continuous phase is 1, there are no bubbles present at that grid location. The void fraction of the continuous phase, a, is based on the sum of the void{s) caused by the bubbles and is given by: a = 1 llib where the void caused by a bubble, ab' is defined as vb ab = o{s) vcell 2.7 2.8 and o{s) is the percentage of the bubble volume located in a particular grid {cell). The percentage, o{s), is based on a fourth order polynomial which spreads the bubble volume over a specified 6
PAGE 15
number of cells in order to obtain a smooth interpolati6n of the discrete phase into the continuous phase. Certain fluid properties (density) for the discrete and continuous phase are determined by the state equation Program Scheme 2.9 The flow of the original code was not altered to include the modelling of bubble interaction. A brief description of the program flow is provided below, including the new features. Figure 2.1 is a fldw chart of the program, with altered and new subroutines indicated. The basic inputs to the code include: material properties for both the continuous and discrete phases, initial particle geometry/position, fluid parameters and grid specifications. The initial continuous and discrete phase densities are first determined in subroutine STATE using equation 2.9. The particle volumes (equation 2.4) and void fractions (equations 2.7 and 2.8) for the continuous and discrete phases are then computed in subroutine PMASS. The drag components for each particle are calculated in subroutine DRAG for later use in subroutine PMOM. The DRAG subroutine calls several new subroutines (see Figure 2.1) to determine wake geometry, wake velocities, bubble placement within the wake, and the overall drag component. 7
PAGE 16
MAIN INPUT Fluid Propert;.a Geometry STATE New Densities MERGE PM ASS Bubble VeJoc:ity, New Bubble Volume. Volume, Oenait)'. Politk)n and radius New Void F'raetion J. OUTPUT Onset at Print Output Sl"'l FJo.? no y .. lime to Quit. _,l es WAKE_LINTER no B u bble Position in Wake and ATea of Intersection DRAG _,l Drag F'orce on each SECWAKE/SECPOS Subbl Secondary Woke Radius Secondary Wak e Velocity PMOM _l PRIMV Now Bubble Velocity New Bubbl Po1ition PrWnory Wake Vloc ity J J CMOM Cont i nuoua PhCH Velocity (;niliol) END l CMASS New Cont inuow PhoN Preaure New or altered l subroutine CMOMFN New Contmuoua Phase Vetocit i n (final) Fi g ure 2.1. Pro g ram fl ow chart 8
PAGE 17
The particle velocity is then calculated using equation 2.1 in subroutine PMOM. The continuous phase velocities are explicitly estimated in CMOM using equation 2.5 (estimated only, since the new time pressure, Pn+l, is not determined yet, and the old time pressure, Pn, is used}. Next, the is implicitly determined in CMASS using equation 2.6, with the voids and velocities in the equation expressed in terms of pressure. The continuous phase velocities are then corrected for the pressure in subroutine CMOMFN. At this point, all continuous phase and discrete phase variables have been determined. The new bubble positions are then checked in MERGE for a possible coalescence, merged if required, and final bubble properties are determined. 9
PAGE 18
CHAPTER 3 NEW CODE FEATURES Bubble Hydrodynamics Bubbles rising in a Newtonian fluid are generally considered to belong to three broad regimes. The first regime is the spherical regime, which usually applies to small bubbles {less than 1 mm in diameter} that are actually spherical in shape or are only slightly distorted. Ellipsoid bubbles, the second regime, cover medium sized bubbles whichare usually oblate in shape, having a convex surface over the entire bubble (viewed from the inside). The last regime, spherical cap bubbles, applies to larger bubbles that, in general, are spherical in the fore geometry and have a flat or skirted tail in the aft. Figure 3.1 depicts the three bubble regimes and the variety of bubble shapes which may fall in those regimes. Three dimensionless groups {Clift, Grace and Weber 1978} are used to establish which regime a bubble falls under: the Eotvos number, the Reynolds number and the Morton number. Re = pdu J.l. 3.1 3.2 3.3
PAGE 19
Spherical Bubbles 0 Elli spsoid Bubbles 1 { ____ J Spherical Capped Bubbles a Figure 3.1. Bubble shapes 11
PAGE 20
Figure 3.2, adapted from Clift, Grace and Weber (1978), shows the relationship between these three dimensionless groups and bubble regimes. The spherical regime includes bubbles with low Re, low Eo and variable M. Bubbles are considered to be ellipsoid at relatively high Re and Eo, and low values of M. For a bubble to be considered a spherical cap, it must have both a relatively high Re and Eo. Spherical bubbles can be truly spherical in shape, or be slightly distorted. Hetsroni (1982) considers bubbles to be spherical when the aspect ratio (width to height) is a maximum of 1.1. Spherical bubbles are dominated by surface and viscous forces (therefore low Re and Eo), while inertia and gravitational effects are negligible. At very low Re, a small bubble may be considered as a rigid sphere and bubble motions are predicted using Stokes law. In an ultra pure system, internal circulation of the bubble reduces the pressure and skin drag on the bubble, resulting in an overall lower drag force of than of a true rigid sphere. In reality, surface contamination and slight deformations tend to increase the drag on the bubble. Surface contamination is generally swept to the rear of the bubble, generating a surface contamination concentration gradient. This surface concentration gradient causes a surface gradient in the shear forces that opposes the surface motion, increasing the drag. With a higher Re, inertial effects become more important, leading to the onset of bubble deformation. 12
PAGE 21
Bubble regimes adapted from Clift, Grace and Weber (1978) Figure 3.2. Bubble regimes 13
PAGE 22
Ellipsoid bubbles are usually considered to be medium sized bubbles, having an equivalent diameter (diameter of a spherical bubble with equal volume) ranging from 1 to 17 mm. The complex interaction of viscous, surface, gravitational and inertial effects causes a variety of "ellipsoid" shapes and motions to occur. Ellipsoid bubbles may be oblate, prolate, or may lack for and aft symmetry and wobble or oscillate while they rise. Ellipsoid bubbles are characterized by the presence of secondary motion. Two types of motion exist; primary rigid body motion, where the bubble can rock from side to side, and have a zig zag or helical trajectory, and oscillatory motion (shape dilations). The two motions may be superimposed leading to extremely complex patterns. Smaller bubbles in the ellipsoid regime tend to be slightly oblate. The increasing oblateness tends to promote the development of an attached wake, and eventually wake shedding. Surface contaminants can also play an important role in wake generation and wake shedding by damping out internal circulation. Wake formation generally occurs at a Re of about 20 and wake shedding at Re of 200. However, if significant deformation has occurred, the of wake generation/wake shedding may occur at significantly lower values. Alternatively, interfacial mobility in purified systems significantly delays theformation of an attached wake and wake shedding. The onset of the secondary motion coincides with the onset of wake shedding . Oscillations disrupt the internal circulation 14
PAGE 23
patterns, leading to an increased drag, which in turn, decreases the terminal velocity. The frequency of the oscillations and wake shedding do not generally coincide. In fact, evidence suggests (Otake, Tone, Nakao, and Mitsuhasi 1977) that when the frequency of the oscillations and wake shedding do coincide, a resonance can occur contributing to bubble breakup. The spherical cap regime (de > 17mm) covers most bubbles of engineering interest. These bubbles are important in fluidized bed, liquid and metal processing and underwater explosions. Spherical capped bubbles may have a truly spherical or ellipsoidal leading edge and usually have a flat trailing edge. Dependent upon the flOid and dynamic properties, the trailing edge may take on different forms: for lower Re and higher M, the bubble may be dimpled on the trailing edge, and for slightly higher Re, a skirt may be formed on the trailing edge (See Figure 3.1). The trailing skirt occurs when viscous forces on the outer rim are sufficient to over come restraining effect of surface or interfacial tension forces. The skirt thickness decreases with the length of the skirt and the length can increase with an increasing Re. The skirts may be smooth or wavy, but has little influence on the drag (and terminal velocity) of the bubble. Surface contaminates have little effect on the drag of spherical bubbles. Spherical capped bubbles may have a closed primary wake at lower Re, or an open, turbulent primary wake at higher Re. 15
PAGE 24
Drag Coefficients and Transition Criteria The previous section discussed different bubble regimes considered for bubbles of different sizes in different liquids. Bubble shapes are influenced by the size of the bubble as well as the physical and dynamic properties of the bubble and its surrounding fluid. One of the most extensive experimental studies conducted in this area was performed by Peebles and Garber (1953). In order to properly model the motion of a rising bubble, accurate predictions must be made of the drag force seen by the rising bubble: The numerical model for the drag coefficients and associated transition criteria were based primarily on the work of Peebles and Garber (1953). Bubbles were divided into the three previously discussed regimes: spherical bubbles, ellipsoid bubbles and spherical cap bubbles. As previously mentioned, cylindrical bubbles, slug flow, annular and vapor flow regimes were not included I in the numerical model. Spherical Regime The spherical bubble regime is actually divided i nto two subregions. Although all bubbles in this regime are considered to be spherical, the underlining assumptions for the two subregions are different. Figure 3.3, Kelley (1989) indicates that the relationship between Re and the drag coefficient, CD, is decreasing for increasing Re. For subregion 1 (known as Region 1.1), the bubble is assumed to 16
PAGE 25
1E+02 4 2 ..... 1E+01 c Q) u \f.'+Q) 4 2 0 u 1E+OO Ol ....., 0 L. 0 2 Ellipsoid Spherical Cop 1 E 0 1 I I I I I I I I I I I I I I I I 0 200 600 1000 1400 1800 2200 2600 3000 Reynolds Number Drag Coefficient (Kelley 1989) vs. Reynolds number, showing bubble regimes. Figure 3.3. Reynolds number vs. drag coefficient
PAGE 26
obey Stokes Law for a solid sphere. Stokes solution is one of the oldest known solutions for a creeping motion. A sphere is assumed to be in a uniform flow, with the normal and tarigential velocity components to the bubble surface are equal zero (no slip boundaries). Inertia forces are considered to be small in comparison to viscous forces and gravity is assumed to be the only external force acting on the bubble. Stokes solution solved the NavierStokes equations and continuity for an incompressible fluid where inertial terms were i.e., grad p = u = 0 where: u = ui + vj + wk 3.4 3.5 The total drag over the bubble was obtained by integrating the pressure distribution and shear stress over the surface of the bubble yielding: 3.6 Using the customary form of the drag coefficient (obtained by dividing the drag force by the kinetic energy (dynamic head) and the cross sectional area), D CD= 3.7 and combining it with equation 3.6, the drag coefficient for a solid sphere was found to be: CD = 24/Re 3.8 The Re in this case is defined as Re = 18
PAGE 27
The terminal velocity of a bubble in this region is found by equating the buoyancy force and 3.6. 3.9 The upper limit for the Stokes region was experimentally found by Peebles and Garber (1953) to be Re = 2. The second subregion (Region 1.2) also assumes that the bubble is of spherical shape, but is tending to become ellipsoid. It was observed by Peebles that the drag coefficient for this subregion was slightly lower than the predicted drag coefficient from Stokes solution. The drag coefficient, similar to the Stokes drag, takes into account both form and skin drag and was experimentally determined by Peebles and Garber (1953) to be: 18.7 c = ,::::= D Re0.68 3.10 The terminal velocity of this region was determined by equating equation 3.10 with the assumed drag coeff .icient for the bubble's terminal velocity. This drag coefficient was determined by assuming that at terminal velocity the buoyancy force on the bubble, was equal to the drag force (equation 3.7) resulting in 8 rbg CD = 3 Ub2 3.11 Combining equations 3.10 and 3.11, the terminal velocity is given by: 19
PAGE 28
3.12 The upper limit for region 1, in terms of a transition Re, was determined based on the drag coefficient relations for Region 2 to obtain a smooth transition. Bubble behavior in Region 1 depends on bubble size as well as liquid and gaseous properties. Bubble behavior is fairly easy to predict, up until a point where the relationship between the viscous, surface, gravitational and inertia forces cause the bubble to cease behaving like a solid sphere, and begin to distort. The drag coefficient used in the numerical model for Region 1 was taken as the maximum of 3.8 and 3.10. Ellipsoid Regime Ellipsoid bubbles are also broken into two subregions. Figure 3.3 shows that in region 2, the c0 increase with increasing Re. In the spherical regime, the drag coefficient was not radically different than the drag on a solid sphere. Figure 3.j shows the the drag coefficient departs radically from the solid sphere drag in the ellipsoid regime. As the bubble distorts, secondary motions begin to occur, and the drag on the bubble is greatly increased. Peebles and Garber (1953) determined the drag coefficient in this region by first making some assumptions about the terminal velocity. The terminal velocity was assumed to be some function of 20
PAGE 29
the bubble radius, liquid viscosity, density, surface tension and gravity. Through dimensional analysis the following relationship was obtained: 3.13 where a and d were determined from experimental data to be a = 1.91 and d = 0.5. Solving for the terminal velocity U = 1.35 (u)0 5 prb 3.14 Equation 3.14 indicates that the terminal velocity in this region is independent of the fluid viscosity. The drag coefficient, found by substituting equation 3.14 into equation 3.11, is then c0 = pu 3.15 To define the transition between region 1 and region 2 the drag coefficient relations were equated and solved for the governing Re: Re = pu3 3 .16 From Peebles work, bubbles in Region 2.2 showed a marked tendency toward a constant terminal velocity, independent of bubble size. Figure 3.3 shows that the drag coefficient in this subregion still increase with increasing Re, but not as sharply as in Region 2.1. The fluid viscosity is assumed to have negligible effect on the terminal velocity in this region. 21
PAGE 30
Using the same procedure as before, Peebles and Garber {1953) determined a relationship for the terminal velocity for this region: u = b p Similarly, the associated drag coefficient is CD= 0.82G1 0 25Re where G = 1 pa3 3.17 3.18 The relationship for the drag coefficient determined by Peebles and Garber {1953) was found to give a terminal velocity which was too low. This was probably because the influence of the walls was not considered in the development of the empirical drag relations {Harmathy 1960). Therefore, the relationship suggested by Harmathy {1960) was used to determine the drag for this region: CD= 0.488 G 125Re 3.19 The drag coefficient for region 2 used in the numerical model was taken as the maximum of equations 3.15 and 3.19. Studying equations 3.15 and 3.19 one can see that the combination of the Re and the dimensionless group G 1 effectively cancelled out the viscosity. Therefore, the motion of the bubble in this region is independent of viscosity. Physically, this indicates that the drag force is governed by distortion and the swerving motion of the particle, and the change of the particle shape is toward an increase in the effective cross section {Ishii and Chawla 1979). 22
PAGE 31
Upon closer examination of the equations in region 2.1, if equation 3.13 is rearranged, it is seen that it is equivalent to a constant ratio between inertial and surface forces which acts to deform the bubble. Spherical Cap Regime In the spherical cap region, the drag coefficient is assumed to be constant with a value of 8/3. The transition to this region was determined by equating the drag coefficient from region 2.2 (equation 3.19) and the 8/3 value, resulting in Re = 5.47G1 0 25 3.20 The dragcoefficients and associated transition criteria summarized in Table 3.1. Bubble Geometry Relationships for bubble shape are considered for use in flow visualization studies. For calculational purposes, all bubbles are considered to be spherical, their volumes calculated from mass continuity (equation 2.4), and equivalent radius determined by Equation 3.21 is used to determine the size of bubbles in region 1. Ellipsoid bubbles are assumed to maintain a constant 3.21 aspect ratio a/b, (See Figure 3.1), although in reality they are distorting from a spherical shape to a spherical cap. From the work 23
PAGE 32
ShaQe Regime Drag Coefficients Transition Criteria 24 Co = Re Spherical (maximum) 18.7 c D Re0.68 Re = 4.02 pa3 C = 0 pa3 Ellipsoid (minimum) c0 = 0.488 G 125Re Re = 5.47G1 0 25 Spherical Cap 8 3 Table 3.1. Summary of Drag Coefficients and Transition Criteria 24
PAGE 33
of Kelley (1989), ellipsoid bubbles distort from an aspect ratio of 1 to an aspect ratio of 4.02 for spherical caps. Therefore, an ellipsoid bubble is assumed to have an aspect ratio of 2, and the equivalent diameter is 3.22 Spherical cap bubbles, having an aspect ratio of 4.02 have and equivalent diameter of 3.23 Bubble Coalescence The coalescence of bubbles is important in determining the size, distribution and interfacial area of bubbles in gasliquid systems. Bubble coalescence influences the bubble induced flow and mass transfer rate. Bubble coalescence occurs in three stages: 1. A trailing bubble enters into the wake, or partially into the wake of a leading bubble. 2. Since the relative velocity of the fluid in the wake is greater than in the surrounding fluid, the trailing bubble experiences less drag in the wake and approaches the leading bubble until they collide. 3. After the collision, only a thin film separates the bubble. A film thinning process takes place until coalescence occurs. Only the first two stages of bubble coalescence are considered for the purpose of this paper. The wake structure behind rising 25
PAGE 34
bubbles is clearly important in the coalescence process. The wakes behind spherical, ellipsoid and spherical cap bubbles were briefly discussed earlier. A single model was established for the wake structure behind all bubbles for all three regimes. The model is based on the wake behind a spherical cap bubble. Spherical capped bubbles were selected due to their importance in many engineering applications (and because of the abundance of data for this regime). Flow visualization studies have clearly indicated the presence of a wake behind isolated spherical capped bubbles. The structure of the wake, is composed of two parts: a primary wake and a secondary wake. The primary wake is an attached, closed (at low Re) wake that exhibits a recirculation pattern, or toroidal movement in the wake. The primary wake is separated from the external flow by a single streamline. In previous studies, the primary wake shape has been assumed to have the shape equivalent to completing the sphere or I spheroid of the spherical cap (giving a spherical or oval shape), or as having a conical or exponentially decaying shape. The velocity of the recirculating fluid contained in the primary wake may have an upward velocity that is greater than the rise velocity of the bubble (relative to a fixed frame). The primary wake structure is important in the coalescence process of bubbles in close range (within a few bubble. radii). The recent study by Bessler and Littman (1987) provided a good description of the wake structure behind spherical cap bubbles. 26
PAGE 35
Bessler and Littman investigated the wake structure behind circular capped bubbles using aspirin traces. They discuss the general wake structure adjacent to the primary wake as a free shear layer which contains large scale vortices (See Figure 3.4) generated near the bubble rim that remain essentially stationary in the lab frame. These vortices dissipate energy, growing in size. In more viscous fluids, the shear layers develop more rapidly and the size of the primary wake decreases. Bessler and Littman (1987) found the shape of the primary wake to be an airfoil shape with a cusped tail due to the growth of the large scale vortices in the free shear layer. They found the length of the wake to be about six radii of the leading bubble for a glycerin solution and 11Slightly longer" for water. The secondary wake represents the general bulk movement of the continuous fluid, due the momentum defect caused by a rising bubbles passage. It is the secondary wake motion which is responsible for the coalescence of bubbles with large separation distances. Wake Structure Primary Wake Two primary wake mQdels were investigated. Both models were based on the work of deNevers and Wu (1971). Two wake shapes were suggested in this study; exponential and conical (see Figure 3.5). The exponential wake shape was of the form 27
PAGE 36
Internal Circulation Primary Wake Length Spherical Copped Bubble Primary Wake Boundary Wake Figure 3.4. Wake structure behind a spherical cap bubble (Bessler and Littman 1987) 28
PAGE 37
where: rpw = rb e kx R = primary wake radius Rgw = bubble equivalent radius x = distance between bubbles k = decay constant, determined 1 experimentally to be 0.28 em3.24 The conical primary wake assumed a wake length based on the equivalent radius of the leading bubble. The length of the wake, Lw, was determined using a dimensionless wake length, Lw which remains constant for a particular fluid. For water, Lw = LwfRb = 10, as determined by deNevers and Wu (1971) (Lw = 7 and k = 0.37 for glycerin). Two models were considered for the velocity structure in the primary wake: 1) the velocity throughout the primary wake was identical to that of the leading bubble and 2) the toroidal vortex motion, described earlier, was assumed to give a center line velocity that initially rose to a value greater than that of the leading bubble (relative to a fixed frame) then decreasing to the velocity of the leading bubble. The toroidal motion model was based on measurements of the primary wake velocity by Bhaga and Weber (1980) where the vortex center was found to be located approximately at a distance 60% of length of the primary wake; The velocity in the wake was assumed be the same as the center line increasing from the tail of the primary wake (which is moving at the same speed as the leading bubble since the wake is attached) to a value of 1.5 the bubble 30
PAGE 38
velocity, ub, at the vortex center. The velocity then decreases back to ub (see Figure 3.6). Secondary Wake The secondary wake model was based on the asypmtotic wake solution described by Batchelor (1967) where a long distance downstream from a body, the wake obtained is independent of the geometry of the body and is determined by the drag on the body. The wake velocity along the center line decreases as the reciprocal of the distance behind the body; and is distributed as a Gaussian curve in the radial direction. The wake induced flow is taken from Stuhmiller, Ferguson and Meister (1989) and is uw(x,r) = uwc(x) e(r/r*)2 where: u (x,r) =wake velocity at r,x. uw (x) line wake velocity c x = axial distance t = radial distance r = wake radius at x 3.25 Several relationships for the center line velocity of the secondary wake were investigated. The first center line velocity relationship considered (Stuhmiller, Ferguson and Meister 1989) was 3.26a The relationship was based on the collection of center line wake velocity data from several publications. 31
PAGE 39
Normalized Wake Velocity, .!:!wake 1 5 vortex center = 0 6 L. o.s 1 vortex center X Normalized Wake Length, L,. Figure 3.6. Primary wake velocity 32 1 0
PAGE 40
Additional center line velocity relationships (listed below) were given by Crabtree and Bridgewater (1971), Bhaga and Weber (1980) and Batchelor (1967), respectively. 3.26b ubB uwc = 3.26c ( + )1.43 3.26d The radius of the secondary wake was determined by using the conservation of momentum flux. Far downstream of the leading bubble, the total momentum is related to the total drag on the bubble. It is assumed that a long distance downstream of the bubble, pressure differences in the longitudinal and traverse directions are small. Therefore, the loss of momentum due to friction J = p f ucw(ub ucw)dA is equal to the total drag on the bubble 3.27 3.28 The drag coefficient was determined by assuming the buoyancy force and the drag force were equal (Equation 3.11). Integrating equation 3.27, assuming that ucw2 is small, yields 3.29 33
PAGE 41
Equating 3.28 and 3.29 gives a relationship for determining the secondary wake radius: 3.30 Therefore, at distances far downstream, the secondary wake velocity may be determi"ned by using equations 3.25, 3.26, and 3.30. At axial locations closer to the leading bubble, the fluid motion is extremely complex and the assumption of negligible pressure differences may not apply. Equation 3.30 cannot be applied to determine the secondary wake radius . The axial location above which the conservation of momentum flux is not applicable was . assumed to be at the point where the primary wake was of no influence, i.e. at the tail of the primary wake. A simple model for the secondary wake radius above this point was used, based solely on the geometry of the leading bubble. The radius of the secondary wake at Lw was first determined using equation 3.30, assuming pressure differences were negligible. The radius of the secondary wake above that point was then determined assuming a linear relationship between the radius and Lw and the bubble radius. The velocity structure of the fluid located in the secondary wake above Lw was considered separately for the primary and secondary wake. No attempts were made to model the vortices in the shear layer. The velocity along the streamline of the primary wake is equal to the velocity of the leading bubble. The velocity at the 34
PAGE 42
"edge" of the secondary wake approaches the free stream velocity. therefore, the velocity in the secondary wake above Lw was assumed to have a profile similar to equation 3.25, decreasing outwardly from the velocity of the leading bubble. The complete model of the wake structure, assuming a conical primary wake shape, is shown in Figure 3.7. Drag on a Bubble in a Wake The total drag on a bubble caught in, or partially caught in, the wake of another bubble is determined based on the various fluid velocities "seen" by the bubble. The location of the following bubble within the wake or wakes is determined first. The area of intersection of the bubble frontal area and the wake are then calculated. .Finally, the area ratio of the area of intersection to the bubble frontal area is calculated (the sum of all area ratios is equal to one for each bubble). Area ratios were used in an attempt to smooth or eliminate discontinuities in the total drag on the bubble as it entered the region of the wake(s). If a following bubble is effected by the wake of more than one leading bubble, only the leading bubble with the highest (combined primary and secondary) area ratio is considered. If a following bubble is completely within the wake of one or more leading bubbles, the bubble closest to the following bubble is considered. 35
PAGE 43
Internal Circulation Primary Wake Length Secondary Woke Velocity Profile i/ Spherical Copped Bubble Primary Wake Primary Wake Boundary Moving at ububble Velocity Profile near Primary Woke +Secondary Wake Figure 3.7. Wake structure behind a spherical cap (code model) 36
PAGE 44
The area ratios are used to determine the components of the total drag force seen by the bubble. A drag component is calculated for each fluid stream seen by the bubble . A drag coefficient for each component is calculated using methods described earlier, and based on the particular fluid stream velocity. The fluid velocity seen by a bubble intersecting the primary wake is assumed to be that 6f the center line velocity (Figure 3.6). The fluid velocity seen by a bubble intersecting the secondary wake is a calculated average velocity based on the frontal area of the bubble in the wake and its location in the wake. Drag components for a bubble partially within the primary wake and secondary wake of a leading bubble, as well as partially in the free stream, are calculated below. Fif = 0.5pC01ub ufiAf FIP = o.5pc021ub uwpiAp Fis = 0.5pC031ub uwsiAs where the subscripts f,p and s refer to the free stream, primary wake, and secondary wake, respectively. 3.3la 3.3lb 3.3lc The drag components are used in the momentum equation (2.1), in subroutine PMOM to solve explicitly for the bubble velocity. Bubble Merging The merging of two or more bubbles is assumed to occur if two bubbles touch. The film thinning process is not modelled. The 37
PAGE 45
equivalent radius, equation 3.20, was used to determine if a bubble touched or overlapped another bubble. The new, "merged" bubble required as estimate of its new mass, volume and velocity. A simplistic model of the conservation of linear momentum and mass was used to make an initial guess of the new velocity: 3.32 The size of the merged bubble is determined by assuming the new density is a weighted average of the individual bubbles and total mass is constant 3.33 Although equation 3.32 assumes a constant linear motion, a new horizontal (y) position for the merged bubble is calculated. The new bubble position is calculated based on the relative volumes of the two individual bubbles 3.34 38
PAGE 46
CHAPTER IV RESULTS AND DISCUSSIONS Comparison between calculated values and experimental data are presented first to show the accuracy of the numerical models added to the existing code. The drag coefficient calculations are tested by comparing the calculated terminal velocity of individual bubbles and experimental data. Second, wake structure models are tested by comparing calculated bubble coalescence times with those from experimental data. Finally, a flow visualization of the coalescence of bubbles rising vertically in a pipe is presented. Terminal Velocity The terminal velocities of bubbles with varying radii were calculated using the drag coerficient relationships and transition criteria discussed in Chapter 3. The discrete phase for the calculations was water and the particle (bubble) phase was air. Figures 4.1 and 4.2 provide a comparison of calculated terminal velocities with measurements of Peebles and Garber (1953) and with that of Haberman and Morton (1953). Also shown in the figure are the terminal velocities calculated by Peebles and Garber (1953). The calculated values agree well with the calculated values by Peebles and Garber (1953), as well as the measured data. Note that
PAGE 47
0 U) E 0 >, ....._. 6 4 3 2 .... ,..u1E+01 0 Cl) > a c E L Cl) I/ 3 x measured o calculated 2 ecode calculated 1 E +00 I I I I I I I I I I I I I I I I I II I II 1E02 2 3 4 1E01 2 3 4 1E+OO 2 3 4 Bubble Radius, em Comparison of code calculated terminal velocities of bubbles in water and measured velocities by Peebles and Garber (1953). Also shown are terminal velocities calculated by Peebles and Garber (1953) using equations in Chapter 3. Figure 4.1. Terminal velocity of bubbles in water (Peebles and Garber 1953)
PAGE 48
"'" 6 4 en 3 .......... E 2 0 >. .......... u 1E+01 0 Q) 6 > a 41 c 3 E 2J XX L Q) tv 1E+OO I 1E02 X x measured ..x code calculated I I 2 3 4 1 E01 2 3 4 1 E+OO 2 3 4 Bubble Radius, em Comparison of code calculated terminal velocities of bubbles in water and measured velocities by Haberman and Morton (1953). figure 4 2. Terminal velocity of bubbles in water (Haberman and Morton 1953)
PAGE 49
the calculated data lies above the experimental data for bubbles with equivalent radii greater the 0.3 em. These values were considered to be low, as noted in Chapter 3, and a different drag relationship was used for bubbles in this regime. Additional comparisons to measured data are shown in Figure 4.2. The measured data was extracted from Haberman and Morton (1953). Excellent agreement is shown for all bubble radii, except where the transition from a spherical bubble to an ellipsoid bubble occurs. A possible explanation is that the water used in the Peebles measurements was contaminated, resulting in lower terminal velocities Coalescence Times The wake model was compared to data gathered by Crabtree and Bridgewater (1971). Crabtree conducted experiments on the relative motion of vertically aligned air bubble pairs, having volumes from 10 cm3 to 40 cm3 in a 67 percent wt solution of sucrose in water. The bubble pairs were tested up to 70 em apart. The density and dynamic viscosity of the sucrose solution were given, but the surface tension and speed of sound were not provided, so estimates of those properties were made. A few test runs were made to determine the sensitivity of the estimate properties. These runs showed that varying the surface tension and speed of sound did not effect the bubble velocities drastically. 42
PAGE 50
Original Model The first wake structure model attempted, used the following 1. Conical/exponential primary wake shape (L = 10, k = 0.28) 2. in the primary wake = U 3. Velocity in the secondary wake in Eqn 3.23a Figure 4.3 shows the coalescence of a bubble pair with an initial separation distance of 15 em. The leading bubble volume was 30 cm3 and the following bubble volume was 25 cm3 Results from both the conical and exponential primary wake model are presented. The leading bubble is unaffected by the following bubble. Its rise velocity compares well with the experimental data. The following bubble is initially affected by the secondary wake only, entering the primary wake at around t = 0.3 seconds. The following bubble data compare fairly well with the experimental data; however, near coalescence the calculated velocity for the following bubble is somewhat less that the velocity indicated in the experimental data. This could suggest improper modelling of the primary wake. Although not shown, additional runs were made assuming the drag component on the following bubble was zero for the area of the bubble considered to be within the wake. The results of this change were small, although calculated coalescence times were shorter. The conical model appears to given better results than the exponential model. The coalescence times of bubbles of the same size and for a smaller leading bubble and larger following bubble are presented in 43
PAGE 51
0.80 0.70 0.60 ,.... E u 0.50 .40 :.;:; "en 0.30 0 n.. 0.20 0.10 0.00 .00 p ,. ............ ,."' ,.x',.... . ... ;' ./ ; .... calculated (exponential} )(measured measured code calculated ......code calculated (conical} I I I I I I I I I I I I I I I I I I I I I I I I I I ITITnl .1 0 .20 .30 .40 .50 .60 70 .80 .90 1 .00 Time, (seconds) Bubble position vs. time of a trailing bubble {25 cc) approaching o leading bubble {30 cc). in sugar water. using the original code model. Measured data from Crabtree and Bridgewater { 1971 ). Figure 4.3. Bubble positions during coalescence original model
PAGE 52
Figures 4.4 and 4.5, respectively. Again, the conical model gave better results than the exponential. In both cases, the model predicted the coalescence times closely to the experimental data for bubbles with an initial separation distance of 40 em or less. For bubbles with an initial separation distance of greater than 40 em the calculated coalescence time were longer than the experimental times. Figure 4.5 gives a more favorable comparison with experimental data. This was expected since the higher rise velocity of a larger following bubble would play a larger role in the coalescence than the wake of the smaller, slower leading bubble. Final Model Due to the relatively poor agreement of the data for both the trailing bubble positions and the coalescence times, variations of the model were made. First, a comparison of the available secondary wake velocity relationships was made, selecting a consistent bubble size and terminal velocity for comparison. The asypmtotic wake velocity clearly gave too high a wake velocity (this is reported by several authors). Equations 3.26b and 3.26c gave a similar velocity profile, while equation 3.26d yielded a slower velocity. After many comparison runs, Bhaga and Weber's (1980) relationship appeared to give the most reasonable results. This relationship appeared to be more appropriate than the others, since the velocity is based on properties of the fluid and the bubble size, not just 45
PAGE 53
.,.. en 7
PAGE 54
...... 6 0'5 Q) (/) .._...,
PAGE 55
the bubble size. The conical primary wake shape was judged to give better results and it was retained for the final model. A value of Lw = 7 for the sugar water, as opposed to a value of 10 for pure water as determined by deNevers and Wu (1971), was found to geve better results. The primary wake velocity was assumed to follow the recirculation profile given in Figure 3.6. Figure 4.6 shows the coalescence of the same bubble pair as in Figure 4.3, but using the final model. The calculated rise velocity of the leading bubble compares well with the measured rise velocity of the leading bubble, indicating the drag coefficient used for the spherical cap was adequate. Additionally, the linear rise indicates that the following bubble has no effect on the rise velocity of the leading bubble. The calculated positions of the following bubble agree quite well with measured data throughout its approach toward the leading bubble. The trailing bubble enters the primary wake of the leading bubble at about t = 0.6 sec., after which the rise velocity of the trailing bubble increases somewhat, but not to the extent of the measured data. Near coalescence, the calculated position of the trailing bubble is slightly less than the measured data. Crabtree and Bridgewater (1971) observed that as the following bubble came within two or three radii of the leading bubble its aspect ratio increased, elongating the bubble in the vertical direction. An increase in the 48
PAGE 56
\D 0.80 0.70 0.60 .r.. E 0 0.50 ..._ .40 "en 0.30 0 0.. 0.20 0.10 0.00 .00 *" measured *me asured code calculated ...code calculated .1 0 .20 11I'TIJ! .30 .40 .50 .60 70 .80 .90 1 .00 Time, (seconds) Bubble position vs. time of a trailing bubble (25 cc} approaching a leading bubble (30 cc}, in sugar water, using the final code model. Measured data from Crabtree and Bridgewater ( 1971 ). figure 4.6. Bubble positions during coalescence final model
PAGE 57
aspect ratio would allow more of the bubble to enter the primary wake, its overall drag, and increasing its velocity. Figures 4.7 and 4.8 are for bubble coalescences in sugar water over larger separation distances. For bubbles of equal volume (Figure 4.7) the calculated coalescence times are much improved over the previous model. The coalescence times match the measured data quite well within 60 em separation distances. The model predicts slower coalescence time at larger distances, approximately 20% slower than measured data at 70 em. The coalescence of a larger trailing bubble and smaller leading bubble (Figure 4.9) were also improved using the new model. Coalescence times are predicted reasonably, even up to 80 em separation distances. Figures 4.6 4.8 gave comparisons of calculated vs measured data from bubble coalescences at fairly large distances with reasonable success. At closer distances the importance of the primary wake is increased. Few experimental data were found for comparisons of at close distances. Experimental data from the work of deNevers and Wu (1971) provided data in two liquids for initial bubble separation distances of five radii. The two liquids considered were glycerin and pure water. The mechanism used to release the bubbles was a horizontal sparger with a single circular hole. As bubbles were released into the stagnant fluid, the motion of the bubble chain tended to cause the general bulk movement of the fluid to have an upward motion. This movement caused the terminal velocities of the bubbles to be 50
PAGE 58
U1 ...... 7 ()6 Q) (/) '::'5 Q) E i4 Q) g3 Q) 0 0 81 =1 10 X X X X X X X +code calculated x measured 20 30 40 50 60 70 80 Initial Separation Distance, (em) Coalescence time vs. initial separation distance for bubbles of equal volume (25 cc), in sugar water, using final code model. Measured data from Crabtree and Bridgewater ( 1971 ). 90 100 figure 4.7. Coalescence of same size bubbles final model
PAGE 59
. U1 N 6 us Q) Ul ......._.,
PAGE 60
increased by as much a 30%. The upward motion, as it will be discussed later, could cause discrepancies between calculated and measured data. The first coalescence event considered was for air bubbles in glycerin (Figure 4.9). The calculated motion of the bubble is somewhat different than the measured, although both coalesce in approximately the same time. The velocity of the measured bubble increases as the separation distance decreases. The calculated velocity increases initially, but then appears to be constant. Effects on the bubble motion by the release mechanism are unknown. However, the primary wake model appears to have adequately predicted the bubble motion in this case. The available data on coalescence times for pure water, as opposed to sugar water or other viscous fluids, were somewhat limited. DeNevers and Wu (1971) obtained data from the coalescence of two bubbles with equal radii in pure water. Attempts were made to match the data,but with limited success. DeNevers and Wu (1971) pictured the bubbles schematically to be spherical capped bubbles. However, the equivalent radii of bubbles for the given bubble volumes, fell in the ellipsoid range. In addition, the author discussed a rocking motion of the bubble (causing the data to scatter) which is characteristic of ellipsoid bubbles. The rocking motion, as discussed next, causes difficulties in predicting the coalescence. 53
PAGE 61
U1 25 L.. '.. ......... :::> 11 20 X Xx xx ......... X lisa< Q)15 >xx x 1.code calculated E XxX iX X I X measured fJ) fll10 Q) c 0 fJ) 5 c Q) E 0 0 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I"T I I I I I I 0 1 2 3 4 5 Dimensionless Distance, X* = del x/ r Time vs. bubble position for a coalescence in glycerin. Measured data from deNevers and Wu (1971). Figure 4.9. Glycerin coalescence 6
PAGE 62
Figure 4.10 shows the calculated and measured coalescence data. According to the author, the trailing bubble was rocking back and forth, moving in and out of the wake of the leading bubble. The bubble was considered to be out of the wake of the leading bubble * from about t = 15 to t = 20. The scatter of data near t = 0 also suggests that the bubble may have been partially out of the wake. If the general motion of the bubble at t = 15 were continued, it appears that the bubble would have coalesced near t = 14 or 15. The calculated data indicated a faster coalescence time, occurring at approximately t = 9. The assumed time of 14 or 15 gives a coalescence time approximately 55% slower than predicted. Part of the discrepancy may be due to the release mechanisms problems, but more likely it is due to the oscillatory motion of the bubble causing it to move in and out of the wake. Flow Visualization Figure 4.11 shows a flow visualization of bubble coalescence in a 20 em diameter pipe over a twelve second time period. At the initial state, t = 0 seconds, three spherical bubbles exist in the lower end of the pipe. Bubbles of various sizes and shapes are released over the first five seconds. The stages of coalescence of the bubbles are shown at times t = 4 seconds, t = 8 seconds, and t = 55
PAGE 63
U1 Ol L 25 +J ::::> X X 11 20 +J Q)15 E i(I) (1)10 Q) c: 0 (I) c: 5 Q) E 0 0 I 0 I_._ code calculated X X X X I X X bubble not in wake measured X XX X X X X X X X X X 1 2 3 4 5 Dimensionless Distance, X* == del x/ r Time vs. bubble position for a coalescence in water. Measured data from deNevers and Wu (1971 ). 6 Figure 4.10. Water coalescence
PAGE 64
0 seconds 6 seconds 12 seconds 10 m 10 m 10 m 20 em 5 m 5 m 5 m 0 0 0 0 m 0 m 0 m Figure 4.11. Numerical simulation of bubble coalescence in a pipe 57
PAGE 65
12 seconds. t = 12 seconds the bubbles have coalesced to the point where the onset of slug flow has begun, and the problem is terminated. Discussion The empirical relationships derived by Peebles Garber (1953) and by Harmathy (1960) resulted in terminal velocities which matched measured data very well. Accurate predictions of the drag force (and terminal velocity) on a single bubble were essential to the accurate prediction of bubble coalescence. The wake model performed well in predicting the coalescence of bubbles from both long and close separation distances. The primary wake model, which at first was assumed to be a major contributor to the coalescence effect, was modelled (and remodelled) _in great detail. However, in the review of the results and the experience gained in the numerous test runs, it is not the dominating process causing bubble to coalesce. The primary wake affects only the end of the approach of a trailing bubble to a leading bubble. It is the secondary wake which plays the dominant role in bubble coalescence. The structure of the primary wake, the internal circulation and the primary wake velocity are, however, an interesting aspect of rising bubbles which help determine the flow structure around the bubble. Many authors have produced excellent photographs of the wakes behind spherical capped bubbles. 58
PAGE 66
In any numerical model one must make a decision on the most important effects of the actual physics, and the level of detail which is necessary to adequately predict reality. Dependent upon the required accuracy, it is possible that the primary wake model may be sacrificed to decrease the complexity of large problems. The modelling of each bubble in the Lagrangian method, as opposed to a homogeneous twophase model, increases the complexity of the problem and requires extra computational time, especially when large numbers of bubbles are tracked. Therefore, any method of decreasing the complexity of the model would clearly be advantageous. The difference in scales between the continuous phase and discrete phase can cause additional problems in modelling twophase flow problems. The level on a more microscopic level of an individual bubble required a different scale than for the continuous phase (failure to recognize this fact or properly account for it can give rise to executional errors in the program which may only be solved by many late night test runs accompanied by caffeine stimulants). Additions to the code could be made to extend not only the range of bubbles regimes (include slug flow, for example) but also to handle more realistic scenarios. Breakup of bubbles could be included in the model, as well as 20 effects of vessel walls. The majority of the available data on bubble coalescence deals with fairly viscous fluids, such as sugar water and glycerin and the empirical relationships for the secondary wake velocity are laminar 59
PAGE 67
in nature, decreasing by x1 Problems for bubbles in turbulent regimes have not been explored. 60
PAGE 68
REFERENCES Batchelor, G.K. 1967. An Introduction to Fluid Dynamics. London: Cambridge University Press . Bessler, W.F., and H. Littman. 1987. Experimental Studies of Wakes Behind Circularly Capped Bubbles. Journal of Fluid Mechanics 185:137151. Bhaga, D., and M.E. Weber. 1980. InLine Interaction of a Pair of Bubbles in a Viscous Liquid. Chemical Engineering Science 35:24672474. Clift, R., J.R. Grace, and M.E. Weber. 1978. Bubbles, Drops and Particles. New York: Academic Press. Crabtree, J.R., and J. Bridgewater. 1971. Bubble Coalescence in Viscous Liquids. Chemical Engineering Science 26:839851. deNevers, N., and J.L. Wu. 1971. Bubble Coalescence in Viscous Fluids. AIChE Journal 17:182186. Haberman, W.L., and R.K. Morton. 1953. David Taylor Model Basin Washington. Report number 802. Harmathy, T.Z. 1960. Velocity of Large Drops and Bubbles in Media of Infinite or Restricted Extent. AIChE Journal 6:218288. Hetsroni, G. 1982. Handbook of Multiphase Systems. New York: McGrawHill. Ishii, M., and T.C. Chawla. 1979. Local Drag Laws in Dispersed TwoPhase Flow. Argonne National Laboratory NUREG/CR1230. Kelley, J.M. 1989. Towards a Mechanistic Model for Multiphase Flow: The FourField Approach. Batelle, PNL (Draft Report). Otake, T., S. Tone, K. Nakao, andY. Mitsuhasi. 1977. Coalescence and Breakup of Bubbles in Liquids. Chemical Engineering Science 32:377383. Peebles, F.N., and H.J. Garber. 1953. Studies in the Motion of Gas Bubbles in Liquids. Chemical Engineering Progress 49:8897. Stuhmiller, J.H., R.E. Ferguson, and C.A. Meister. 1989. Numerical Simulation of B _ubbly Flow. Electric Power Research Institute NP6557. 61
PAGE 69
APPENDIX
PAGE 70
DlSCIII c C !IllS SDLVU 1111: EIIUATlOIIS FOil A MIXTURE aF A COIITlNUOUS C AIUI DISCRETE PIIAS(. lliE COMTliUIUS PHASE IS OM A FIXED IESII. C 1111: DlSCIITE PIIASE lS UQIIAII. 1111: COD IS A5 EliPlltlT AS C I'OSSIBLE VITIWT A SCliUC LIMIT. c c c C liCIIEX STRUCTIJII! A PlCTIJII! c c C VOLliES LST LSTPI LSTPZ c C .AJICTlOIIS LST LSTPI LSTPZ c C ac IAAL VOLliES AIUI .AJICTIOIISI ac I c c ... Cl...uT DATA c ... c CALL l...uT c c ................................................................................................................... .. t ... PRELI"liWI'f lNlllAI.lZATlOIIS c ... c IICYC I TilE. 1.1 lPRlNT I c c ........................................... ............ ... t L11ET1.1111 TO ClliiTlTIME C"fCU c .. c Ill CIINTlC IF INCYC.EQ.UCO TO 7'11 c: ... t "" IIERCE IIJIIBLES !1IA T TOUCH OR O'IDIUP c ............................................................................................................... .. CALL IIEIIGE c ... .. ............................................................................... ...... ...... C "" STATE PRIIPEIITlES c ................................................................................................................. .. c 791 CALL STATE IF IIICYC .EQ. II n DO 150 I 1, fill IFIIBUB!Il.ED.IJCO TO lSI 111111AUl \IOLPNUl . \IOLPUl 151 ClliiTIEMD IF t c c ................................................................................................................ .. C PARTlCl liAS$ EIIUATIOIIS SOLVED Flll YOLIIE C !ALSO filiAL SHAPE DlSTRIIUTlOIIS .uaJ VOIDS CALCULATED! c ... c CALL PM$$ t C If INC"II:.EQ;OICO TO c t ... IIEIIGE WUI.ES THAT TOUCH OR OVUIUP MOT FIRST TIME c .................................................................................................................. .. C CAI.L IIEIIGE COIITIC .......................................................... t ... OUTP\IT EDIT IF CALLED FOR c ... c CALL OUTP\IT IlliTE 16,90211 NC"II: DO 8Z9 N 1 .. TO BZ9 IIRlTEI6, 91Zllll,liPIIII, YP!IIl, VELPIIII &Z9 CIINTlC c ... '"' t ... LOGICAL IES:INNIIC Of 1111: 1111E STEP C'ICI.E c ..................................................................................................................... .. t c MCY'C NCYC 1 TillE TIME + DT XIII 1.1 If IIICYC .LT. I.cYI:l 111 ... 1.1 c c ... ... ............ .. t SET 1.1' OlD Til VALUES c: ........................................................................................................................... ... c: t C DO Sll I 1, LSTI'Z PRXl PRES!Xl 1111111( Xl All)( Xl VOIIIIIKI VOlDIU DO Zll I 1, lfllBUI!Il.EG.IlGO TO 211 IK,Il VOlOP!K,Il 11CTN!K,ll IIIICT 11,11 Zll ClliiTl311 CIIITlDO 411 J I, LS11'Z YELIIIJI YELIJI YEUFII!Jl YELSF!Jl 411 COIITlDO 501 I 1, ... IF!lBUI!ll.EQ.QlCO TO Sle 1111111. l(l) YELPN!ll YELPUI XPII!ll liP! ll VOLPNU I 1/DlP!ll 511 ClliiTJc: 63
PAGE 71
c ................... ... C !;lT IIITDPOUTED CCICTII&IOUS fiELD VAAIAILES FOil PARTICLE EQUS. e .... .. .............................. ... c DO Ill I 1 .. IF!IIUIUl. EQ.DlGO TO IDI RHOAVEIIl D D DO WO It l, LSTPl lllllAVEIIl RHOAVE!ll PERCTII,IlRIIllll 600 CIINTIMI c 700 c III:LAVEUl D.D .U:CAVEUl D.O liD 700 J Z, LSTPI VELAVE!l l III:LAVEII l ACCAVEIIl ACCAIIEIIl CUNTliiUE PEIICTJIJ,IlIIELS11Jl PERCTJIJ,Il.U:CEL!Jl 800 CCICTIIIUE c c c !;lT ALL G FORCES c :c CALL IIIUC c c ............................................................................................................................. .. C lllVAICf PARTICLE IIOIEliT1M EQUATICICS fOil VELOCITY AIID I'OSITIOII c ...................................................................... c CALL Ptllll c c ...................................................................................................................... .. C A/11/AIU CCICTINUOUS lDEIITIJI EQUATICICS FOR VELOCITY c .............................................................................................................................. .. c CALL QGI c c ......................................................................................................................... .. C on CCICTIIIDIS !lASS EQUATICICS SOLVED FOR PIIS$UI c .......................................................................................................................... .. c c CALL eliAS!! c ...................................................................................................................... ... C ... filiAL CCICTIIUIUS IIOIENTUII EIIIATIOII A/IVAIICEJ)IT FOR IIELOCITY c ............................................................................................................................ .. c c e ........................ ... .. c CCICTINUE Til CYCLE LOOP c ....................................................................................................................... ... c ;o TO ltl c c ................................................. e ... FOII!Jt&TS c ............................................ c 90ZG FCAI.AT(/,3X,'f
PAGE 72
RICIIT BIUIDAR'I PRSCLSTPll PR II'I.ICIT CORRECTION TO THE CONTIIUlUS PHASE VELOCITIES IDlE TD THE III'LIClT PRESSIJ! 111 c INCLUDECTTCOf!lllllll DO HI J z, LSTPI l J VELCJ) VELCJI + V11PCJllll'lliII'Ulll VELSFCJI VOID.IIJIVEL!Jl DO 100 1 1, .. IFCIIUICtl.EQ.OIGO TO Ill VELSFIJI VELSFCJI + VOIII'JCJ,tlVELPCII CDNTII&JE loCCLCJl CVELSFCJIVELSFli!Jil/DT VELSFNCJICVELSFNC.I+llVELSFNCJlll/C Z. oDXI c c c c Zll CONTII&JE IIUCDAR't COIIUITIDIIS LEJ'T IIIUITRY VELCZI VL VELUI VELCZI VELSFCl I VELSFC Zl RICIIT IICUIDARY VELCLSTPZI VELCLSTPII VELSFCLSTPZI VELSF
PAGE 73
c DO ZDD I 1, II' IF CALL DAAaZIWAlES,CSUID3,K,AY3l IFIAATID,GIE.I.99999)TIJI AFRI*TO.O FUll 1.0 FIZIKI 0.5 AHIIAYEIU CSUIDZ YliiZ AFRI*TZ FI31U 0.5 IIHOAYEIU CSUBDS VLIIEl AF!3 ELSE AFRI*T 11.0 AATIDl AFIIIIIIT FUll t.S RHO&VEfl.l CSUID Vl.JME AFRON'T FIZIK) 0.5 o AHIIAYEIU CSUIIDZ YliiZ AFRI*TZ FI31Kl 1.5 o AHDAYiiiKl CSUIDS VLIIEl AFRI*T3 c ElllliF ELSE CONTINI END IF SOZ CDITIC c RETIAl END SIIIROIITIIE IIIAGZI VEL VAlE ,CSUIDZ ,I,IIEYliOLDl C THIS SUIAIIIITII IS USED TD CALCUlATE THE DRAG cDUFIClE'N'T, CSUIDZ, C THAT A IIUllllE IIGULD SEE IF IT ENTIRELY IN Tt WAlE OF l LEAlie IIC IIUIBLE c c c REPII JZ. ORADUSPI I lOAISIVELPNI I>IIELVAlEl RHDAY!il II/VIS AEYIIIIL!JoiiEIIIll GlDABSI GAAV lVISoo4/IIIHOAIIEI I l SIItTEJIOSl TR&IUZ.D TRAHZ 4.0Z ;1(.214) TRANS1.1 GltO.H) TR5.47 o ClIO.ZS) CW ADD NEll DRAG COEFFICIENTS IF I REP Ill. LT, TAAHZl Tlll C XL OIII. OE XLOIIl.O IFIAEPUl .LT .XLOIIITIO 66
PAGE 74
ClZ4. CZZ4 GO TO 5'1'12 ENDIF ClZCt. O< 1.13/A[P(J) J C2. 7'P< 1 J .. C .61l 5'1'12 CSIJBDZ DIIUICCl,CZl ELSEIFIAEPIIJ.Ll. TIWito .lHil.REPI ll.c:E. TAANZ l TIIEH C 1, GZ75ClREPt 1 )ft CZO.ft&aGl( 0 ZSJaRP(J) CSIJIDZ DIIIMII Cl,CZl ELSE CSIJIDZI .Il. END IF AETURH EMD SUBROUTIII UI'UT c C THIS SUIAOUTIIE III'IITS .lHil PI!DaSSES ALL TIE PROBL DATA, C INITIAL CONDITIONS AND BDIJjDAIIY CDNDITIDNS c INCLUilf I TTCDIIIII l .... LS 0 C BASIC NUIEAICAL DATA C DT O.OOl DT o.n 01 0.05 C OX O.OlQ IUXCYC lD C MXCYC 1001 C taLLS 191 taLLS 1'1'1 C taLLS 491 PIPELEMDXFLOATIIICELLSI NPZS !NT I ICRAPifwl NOUT I NDEBUG I C INITIALIZE WilLE FLAS DO 57 I l,NP IWBIIll 57 CONTIIA c C PEAIHAMT IICII[)IES IEEDEII FOil LSTRUCTIIIE LST NCELLS 1 LSTPl LST + 1 LSntZ LST Z LSTM1 LST I c C EGUATIDII PAIWTEIIS PIE 3.14159 CW IWIE PIPE DIAIInER LARGIE EIICUGit MOT TO EFTECT BUBBlES W/ D6 01 C "ARA PIEl.O .. C P1PDIA 1.1 PIPDJA O.Z AREA PIEPIPDUZ/4.1 YOL ARE'ADX C CALC TIE IIAX EQUIY RADIUS FOil AIOYE PIPDIA W/ SPEHR1CAL CAPI'II C IUIIl[S .. llH ASPECT RATtO 4 0Z, CA9 VOL Cl/6liClAZI .. ZJ C RIWC O.SD+OOc PlPOlA/((S.D+00/4. 0ZD+00)+( Z C C3.0+00/lt 02D+Ot) + U .0+00/4. 0Z.O.OOlJ >U. D+OO/l.Dt00) RJWC 0 51>600 PJPDIA l lllllTE16,39ZlAEIWI 39Z FDRIIATI' AEIWI ',ElZ.S,' IIElDIS'I c SRAY 9.1 CS\1110 0.4 CADDI!S 0.5 CIMERA 1.0 XlEM 1.0 IIIOC'IC ZO C INITIALIZE SHAPE FACIDRS DO 94 M l,NP ISIW'EIHI 0 94 CONTIIA CW IMPUT WillE D1AIIElDIS tiTOSI C D1AlAEIIU.99'l'lll+OI DUl.005 DUZO.OlOO DI.UWO. 0050 DU4.0108 DIASo. oua D1A6a0.0101 DU7.01BO DUI.0101 DIA,.O.OIOO DU10.0100 DUUO. OUO DIAIZ.O. 0400 DIAll0.0101 DIA14 ,0100 DU130.ll00 DIA16al. 0601 DU17. 0101 DIAliO.OSOD. DU19wl.1001 DIAZOO. 0050 DIAZlO.OlOI DIAZZ.I.005 DIAZl1.0900 D1AZ4.0lOa DUZSwl.lOI WAITE !6,> 'BUB DIA 1 ' ,DIAl WRITE (6,) 'lUI DIA Z " ,DIAZ WAITE !6, 'BUB DIA 5 ' ,D1.U WAITE 16, 'BUB DU 4 ' ,DU4 WAITE 16, 'BUB DIA 5 ' ,DIAS WRITE (6,) 'SUI DIA 6 " 01&6 WRITE (6,) 'lUI DU 7 ' ,OU7 WRITE 16,) 'WI DIA I ' ,DIU WRITE t6,) 'lUI DIA 9 ' ,OU.9 WAITE 16,l 'BUB DIAIO ' DIAIO WRITE 16,) 'lUI DUll ' ,DlAU IIIIITE !6, 'BUB DU1Z ' ,DIAlZ WRITE (6,) 'lUI DIA13 ' ,01113 WRITE (6,J 'lUll DJ.&l(t ' ,DUlte 67
PAGE 75
WRITE t6,J WRITE (6, WRITE l6,) WRIT (6,J WRITE (6, WRITE t6,J WRITE t6,J WAITE (6,J WRITE l6,) WRITE (6,J WRITE t6,) WilE (6,) WRITE IIIRtTE t6,J WRITE (6, WAITE l6, WAITE (6, WRITE t6,J WRITE l6,) WRITE (6,a) WAITE l6, J WRITE l6,J WRITE (6,aJ WRITE C6,a) waiTE <6,J 'BUB DIAlS 'BUB DIAI6 ' ,DUI6 '808 DIAI7 ' ,DIAI7 'au1 DUll ' ,DIAlS '808 DIAI9 ' ,DIAI9 'BOB O!AZI ' ,DIAZO '808 D!AZI ' DUZI '808 DUZZ ' ,OUZZ '808 D!AZ3 ' ,01&23 8UI OUZ4 ' ,DUZ4 'lUI DIAZS ' ,OUZS P!Pt: LENCTH ' PIPt:LEM 'P!Pt: DIA P!PDIA 'AREA ', &aA 'YOL ', VOL ', ;a.w 'CSUBD ' CS\JID CAD liltS ' CADDIIS 'CIHEIIA ', CIIERA 'X1.EM ', XLEN 'UICYC ' UlftYC 'IWICYC ' IWICYC ox ox 'DT',OT 'llttLLS ', llttLLS c CW MATERIAL PIIOI'ERTIES C VISI.OOSE3 C SURTEM7 350EZ C fOil CL YCERIN C VlS. C SURTEMoO. 063 C FOil SUGAR WATER VIS O.Z6Z6 SURTEM0 071 CW ABOVE PROPS: VISIM5/IIIt"ZI,SIIITEMUCIIII AT ZD m:; C C STATE PAIWTERS PIIEF I. DES C WATER C RRF I .OEJ C GLYCERIN C RRF l .zsa C SUGAR WATER RRF 1.333 SOUND I. OE3 DROP 1.0/ISIUIJIOSIUCIII RPIIEF 1.0 S l.OEZ OWDP 1 0/(sotJNDPSOl.ldJPl WRITE t6,a) '$DUd) ", SOLIGJ WRITE ( 6, J SCIICDP SOlJGII IIRITE 16, I '.JIJICTICIC POSITICICS' DO ZOO J I, LSTPZ VELIJl 1.0 VELSF
PAGE 76
YI'C25l 0.12 IIIIITEC6,703l 783 FCIRMH' BUill X ANDY POSITJONS',/,SX."IIJ& t, 't1' I YP "I llO 733 I 1,11' VRITE<6 57lt U ,XI'( I) ,Yilt 1) 574 733 CDNTIIIJE llO 400 I I, II' IFCIBUICll.EQ.ou;:a TO 400 C XPCII 1.5 CW SET INITIAL IIUBBLE POSitiON TO 4TH .IJNCTIDN !FOil DXO.OIJ c XP
PAGE 77
ot: I(Kl AllllPPIKll XPPU en J/DNOtt XllQUl IDIKlAIKlXIICIIKlll/ot:1101 CONTINUE IIPILSTl IDILSTlAILSTlXIICIILSTII1l l/IIILSTlAILSTlXPPILSTIIll l DO lZOO ll t, LSllU1 c I LST I DPill XPPillIJPil+ll + Xlllllll Ull CONTINUE RliJIII END SlJaaouTII FIACT CXL,XP,XR,JCLDI,PERCT) c C THIS SUIROUTINE ;TS THE PERCENT IIF TliE Y!liD FOR TNE PARTICLE C LOCATED AT 7/Sl llfTIIEEII XL AND XA. A FDIJITH OIIDER SIIOCITH C IS USED. c c c nLitlT DIJIIILE PRCISIOM IAH, 0Zl XMIN liP O.aXLEN X1W1 XLEN XLL XL JOliN XRIIJIRXIIIM IIIJIPII IWIIIIIMIXRII,XLENl,D.DDDliXLEM XLIM:R "IMIIWIIXLL,D.DDOl,XLENliXLEN IAJI 6.1nPPRliZ lS.OnA' + 10.1 MLW XLIMRSIIC 6. IXUMRZ 15. DXLOWER + 10. D PRCT Ill. RTIIIII END SIIIIIOUTINE WTl'UT c C THIS SUIRIIUTINE DOES ALL THE OlJlliUT EDITS WliEll IIEILSTED c c lfllPIIliiT .liE. IC'/Cl (;0 10 Sll WRITE (6, 9010) MCYC, Ttl, DT, DX WRITE 16, 901Dl PRSHDILSTP1l PRS(LSTPll DO 50 J LSTPI, z, 1 K J PRSHDIIll PRSHDIIl CRAYHD!Jl 51 COMTINUE DO 100 I 1, LSTPZ, WRITE (6, 9DZCU l, VElCU, ..scu, VOlD(IJ ,PRESClJPIHDtiJ, PUSHDllJ, YELSfCIJ lDD CONTINUE IIRITE 16, 91511 DO zoo I 1, ._. IFIIIUIIll.EQ. Dl(;O TD ZDI ZDI IIRITE 16, 91211> I, liP II l VEL PI ll RAIJUSPII l, Y!li.PI ll, RHIIPAIIl, SC.&LEPI ll c cw CONTI !&I( IPIIINT IPIIINT NPIIINT IF IIC'IC .SE. IWICTCl STOP IIRITE 16,90181 9011 FORMATl/,lX,'MIII IUallE ItFO, CMCYC,XP,VELP,Il,MPJ') cw 500 CIIIITIU c C F.U.TS c 9011 fCJRMATC/ ,IX, 'ICYC' ,IS,SX, nfE' ,lPElZ.S,SX, "DT', 1PZ.5.sx. ox ,1Pz.:n 9010 FORMA Til ,IX, 'CELL JtFO' ,/,ZX, 'NODE' ,4., 'VELOCITY' ,r.x, 'MSSlJIE' ,ax, 'VOID' ,3JC, 'PRS..sHD" ,SX, 'PRSHD' ,ax, 'YELSf'' > 91ZD FOIIIIATI1X,l5,9UPE1Z.5ll 9031 FORIUTC/ ,zx, 'lUll' ,41, 'POSITICit' ,4X, 'VLOCITV' ,6X, 'RADIUS' ,6X, 'YOLliE' ,SX, 'DENSITY' ,6X, 'SCALEP' l RTIIIII EHO SUIIIOIITIJE c C THIS SIJIIIIlUT11II! SOLVES FO. THE PARTICLE Y!lLLK FliOII THE liASS C CONSERVATlOII EIIUATlOM. lT ALSO SElS TNE FINAL Y!l!DS AND IIEIGHTIICS C ASSOCIATED lnt THE PAIITICLES. c c c c c c DO lDI I 1, LSTPZ VOIDIU 1.1 111 COMTIU DO 150 J 1, LSTPl IIDIDJ!Jl' 1.1 151 COMTIU DO 411 I 1, ,. IFIIIUI!ll.EII.OIGO TO 400 VOLPIIl IIMIIPANIIJY!llPMIIl,_&lll 511 IIADUSPill 10.7SVOLPIIl1PIEIHil.DI5.DI 511 SCAI.EP( l> III.ENoDII DO ZOI I 1, LSTP1 J l CALL fRACT(XJ(Jl,XP!ll,XJIJOll ,SC.lLEPI ll,PRCT(l,lll VOIIJP(l,Il PRCTII,IlVOLPI II/VOL IFIVOIDPII,I l .GT .1.01>+00 IYOIDPU,l lO. 9'1'10 YIIIDII) VOIDill \IOIIIPil,ll ZDO COMTII&I( DO SOl J z, LSTP1 I J CALL FRACTIXIIll ,XPill,XIll ,SC.lLEPII l ,PRCT JIJ,Ill YOUJIIJCJ,U PERClJ(J,IJVOlP
PAGE 78
300 CIIIITIIIJE c 410 CIIIITIIIJE c RENIN SUIIIOIITII PIOI c C THIS SIIIIIOUIII SOLVES 11 PARTIClE LIHEAA IIOIIEHTU!I ECIUATIOHS FOil C THE VELOCITIES AHD THEN II'DATES IHE PAAIICLE POSITIONS. IT ALSO C G(IS II INTERFACE TEAIIS IN THE PARTICLE MCIOENIUII EIIIIATIOOIS fOil USE C IN Tl CONrtHUOUS PHASE IIOIIEHTIJI QIJATIONS. IT ALSO CALCULATES C nPARTlCLE DENSITY TEIIJI5 NUllED IN THE COIITIIIJOUS !lASS EQUATIDHS. c c c IHCLUilEITlCCIIDfl DO lDD J I, LSTPZ lfACEPIJl 0.0 ACCl:LPIJl 0,0 CAAVP(J) 0,1 !DO cONTIIIJE DO ZOO I 1, LSTPZ OAIIPIIl a.o ZDI COIITIIIJE 007001, .. IFFIZlll>FI311ll DToFPIIALLill VELPU l 5111/CDEF 903 fmtiiATl IM PMJM, Ft1t IUI8LE' ,13,1, I VELAVE,YP,VSO' ,3Z .Sl fll,Z,3' ,lElZ.S,J, c cc c c .. t t c c SOQ c &DVAIICI! POSITION CILIA Tlllll XPUl XPMUI O.S.ODTlVELPitlVELPIIIIll HPLOCJUl IDINTI !IPUll.SDXl/DX l l HPLDCVIll IDIHTI IXP!Il DX>/DX l l GET INTEGRATED Df PARTICLES FOil CIIIITIIIJOUS ECIUATllliiS DO 500 J z, LSTPl l J ACe !VELP!llVELAiUlliDT XFAI:e'!Jl XFACEP!Jl PERCTJIJ,ll< FllllIVELPUlVEUVEIIll ADDIIAS! ACCACCAVEUll l XFAI:e'IJl Xf&CO'(Jl I PERCTJIJ,Ill Fl!IlFPWALlliloVELPUI TEAll PERCTJ!J,llAHOPAIII IlYOLPIIII l ACCELP!JI lCCl:LP(JI TACC l&lVPIJI GRAVPIJI TEJIII 100 CONTIIIJE DO 500 I 1, ,. IFIIIUBUl.EQ.IlGO TO 300 AHOPAI 11 0 t DO ZOO It 1. LSTPI RMDPACU RHOP&Ul PERCTCI,JlRHOPCil 20t CONTIIIJE 301 CONTIIIJE 71
PAGE 79
RETI.RN [MD SUIIIOUTJIE W&XEll,RUIOI RUIOZ,IHiliXI C THIS SIJIIROUTIH[ CALCULATES MW DRAG CCt:FFICIHl !CSUliDI FOR BUIILES C THAT lR JHFLUEHCD IV TIE WAllE DF liiOTHEll BUBILE. c IHClUili!TTCOIIIOONI C SET AREA RATIO TO ZERO FOil EACH IIUIBLE c DO 603 l.l,NP IFIL.EQ.l.Oit.IBUBILI.ED.II ;Q TO 60S lRATIO SECOICD VAllE LENGTH SUP OUT C VAllE LDIGTH IIASEil ON STJI.LO PRCENT 0 .01 SOLU ( (O IZPCEHTIIlSIIIITI IO.lZOPCEHTIUZ I I II WAIIELS SOLU RAOUSP(Jl IUIIDIS XPIJI XP!ll IFIIUIIDIS.;T.V&XELSI THEN C WRITE16,705ll,J,WlllLS,IIUIDIS 705 FORIIATI' BUBBLE' ,13,' IS HOT EFFECTED BY IIU88LE' WAllE USED 1 Dll TIEIR IIERTICll DISPLACEIIIT : ', I, Z SECONDARY WAllE L.15Tll"' ,.5,' BUIIDIS ',.51 ;Q TO 100 Elllllf c ... .............. ... C 11011 CHECl If ITS IN TIE W&XE DUE TO HORIZONTAL POSITION C Cltl TO SEE IF TIE BUBBLE TOUOIES OR OVERLAPS TIE WAlE DF BUBBLE J C AT XP ( COIITIII.IE IF TIY 0011' T TOUCH) c C ClEO: TO SEE IF BUBBLE IJIIEIILAPS SECOMIIAIIY VAllE fliiST: CALL SECVAII(J,l,IIELSW,WAIIRDSl YYU YP!Il RAOUSPIIl IFIYYU. L T. O.OIYYUO,I YYS! YPCI) RAIJUSP(I.) YYSSo YP(Jl WAIIRIIS IFIYY55 .LT .I.OlYY5Sol.l YY66o YPI J l V.UCRIIS IFIYY66 .;T .PIPDUlYY66ofiiPIIJA JFIYYU .GT. YY66.DII. Y'lll. LT. YYSSITHEN APRDJI PIE IIAQUSP(ll Z ;o TD 101 EMDIF c C IF IIUI DOES OVERLAP SEC VUE (AJIN[), CHECl TO SEE IF 111111 IS IELIJII C THE ESTIIIAT VIDTII DF TIE PIIIIIAIIY VAllE xvs 7 1 XII RAIJUSP(Jl IFIIIU8DIS.5T .XV);Q TO TTT7 C IF IIUI CAll OVERLAP PIIIM Vlll AND SEC VAllE, CALC PIIIN VAllE RADIUS ALPHA DATAIIIRADUSPIJJ/1(11) TRILEN XV BUIIDIS WlllRDP TAILEN o DTANIALPMAl CH............ c C> PLACEI!EliT DF IIUIUE TO TIE VAlE AT XP!ll, c:> TO SEE IF BUBBLE C OVERLAPS THE VAllE ON TIE RIQIT 011 TIE LEFT DF THE VAllE OR IF THE C BUBBLE IS ENTIRELY VITIIIN TIE IIAIE C IIUBBLE IS LEFT DF VAllE C SUP OUT IF ENTIRELY IN WAllE OR PIIOJ VAllE AREA W/IN BUIIBLE APRDJI APROJa PIE RADUSPU. ) 2 C VIIITE16,Z4lll,J,VlllRIP,VllRDS 243 FCRIIATI' Fill IIUI' ,13,' IHFLISIC BY VAllE' ,13, 'IWlPIIIADS' ,ZZ.Sl YYZ YP!Jl WAIIRIIP IFIYYZ.LT O.OIY'IZO,I n VP(Jl + WAI.ADP IFIYY4,;T .PIPDUIYY4.PIPDU YYJ YPIK l AAOUSP!Kl IFIYYl .LT .O.OlYY1.0 YP!Kl IWIUSPill YY5 YP!Jl W.UCRDS IF!YYS.LT .O.Ol'IYSol.l Y'Y6 VP(J) + WAIRDS IFCV't6.;T .PIPDI&lYV6PIPII& &IIAIEZ PIE W&IIADS""Z AWAKE! PIE WAIUIJ)II4IoZ c C lliFIIE BUBBLE POSITION IN VAllE 111111 C 1 BUBBLE ENTIRELY IN SfCDIIDAIIY VAl[ Allll COII'I.ETL Y OVERLAPS PIIIIWI C Y Wll[ C Z IIUBBLE EHTIREL Y IN PIIIIWI'f IIAIIE C 3 IUBILE ENTIRELY IN SECONDARY WAKE BUT DOES NOT OIIOILAP PIIIIIAIIY C 4 IIUBILE INTERSECTS SfCOHDAIIY WAllE, IIUT NOT PIIINAIIY C 5 BUBBLE INTERSECTS BOTH PIIIIWIY AND SECDIIDAIIY WAXE C 6 BUaiLE INTERSECTS SECONDARY lHD OVERLAPS PIIUWIY WAll[ C 7 I!UaiLE ENTIRELY IN SECONDARY AND INTERSECTS PIIINARV C 8 SECOMDAIIY WAllE "ENCLOSED" Ill IIUBILE RADIUS c IFIYYS .LE. YYl.AND. YY1.LT. YVZ.AIIIl. YY4,lllll. 'Nl. LE. YY6ll1 72
PAGE 80
lffYVZ.LE. YVI.AMD. YY3. LE. VY4liWZ IFI'I'YS.LE. YYI.AHD. JYl. LE. VYZ liW IFIYV LE. VYI. AMD. JY3. LE. VY6 l 1S IF!'I'YS.L T YH.&MD. YVI.LT. \'Y5.ANII. YY3. LE. YYZl 1IFr'I'Y LE. VYI.AHD. VYI.L T. VY6.AMD. VY3.;T. YY6 ll.VYZ.AHD.YY3 .LT .VY AIID.VYI.LT .VYSll5 VYZ.AIID. YY1. L T. VY AIID. VY5.;T. YY6JIS IF!YYI.GT. VYS.AHD. YY1.L. YYZ.AIID. YY3.&T. YY6liWW IFIYY5.G. YY.AHD. YY5.L T. YY6. AND. YYI.LT. YYSl IW IFIYY1.1l. YYS,AIID. YYI.LT. VYZ.AIID. YY3.CT. VYZ.AND. YYS.LT YY4Jiw7 IFIYVl.CT. YYZ .AND. YYl. LT. YY4.AND. YY5.CT. YY4 .AND. YVS.LE. YY6Jiw7 IFIYYI.L. YYS.AND. YY6.LE. YY5JII IFIHCYC .E0.517 .OR.lCYC.EQ.SlllTIM WIIITEI6,15ZlNCYC ISZ fORIIATI' FOR CYCLE' ,131 ENlllf C SET ARU RATIOS aASED Ill IIIJIBLE POSITIIII IN Tl WUE!Sl c IFI IW.U.I lTl ARATID!Jl AWUU/oii'IIOA ARATIOZIJI IAPI!OJBAWUEII/.IPIIOA ZlTIM ARATIO!Jl 1.1 ARATIOZIJl 0.1 lw.EQ.51Tl ARATIO!Jl 1.1 ARATIOZIJI 1.0 ARATIOIJl 1.1 CALL INTER!l,J,YY5,YY6,AINTZI AIIATIOZI J l AllfTZIIPIIDJI ELSEIFI lw.EG.5lTlN CALL JMTER(It,J,YTZ,YY4,ADIT1) CALL INTER
PAGE 81
C IIRITE16,39
PAGE 82
fORIIAT(' BUilL' ,I], OUT Of' TOP r:1F PIPE AT TillE lSI' ,f& .41 llluaUJ a Gil TO 101 DCIIIF C CHEOI Fill SLUG FUN 8150 1111 IWC ECIUIV DUIElO STOP 409 1 z 3 IFIRACIUSPIIJ ,40, 11UlTIEII lEND I 1111Tl:i6,409 J NCYt,I ,YPill,liPUI,R&IIUSPUI Ill IIIITl:IZ a9JIICYt; ,I,YPI I),XPII I RIDUSPII I ,I SHAPE I I I FORIUT(" lllllllttllltlllllllllllllllltlllllllllllllllltllfl", 1," AT N'CYC'.ts, atJ8BLE',I3,"'RACHED SLUG flOW SIZ', t,' .,.. .. nz.s. JCPo' ,e1z.s. .nz.s, SHAPE' .n. ..................................................... ) Gil TO 541 ENIJIF C CHEOI !0 S'E IF IIUIIL MIIUPS Slot: Of' PIP IFUSHAP .EQ I l THlll RIGHTYPIIl R&IUSP!ll XLEFTYPUl IWlJSPtll ELSFI ISIW' I l ,Q,2)1'10 Ill 0 6Z9960>01 o Z. II> II o IWliSPU l Ill AZZ/Z. II> II RIQ!T YPU l + I2Z ILEFT YPIIl I2Z ELSE IF I ISIW' I l .3 lTIEII Ill IIZ.!I>OOORAIUSPI lll:VIU. IJOIOt' OZIJOOil+ Cl .O>OOt'.OZO>OID3lll H ti. IJOIO/l ,IJOOOJ AZZ. AZZ/Z .IJOOI IIQ!T YPU l 12Z XLEFT YPU l IU ELSE IIIITl:l.,59tliC'tC ,I S9t FDIIIIATI' UTTTTTT ISHII' IE l,Z,S AT NCVC',IS,' BUI',IS,'TTTT'l CDIITIIIE Ellllf 1Ft RIQHT. CT. PIPDIAIYPU lVPI ll( RIQ!T PIPDIAl IFIXLEFT .LT .I.IJYPUlYPUl+IJIISIXLEFTl 151 COIITIIIE c DO 111 ,_1,11' IFIIBUIIJI .EO. I)GII TD Ill lfiJ,fQ,IJGII TO 111 C CHEOI TO Mil: Sill IUIIIlE IS RELEAS IFIRELTIIIIJl. CT. TliElQII TO Ill c RSIII RACIUSPIJ) + IWliSPUl T IJIIS( JIP(JJ1!1'111 l TRIIZ DAliSI VP(JlYPill l Rt;AJ.t; IISIIRTI T RIIl .. Z + TERIIZHZ l lf!Rt;AJ.t; .CT RSIJI)QII TD 111 C CAI.Cili.ITl: lEV VELOCITY, IIRIUS !SIZE I oiiiD POSIT IIIII Of' lEV IUNI. 711 IOIASSl IIHOPIIIIll o VOLPNUJ IOIASSZ III
PAGE 83
IFIPRIIUC. LT.OLTUJGO TO AI.PHA DATAIIIRADIISPIIJIPRII!Xl TRIL.EN PAIMX DEL T.U AAOPRl TAILEN DT ANC ALPMA J C PRIIOI li.OIDECAYI D\DGIO$CII1(0 .01IWIIJSIIIIlZliR.lOUSPllll SID RHOREF AMOP.uoll l IF I VEL Pill li. L I .I.OElll0 CD Z4.0 GO TO 10 END IF CD ti. O RHOAFIIIHOAVEIIll DA!SlCRAVl a.o RliJUSPlll I 1 CVELPMCI>Z .S.OJ c IFICD .GT .Z4.0lCJloZ4 1 lD COHTIC CAI.C SECONDARY WAI VELOCITY IIASEll ON STIMllLER C UDUI 1.0 I ( O.Z U.IZot:LTAXIRADlJSPIIll + 1.01 o C 1 I Ot:L 1AXIR.llllJSPll ))ooz I C UDI UDUI VELPIIlll C UDIPRIIIXI./1 o.z U.lZoPRIIOV1WlUSI'Ull O.U C I IPIIIIIXIII.lDUSPIIll""Zl C UDPRIJOCI UDIPRIJIX VELPIIIll C CALC SECOICIIARY WAlE VELOCITY lASED ON CR.li1REE UDUI 16.2 R.lDUSPIIll I DEL TAX UDI UDUl VELPIIII I UDIPiiiiiX C6.Z R.lDUSP(J)) I PRIIII UIIIXI UDIPIIIIII VELPIIIIl C CALC SECIIICDARV WAlE VELOCITY IIASEll ON aHACA C 111011'1 DIVIDE IY ZERO ON FIRST TIIIESlEPl IFIVELPIII Il.LE.l.DEHDI UDZO.I UDPRIIOIZOO. I GO TO ZD END IF IVALIDAASI!;UV)o( Z. l IZ IIMDAVEII Ill IVELPIU llOVISl EVAL I. Ill. 7 C WRI1E16,7191BVAL,EVAL,DEL1U ,PRIIII,IWIUSPI II C 7&9 FORMAT( 1M S[CVME : a ,ElZ.S,' EVIL' ,ElZ.S,/ ,SX,' DELX' ,ElZ. S,. C I PIIIIOI' ,ElZ.S,' ,.51 C ECIUATION IS lASED ON DIS1AIIZ FROII IIIEAII S1AQIATION POIJfl, NIJII, TO C Tl NOSE OF THE FOLLDWIIC ..aLE DIFfER IDELTAX R.lDUSPIIll PRIIII IFIDIFFD.Ll.O.OlDIF'FER 1.1 C UllU lofVAI.l UDUZ BVALIIZ4. 1 ((!YAL/Z4.0J.7
PAGE 84
W VIII IDSIIIITI U.Go1 ID/0.6oGJ .. ll G.SOOGZl l.GoGJ ELSE )()( IXTN XVCliWAlL IFIQ.;T .G.40>0GllOIG.4D>U W VU (0SC:.f( tl.OD+I .. ClOC/0.4D+Ol .. Z) 0.51)1DU t.DD+I) END IF RETIIIII END SUBROUTINE SECI'OSIYI, YJ,IIII,V(LC ,DEL TAX ,PRIHX,RI, YELl ,PIPDUl C SUIIIIOUTIN TO CALC VELOCITY DF SEC WAitE FOR IIUULE IN VEL PROFILE IMPLICIT DOUILE PRECISION IAH,OZl IF I DEL T.U .LT .PRIIIXlTl YELl IVELC VELtODCPIl. ID>OOll I Z.G DID IF D DAIS IYI YJl POSl Y1 Rl POSR Yl Rl WAlL YJ Ill MAKR YJ W IIIL Ill 11 I IR I t IIAI 5Ift SEC WAlE DDES NOT DVERLAP PIPE IOIICD.UI\' IFIVAlL.LT .I.OlTHEII c VAlL 1.1 IL l DID IF IFI II.IU .liT .PIPDIAlTHEII WAU PIPDIA IR Ellllf C IF BUill aJ1TE11 IN SECII.&I T1tll VELOCITY IS T lUI COlER SEES IFIYI.IIT.V.utL.AICl.YI.LT.V.utRJTHEII IFIYI,LT.YJ.AICl.IL.EII.llllll YJ IFIYI.IIT.YJ.AIIIl.IR..11 PIPDU YJ WR DIIIIIUIIIL,IIIIIIl YEll VELC ODCPIIIl/IIIUNZl ElSE IFIPIISL. LE. YJ.III.POSII.IIE. YJlTHEII YEll IVElC VELCOOEICPIl.OoUl l/Z.I ELSE RDA&SI ... RU IFIYI.LT.YJ.AIIIl.IL.EII.lliiiL YJ IFIYt.;T.YJ.ANil.tR.ED.llllllll PIPDU YJ WR DICINIIWRL ,1111111 VElC OElCPIIRIIIRlZl ll1liF RETIIIII DID 77
