J Fluids Engineering 2009 Vol 131 N9
J Fluids Engineering 2009 Vol 131 N9
J Fluids Engineering 2009 Vol 131 N9
Fluids Engineering
Published Monthly by ASME
RESEARCH PAPERS
Flows in Complex Systems
091101
Pharyngeal Airflow Analysis in Obstructive Sleep Apnea Patients Preand Post-Maxillomandibular Advancement Surgery
John Huynh, Ki Beom Kim, and Mark McQuilling
091102
091103
091202
091203
091204
PUBLICATIONS COMMITTEE
Chair, B. RAVANI
OFFICERS OF THE ASME
President, AMOS E. HOLT
Executive Director, THOMAS G. LOUGHLIN
Treasurer, WILBUR MARNER
PUBLISHING STAFF
Managing Director, Publishing
P. DI VIETRO
Manager, Journals
C. MCATEER
Production Coordinator
A. HEWITT
Multiphase Flows
091301
091302
091303
091304
Downloaded 03 Jun 2010 to 171.66.16.159. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Contents continued
Journal of Fluids Engineering
SEPTEMBER 2009
TECHNICAL BRIEFS
094501
Analytical Solution for Newtonian Laminar Flow Through the Concave and Convex Ducts
M. Firouzi and S. H. Hashemabadi
Downloaded 03 Jun 2010 to 171.66.16.159. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
John Huynh
Resident
e-mail: jhuynh1@slu.edu
Ki Beom Kim
Assistant Professor
e-mail: kkim8@slu.edu
Graduate Orthodontic Program,
Center for Advanced Dental Education,
Saint Louis University,
St. Louis, MO 63104
Mark McQuilling1
Assistant Professor
Department of Aerospace and Mechanical
Engineering,
and Center for Fluids at All Scales,
Saint Louis University,
3450 Lindell Boulevard,
St. Louis, MO 63103
e-mail: mmcquil2@slu.edu
Introduction
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Computational Models
Sections 2.1 and 2.2 detail the method used to create and analyze CFD models starting from CBCT scans.
2.1 Sample and Imaging. Records of four patients who have
undergone combined orthodontic treatment and MMA surgery
were selected. The sample was comprised of one female and three
males. The mean age was 40.0 years. The inclusion criteria for
these patients included having a history of obstructive sleep apnea, being treated with a combination of orthodontic treatment
and orthognathic surgery, had MMA surgery, and had presurgical
CBCT scans taken within a week prior to surgery T0 and a scan
at least a minimum of 7 weeks following surgery T1. The exclusion criteria were patients with craniofacial syndromes and patients who had surgical expansion of their upper jaws. All surgeries were performed by one of two oral surgeons utilizing similar
surgical techniques. For the mandibular advancement surgeries, a
bilateral sagittal split osteotomy technique was utilized. A nonsegmental Le Fort I osteotomy was utilized for the maxillary advancements. Bone grafts were placed in areas where large surgical
movements were performed and rigid fixation was used to fixate
both jaws. Prior to fixation, the mandibular condyles were seated
in centric relation in all patients.
The CBCT scans of the patients were performed using the same
i-CAT CBCT machine at both pre- and postsurgical time points.
The field of view used was 23 19 cm2. The voxel size of the
CBCT scans taken was 0.4 mm. All scans were taken with the
condyles seated in centric relation. V-WORKS 3D software and DOLPHIN 3D were used view and analyze the CBCT scans.
The amount of surgical movement in the anteroposterior AP
direction was determined using the pre- and postsurgical CBCT
scans. The head position was oriented with Frankfort Horizontal
parallel to the horizontal plane x-axis. The facial midline was
then centered on the vertical plane z-axis. A plane perpendicular
to Frankfort Horizontal, passing through nasion, was created to
represent the y-axis. These three planes are displayed in Fig. 1.
Four skeletal landmarks basion, nasion, point A, and point B
were identified along the z-axis plane; these are shown in Fig. 2.
The AP surgical movement was calculated using the z-axis coordinate for each jaw. Basion was used as a reference point to evaluate the reliability of the position of the other landmarks. The rela091101-2 / Vol. 131, SEPTEMBER 2009
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
between 0.43 and 0.46 for all models. This h-ratio describes the
quality of surface triangulation, with the best h-ratio equal to 0.5
for an equilateral triangle. The tetrahedral volume mesh quality is
similarly described by its h-ratio, which is the ratio of radii of the
inscribed sphere and the circumscribed sphere of a tetrahedral
element. The tetrahedral h-ratio was between 0.26 and 0.27 for all
models in this study, where a value of 0.333 is the best possible
h-ratio. These numbers indicate the meshes used for the simulations described herein are of high quality.
In order to model turbulent flow within the pharyngeal airway,
the RNG k- turbulence model was used, as suggested by Sung et
al. 23 for simulations on similar geometries. Therefore, no other
turbulence closure models were compared. The inlet boundary
velocity was set as uniform and was determined for each model to
IMAGEJ
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 5
PRO/ENGINEER
correspond with the volumetric flow rates of 340 ml/s, 400 ml/s,
and 460 ml/s using the equation Q = AV. The inlet areas, A, were
computed from the solid modeling program. The outlet boundary
was set to a static pressure of 0 Pa for incompressible flow. A no
slip condition was set for the walls of the airway, which were also
assumed to be rigid and noncompliant. The simulations were performed on a PC with an AMD Athlon64 X2 Dual Core Processor 6400+ 3.21 GHz with 8 GB of RAM. Typical simulation
times were between 5 min and 8 min.
Results
others, where the MMA surgery creates more constriction in several locations of the airway by reducing the effective flow area,
represented here by the negative change in hydraulic diameter.
The other three patients experience an increase in hydraulic diameter throughout the airway. The authors suspect the local decreases in hydraulic diameter for patient 3 are the primary reason
for the difference in flow behavior when compared with the other
patients.
The Reynolds number is a dimensionless number used to characterize different types of fluid flow as laminar, transitional, or
turbulent. For laminar flows, the Reynolds number is 2000,
while turbulent flows have a Reynolds number 4000; in between
those two ranges lies the transitional phase of fluid flow 39. The
equation used to calculate the Reynolds number is Re= VD / ,
where is the density of air= 1.2 kg/ m3, V is the mean fluid
velocity, D is the hydraulic diameter, and is the fluid
viscosity= 1.8 105 N s / m2. The volumetric flow rate, Q, is
equal to Q = AV, where A is the cross-sectional area and V is the
mean fluid velocity 40. Substituting these quantities transforms
the Reynolds number equation into Re= 4Q / P. This equation
for the Reynolds number illustrates how for a constant volumetric
flow rate, the flow can become more laminar in nature due to
increasing the perimeter of a cross section in the airway, as done
with the advancement surgeries. Figures 912 show the variation
in Reynolds number throughout the pharyngeal airways for each
patient at T0 presurgery and T1 postsurgery for all three flow
rates. The Reynolds numbers observed for all patients at T0
showed there were significant areas of reduced perimeter, or
Fig. 6 Typical CADTHRU models used for geometry cleanup. a Model after importing into
surfaces removed in cleanup process.
CADTHRU.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Maxillary
Mandibular
Pt 1
Pt 2
Pt 3
Pt 4
Mean
6.1
5.8
4.6
8.5
3.3
5.3
2.0
11.3
4.0
7.7
12
10
2
0
10
20
30
40
Distance from inlet [mm]
50
60
70
4500
4000
Reynolds number
3500
3000
2500
2000
1500
1000
500
10
20
30
40
Distance from inlet [mm]
50
60
70
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
4500
4000
Reynolds number
3500
3000
2500
2000
1500
1000
500
10
20
30
40
Distance from inlet [mm]
50
60
70
4500
4000
Reynolds Number
3500
3000
2500
2000
1500
1000
500
10
20
30
40
Distance from inlet [mm]
50
60
70
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
4500
4000
Reynolds Number
3500
3000
2500
2000
1500
1000
500
10
20
30
40
Distance from inlet [mm]
50
60
70
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Pt
Pt
Pt
Pt
Pt
Pt
Pt
Pt
1,
1,
2,
2,
3,
3,
4,
4,
T0
T1
T0
T1
T0
T1
T0
T1
340 ml/s
400 ml/s
460 ml/s
12.7
1.1
23.6
2.3
2.1
2.4
43.5
1.1
17.1
1.5
31.9
3.1
2.8
3.2
59.3
1.4
22.3
1.9
41.4
4.0
3.5
4.1
77.7
1.8
Conclusions
In this study, CFD modeling illustrated the effects maxillomandibular advancement surgery had on flow characteristics in the
pharyngeal airway of obstructive sleep apnea patients. This paper
Table 3 Percentage change in airway resistance
Pt
Pt
Pt
Pt
1
2
3
4
340 ml/s
400 ml/s
460 ml/s
91.2
91.2
+11.1
97.5
91.3
91.3
+12.5
97.6
91.4
91.4
+13.7
97.7
Acknowledgment
The authors wish to thank George Huang of Wright State University and David Welsh of Software Cradle Co. for their guidance with the computational modeling.
Nomenclature
A
D
m
P
p
Q
R
Re
T0
T1
V
cross-sectional area m2
hydraulic diameter, D = 4A / P m
= Q kg/s
mass flow rate, m
perimeter m
total area-averaged pressure drop Pa
flow rate, Q = VA m3 / s
Pa s / kg
flow resistance, R = p / m
Reynolds number, Re= VD /
time denoting one week prior to surgery
time denoting seven weeks after surgery
mean flow velocity m/s
Greek Symbols
density of air kg/ m3
viscosity of air kg/ m s
References
1 Young, T., Palta, M., Dempsey, J., Skatrud, J., Weber, S., and Badr, S., 1993,
The Occurrence of Sleep-Disordered Breathing Among Middle-Aged
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
Advancement Surgery Utilizing Cone Beam Computed Tomography, MS thesis, Saint Louis University, Saint Louis, MO.
Sung, S., Jeong, S., Yu, Y., Hwang, C., and Pae, E., 2006, Customized ThreeDimensional Computational Fluid Dynamics Simulation of the Upper Airway
of Obstructive Sleep Apnea, Angle Orthod., 76, pp. 791799.
Shome, B., Wang, L., Santare, M., Prasad, A., Szeri, A., and Roberts, D., 1998,
Modeling of Airflow in the Pharynx With Application to Sleep Apnea,
ASME J. Biomech. Eng., 120, pp. 416422.
De Backer, J., Vanderveken, O., Vos, W., Devolder, A., Verhulst, S., and Verbraecken, J., 2007, Functional Imaging Using Computational Fluid Dynamics
to Predict Treatment Success of Mandibular Advancement Devices in SleepDisordered Breathing, J. Biomech., 40, pp. 37083714.
Xu, C., Sin, S., McDonough, J., Udupa, J., Guez, A., and Arens, R., 2006,
Computational Fluid Dynamics Modeling of the Upper Airway of Children
With Obstructive Sleep Apnea Syndrome in Steady Flow, J. Biomech., 39,
pp. 20432054.
Vos, W., Backer, J. D., Devolder, A., Vanderveken, O., Verhulst, S., and Salgado, R., 2007, Correlation Between Severity of Sleep Apnea and Upper
Airway Morphology Based on Advanced Anatomical and Functional Imaging, J. Biomech., 40, pp. 22072213.
Croce, C., Fodil, R., Durand, M., Sbirlea-Apiou, G., Caillibotte, G., and Papon, J., 2006, In Vitro Experiments and Numerical Simulations of Airflow in
Realistic Nasal Airway Geometry, Ann. Biomed. Eng., 34, pp. 9971007.
Wen, J., Inthavong, K., Tu, J., and Wang, S., 2008, Numerical Simulations for
Detailed Airflow Dynamics in a Human Nasal Cavity, Respir. Physiol. Neurobiol., 161, pp. 125135.
Heenan, A., Matida, E., Pollard, A., and Finlay, W., 2003, Experimental
Measurements and Computational Modeling of the Flow Field in an Idealized
Human Oropharynx, Exp. Fluids, 35, pp. 7084.
Huang, Y., White, D., and Malhotra, A., 2005, The Impact of Anatomic
Manipulations on Pharyngeal Collapse: Results From a Computational Model
of the Normal Human Upper Airway, Chest, 128, pp. 13241330.
Huang, Y., White, D., and Malhotra, A., 2007, Use of Computational Modeling to Predict Responses to Upper Airway Surgery in Obstructive Sleep
Apnea, Laryngoscope, 117, pp. 648653.
Jeong, S., Kim, W., and Sung, S., 2007, Numerical Investigation on the Flow
Characteristics and Aerodynamic Force of the Upper Airway of Patient With
Obstructive Sleep Apnea Using Computational Fluid Dynamics, Med. Eng.
Phys., 29, pp. 637651.
Ishikawa, S., Nakayama, T., Watanabe, M., and Matsuzawa, T., 2006, Visualization of Flow Resistance in Physiological Nasal Respiration: Analysis of
Velocity and Vorticities Using Numerical Simulation, Arch. Otolaryngol.
Head Neck Surg., 132, pp. 12031209.
Wang, K., Denney, T., Morrison, E., and Vodyanoy, V., 2005, Numerical
Simulation of Air Flow in the Human Nasal Cavity, Conf. Proc. IEEE Eng.
Med. Biol. Soc., 6, pp. 56075610.
Patanker, S., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere,
New York.
van Leer, B., 1977, Towards the Ultimate Conservation Difference Scheme 4:
A New Approach to Numerical Convection, J. Comput. Phys., 23, pp. 276
299.
Rhie, C., and Chow, W., 1983, Numerical Study of the Turbulent Flow Past
an Airfoil With Trailing Edge Separation, AIAA J., 21, pp. 15251532.
White, F., 2006, Viscous Fluid Flow, McGraw-Hill, Singapore.
Fox, R., McDonald, A., and Pritchard, P., 2006, Introduction to Fluid Mechanics, Wiley, Hoboken, NJ.
Collins, J., Shapiro, A., Kimmel, E., and Kamm, R., 1993, The Steady Expiratory Pressure-Flow Relation in a Model Pulmonary Bifurcation, ASME J.
Biomech. Eng., 115, pp. 299305.
Isono, S., Feroah, T., Hajduk, E., Brant, R., Whitelaw, W., and Remmers, J.,
1997, Interaction of Cross-Sectional Area, Driving Pressure, and Airflow of
Passive Velopharynx, J. Appl. Physiol., 83, pp. 851859.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
C. L. Chen
Teledyne Scientific Company,
1049 Camino Dos Rios,
MC A10,
Thousand Oaks, CA 91360
Laminar air flow through a curved rectangular channel with a variable cross-sectional
(c/s) area (diverging-converging channel) is computationally investigated. Such a flow
passage is formed between the two fin walls of a 90 deg bend curved fin heat sink, used
in avionics cooling. Simulations are carried out for two different configurations: (a) a
curved channel with long, straight, constant c/s area inlet and outlet sections (entry and
exit lengths); and (b) a short, curved channel with no entry and exit lengths. Formation
of a complex 3D flow pattern and its evolution in space is studied through numerical flow
visualization. Results show that a secondary motion sets in the radial direction of the
curved section, which in combination with the axial (bulk) flow leads to the formation of
a base vortex. In addition, under certain circumstances the axial and secondary flow
separate from multiple locations on the channel walls, creating Dean vortices and separation bubbles. Velocity above which the Dean vortices appear is cast in dimensionless
form as the critical Dean number, which is calculated to be 129. Investigation of the
friction factor reveals that pressure drop in the channel is governed by both the curvature
effect as well as the area expansion effect. For a short curved channel where area
expansion effect dominates, pressure drop for developing flow can be even less than that
of a straight channel. A comparison with the flow in a constant c/s area, curved channel
shows that the variable c/s area channel geometry leads to a lower critical Dean number
and friction factor. DOI: 10.1115/1.3176970
Introduction
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
recirculation cell as shown in Fig. 2c. For higher velocity conditions, in addition to the main secondary flow, multiple secondary flow structures are formed, commonly known as Dean vortices, after the first analysis by Dean 4. The nondimensional
number characterizing the onset of Dean vortices is called the
critical Dean number Decrit.
The curved channel/duct flow is widely encountered in many
engineering applications, such as bent tube heat exchangers, turbine blades, air conditioning ducts, the condenser area of solar
collector heat pipes, passages in biological systems, etc. A significant amount of prior research has investigated the flow and thermal fields of a curved channel. Since this article focuses on the
flow characteristics, we only cite the flow related literature for the
sake of brevity. Cheng et al. 5 experimentally visualized the
secondary flow pattern in curved rectangular channels, 25 mm
wide b, curvature ratio CR of 5, and the aspect ratio AR
varying between 1 and 12. CR is defined as the ratio of radius of
curvature R and the channel width b, while AR is the ratio of
channel height H and width b. Sugiyama et al. 6 reported
flow visualization at the end of 180 deg turn in channels, 20 mm
wide, CR varying between 5 and 8, and AR varying from 0.5 to
2.5. Ligrani and Niver 7 experimented with a nearly 2D channel
AR = 40 for De varying between 40 and 220, and showed the
formation of secondary flow patterns at various angular positions
along the curved section. Effort on the numerical front dates back
to the early 80s. Cheng and Akiyama 8 reported a 2D model of
fully developed laminar flow entering a curved rectangular channel, while Ghia and Sohkey 9 and Humphrey et al. 10 investigated the developing flow. In these calculations, AR varied in the
range of 0.25, and De was limited to 100. Hwang and Chao
11 investigated fully developed laminar flow in a curved square
duct under isothermal condition. More recently, Chandratilleke
and Nursubyakto 12 presented a 3D computational model of
flow through a curved rectangular channel with uniform heating at
the outer concave wall. They considered parametric ranges of
1 AR 8 and 20 De 500, and included the effect of buoyancy. To summarize the findings of these research worksThe
Decrit value; the main secondary flow and the Dean vortices for
De Decrit condition all depend on several geometric and thermal
conditions, such as the channels, AR and CR location along the
curved section and the wall heating pattern.
In spite of 3 decades of research, most prior works report flow
pattern at a specific cross section. Literature on the evolution of
multiple vortices in a curved channel remains scarce. In an experi091102-2 / Vol. 131, SEPTEMBER 2009
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
the inlet section b is used for calculating all the geometric parameters such as hydraulic diameter Dh, AR, and CR.
For both the long and short channels, air at temperature Ti
enters the channel at a uniform flow velocity of Vi. Ambient pressure boundary condition is specified at the channel exit. The channel bottom and the side walls are maintained at a constant temperature of Tw Tw Ti. The isothermal side wall assumption
essentially means that the fins of the heat sink shown in Fig. 1
are 100% efficient. Indeed a computational model of the entire
heat sink made of aluminum shows that the fins are close to 99%
efficient. In a practical avionics application, often the top of the
fins is covered by a plate to avoid flow bypass. To simulate this
condition the top channel wall is assumed adiabatic.
The two critical dimensionless flow parameters are the Reynolds number Re= ViDh / and the Dean number De
= Re Dh / R. Vi values are chosen on the basis of a practical
range of air flow rates used in the heat sink. Air properties at
average temperature Tav = Ti + Tw / 2 are used to report these
numbers. Re and De vary between the range of 183 to 2017 and
92 to 1008, respectively. The last air flow passage, having a length
Lav / Dh = Lcurv / Dh of 6.8, represents one asymptotic condition
where uniform flow directly enters the curved section. Therefore,
it is picked as the short channel in this study. However, the other
asymptotic condition of fully developed flow entering the curved
section does not occur for any of the channels of this heat sink.
The dimensionless entry length Lfd / Dh for fully developed flow
in a constant c/s area straight channel of AR close to 2 is
0.02 Re 13. At the highest velocity condition of Re= 2017,
Lfd / Dh is 40. Not even the longest air flow passage of the heat
sink, shown in Fig. 1, has this much entry length. In other words
for Re= 2017, Lentry Lfd for all channels. Instead of focusing on
the ideal asymptotic condition, we use the current heat sink geometry and arbitrarily pick a long channel having an overall
average length of Lav / Dh = 39.1 Lav = Lentry + Lcurv + Lexit, with
Lentry / Dh = 16.1 and Lexit / Dh = 16.1. Flow entering the curved section of the chosen long channel is fully developed for Re 805.
Above this value, the flow enters the curved section in partially
developed state. Clearly if a different channel were chosen as the
long channel, magnitudes of the flow and thermal parameters
would be different. However, since this study is not intended to be
a parametric study of the geometric features, such as Lentry / Dh,
Lcurv / Dh, Lexit / Dh, CR, AR, etc., the specific values are not so
important. The primary objective here is to understand the complex flow field for a given geometry. Therefore, for both the long
and short channels, the channel geometry is kept constant, with
Dh = 0.007 m, bmax / b = 1.5, AR = 2.3, and CR = 5.6; and the
Lentry / Dh, Lcurv / Dh, and Lexit / Dh values discussed above.
3
2
Physical Model
Numerical Simulation
Steady, laminar air flow through the passages are computationally simulated using the commercial computational fluid dynamics
CFD code FLUENT. The continuity, momentum, and energy
equations are solved. These equations are coupled due to the use
of temperature dependent air properties. Some of the underlying
assumptions and observations are as follows.
1 Even though all the features of solving compressible flow
are considered, the flow is essentially incompressible. The
Mach number Ma based on the highest inlet velocity Vi
is 0.015.
2 The value of Gr/ Re2 is 9 107 for the lowest flow velocity condition, where Gr is the Grashof number. Therefore, buoyancy effects are neglected. It may be argued that
the local velocity at any point in the channel could be significantly lower than Vi. Therefore, if the above ratio is
calculated with Re value based on the local velocity, the
Gr/ Re2 ratio can actually become orders of magnitude
higher. To address this concern, we carried out a calculation
SEPTEMBER 2009, Vol. 131 / 091102-3
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
carried out a flow visualization study at the end of the bend and
reported the critical Dean numbers Decrit for the onset of Dean
vortices for various channel aspect ratio AR = H / b and curvature
ratio CR = R / b values. Figure 4b shows a comparison of the
predicted values from our computational model with the experimental data, for a constant CR value of 5. Note that this CR value
is very close to the variable c/s area curved channel geometry
investigated in this article CR = 5.6. This experimental paper calculated De based on channel width b, i.e., De= Re b / R, even
though Re is calculated based on Dh Re= ViDh / . We report De
based on Dh, i.e., De= Re Dh / R. However, for the sake of the
present comparison, just for this section Fig. 4b, we follow
their procedure. An excellent agreement between the prediction
and experimental data in Fig. 4b, for the entire range of AR
values, lends significant confidence about the present computational model.
Fig. 4 Validation of the numerical model against the experimental data of Sugiyama et al. 6. a Schematic of the experimental channel geometryconstant cross-sectional area
curved channel with a 180 deg bend. Lentry = 1000 mm, b
= 20 mm, CR = R / b = 5, and H and AR=H / b variable; b comparison of numerical prediction of the critical Dean number
Decrit with the experimental data 6.
The accuracy of the numerical procedure is ascertained by validating the model against existing experimental data from literature. For this purpose, a computational model is developed with
channel geometry replicating the experimental set up of Sugiyama
et al. 6. Isothermal air flow through a constant c/s area 180 deg
bend, as shown in Fig. 4a, is considered. In the model, a uniform
velocity is prescribed at the inlet of the straight entry section
Lentry, whereas an ambient pressure outlet boundary condition is
specified at the end of the 180 deg bend. Sugiyama et al. 6
091102-4 / Vol. 131, SEPTEMBER 2009
A grid independent numerical solution has to resolve the flowthermal characteristics of the overall channel as well as all the
secondary recirculation cells, at the highest flow rate condition.
The following parameters are used to quantify the secondary flow
structures at two angular positions, 45 deg and 90 deg: a the
total number of cells Ncells at that c/s area, b minimum helicity
Hmin, and c the dimensionless maximum velocity at that c/s
e.g., cross-sectional Vmax at 45 / Vi. In addition, two points A and
B are chosen, respectively, at the 45 deg and 90 deg angular
positions located at the core of the recirculating flow region. The
local dimensionless velocity such as VA / Vi, and temperature
such as A = TA Ti / Tw Ti are investigated at these points.
Four other dimensionless parameters are used to characterize the
overall channel transport fields. They are: the maximum velocity
at the exit Vmax,outlet / Vi of the long channel, the overall channel
pressure drop 2P / V2i , the average fluid exit temperature
exit = Tav,exit Ti / Tw Ti, and the Nusselt number Nuo
= hoDh / k f , characterizing the total heat transfer from the channel.
Values of these local and global parameters, obtained from computational simulations with four different grid sizes, are presented
in Table 1. Based on the comparison, a 351,000 cell grid is chosen
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Table 1 Grid independence study for a long VACC for the highest flow velocity condition. Re= 2017 and De= 1008 based on
uniform flow velocity Vi and hydraulic diameter Dh at channel inlet and properties at average temperature Tav = Ti + Tw / 2.
AR = 2.3, CR = 5.6, Lav / Dh = 39.1, Lentry / Dh = 16.1, and Lexit / Dh = 16.1. Location of point A: axial locationon the 45 deg angular position from the start of the curvature, heightat H/2, and lateralat b45 deg / 4 closer to the ICW. Location of point B: axial
locationon the 90 deg angular position; height and lateral locationat the geometric center of the curvature outlet plane.
Grid size No. of cells
Flow-thermal characteristics
126,720
230,400
351,000
620,880
2
44,027
1.42
0.17
0.78
2
51,280
1.45
0.08
0.79
6
51,686
1.45
0.10
0.88
6
52,079
1.45
0.09
0.82
6
24,775
1.62
1.30
0.31
6
28,313
1.64
1.35
0.32
6
28,610
1.64
1.35
0.32
6
28,776
1.65
1.37
0.31
At channel outlet
1.77
0.61
1.82
0.62
1.85
0.62
1.85
0.62
2.97
5.90
2.99
6.03
3.03
6.09
3.04
6.10
As the channel geometry plays a pivotal role in the flow dynamics, it is useful to decouple the effect of variable c/s area first.
Therefore, we start with investigating flow through a 90 deg bend
of a long CACC with entry and exit lengths. The channel geometry matches the schematic in Fig. 3, except for the fact that channel width b remains constant through the curved section. These
results serve as the baseline case for comparison with those of
variable c/s area curved channel. Even though flow through
CACC has been widely investigated in the past, studies focusing
on the space evolution of the flow structures remain scarce.
To follow the complex 3D secondary flow pattern, we investigate the cross-sectional view of the secondary flow, i.e., the inplane velocity vectors, at various angular positions Figs.
5a5e, in conjunction with the oil flow path lines on the ICW,
OCW, and TW Figs. 6a6c, respectively. The in-plane velocity vector plots are color coded by dimensionless local velocity
magnitude =local velocity V/uniform channel inlet velocity
Vi. Since the flow is symmetric heightwise, in Fig. 5 the computational results are shown in only half the c/s, while the other
half is used for explanatory schematics. The oil flow lines are
plotted by releasing particles from the same curved section walls
ICW, OCW, or TW on which the flow is visualized. The lines
are color coded with particle ID. Thus the color codes do not have
any physical significance as such. To retain the clarity of the oil
flow line plots, various markings are shown on a separate, identical plot below each figure.
The secondary flow sets in immediately after the fluid elements
enter the curved section. This secondary flow along with the forward moving axial flow, i.e., combination A, creates a pair of base
vortices, at the top and bottom portion of the channel. The traces
of base vortices in the oil flow lines plot Fig. 6 indicate that they
travel through the entire length of the curved section 090 deg
angular position, although their shape and size change along the
path. At the early portion of the curved section, the secondary
flow is fully conformal to the channel walls, as in Fig. 2a. ThereJournal of Fluids Engineering
fore a pair of base vortices covers the entire c/s area of the channel, as seen in Fig. 5a for an angular position of 11 deg. Flow
from the upper and lower base vortices collide and separate from
the ICW and attach to the OCW, both around midheight, respectively, forming separation line S1 in Fig. 6a, and attachment line
A1 in Fig. 6b.
Further downstream, the secondary flow starts to separate from
the walls, in the corner regions between ICW, and the top and
bottom walls. As the inset of Fig. 5b shows, the secondary flow
separates from the TW and attaches to the ICW, forming a local
recirculation similar to that shown in the schematic of Fig. 2c.
The oil flow lines in Fig. 6c indicate that the separation from
TW starts at an angular position P1, just before the 22.5 deg
location. This local, corner region, secondary flow recirculation
can be traced along separation line S2 on TW Fig. 6c and
attachment line A2 on ICW Fig. 6a, both continuing all the way
to the end of the curved section. This secondary flow reversal,
along with the forward moving axial flow combination B leads
to the formation of another local vortex. Since this is in addition
to the base vortex, following the traditional definition, it is categorized as a Dean vortex. We name it the corner region Dean vortex.
The oil flow lines in the inset of Fig. 6a show another interesting aspect of the corner region flow. Between the angular positions P2 and P3, the oil flow lines are backward and upwards.
The backward axial flow can also be seen from the insets B and C
of Fig. 6c. This indicates an axial as well as secondary flow
separation in the region, or combination D, that leads to the formation of a closed separation bubble. Figure 7a shows a schematic of this corner region axial and secondary flow separation,
and the formation of the separation bubble. The corner region
Dean vortex bypasses this separation region by detaching from the
ICW, squeezing toward the separation line S2, and then coming
back to the ICW beyond the separation region.
After attaching to the ICW along attachment line A2, the secondary flow follows the wall contour for a while, flowing toward
the centerline heightwise. However, before reaching the centerline it bounces off the ICW, causing another local secondary flow
separation and flow reversal, as sketched in Fig. 2b. The trace of
this separation line S3 in Fig. 6a, and Fig. 5c indicates that this
flow separation starts just upstream of the 34 deg position and
increases in size along the curved path follow Figs. 5c5e.
SEPTEMBER 2009, Vol. 131 / 091102-5
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(c)
(b)
(d)
(e)
Fig. 5 Formation of secondary flow at various angular positions along a 90 deg bend of a long CACC. Re= 2017 and De= 1008
based on uniform flow velocity Vi and hydraulic diameter Dh at channel inlet, and properties at average temperature Tav
= Ti + Tw / 2. Lav / Dh = 39.1, Lentry / Dh = 16.1, Lexit / Dh = 16.1, AR = 2.3, and CR = 5.6. Views: a angular position 11 deg, b 22.5 deg, c
37 deg, d 56 deg, and e 90 deg outlet from the curved portion of the channel or beginning of the straight constant c/s area exit
channel.
The angular position where this separation started cannot be precisely located from the current data. We tentatively identify an
angular position P4. Coupled with the forward moving axial flow
this secondary flow separation Combination B leads to the formation of another Dean vortex, which we term as the ICW Dean
vortex. Flow from the top and bottom ICW Dean vortices meet at
the heightwise centerline of ICW. Thus, beyond the angular position P4, the ICW centerline becomes an attachment line A3 Fig.
6a.
A comparison of Figs. 5a5e shows that as the flow moves
downstream, the return secondary flow from ICW back toward
OCW gets stronger. As a result the location of the highest velocity in the channel shifts toward the OCW and the velocity boundary layer there gets thinner. An abrupt drop in the centrifugal force
occurs across this thin boundary layer near the OCW. At some
location between the 45 deg and 56 deg angular position, it causes
091102-6 / Vol. 131, SEPTEMBER 2009
another separation of the secondary flow, this time from the OCW.
A local flow reversal follows, as seen in Figs. 5d and 5e. This
along with the forward moving axial flow combination B creates
the third pair of Dean vortices in the channel, now termed as the
OCW Dean vortex. It is bounded by the separation line S4 and
attachment line A4 in the oil flow lines of Fig. 6b.
Figures 5d and 5e show an additional pair of recirculation
cells inside the c/s, adjacent to the base vortex. Following the
conventional definition, this should be another Dean vortex. The
origin of this vortex is not clear. It is possible that this vortex
arises from a local split of the base vortex, somewhere downstream of the 45 deg angular position. This Dean vortex weakens
shrinks in size and gradually shifts toward the center of the
channel in the downstream direction.
To summarize, Fig. 7b shows a schematic of the base vortex
and all the Dean vortices in the top half of the channel, close to
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 6 Oil flow lines on the walls of the curved section 90 deg bend of a long CACC. Re= 2017 and De= 1008, Lav / Dh = 39.1,
Lentry / Dh = 16.1, Lexit / Dh = 16.1, AR = 2.3, and CR = 5.6. a ICW or concave wall, b OCW or convex wall, and c TW. All views are
looking into the wall from front.
the curvature exit. The origin of the multiple vortices and their
evolution in space are also summarized in Table 2.
used for the CACC, there was no asymmetry. Clearly the temperature dependent property effect is more pronounced for the decelerating flow in the diverging section. Since our focus is on understanding the overall flow field, for all practical purposes we will
ignore the slight asymmetry in the subsequent discussions.
The main/base secondary flow at the early part of the channel
close to the curvature entry is very similar as before. Fully conformal secondary flow, in conjunction with the forward moving
axial flow combination A, creates a pair of base vortices at the
top and bottom half of the channel Fig. 8a. They meet at separation line S1 on ICW Fig. 9a and attachment line A1 on OCW
Fig. 9b. As before, the secondary flow separates from the corner regions between the ICW and the top and bottom walls somewhere between the 11 deg and 22.5 deg angular position Point P1
SEPTEMBER 2009, Vol. 131 / 091102-7
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 7 Schematics of various 3D flow patterns in the curved section 90 deg bend of a long
CACC. Re= 2017 and De= 1008, Lav / Dh = 39.1, Lentry / Dh = 16.1, Lexit / Dh = 16.1, AR = 2.3, and
CR = 5.6. a Axial flow separation from ICW near the corner region and formation of a separation bubble, and b base vortex and other Dean vortices, such as corner region Dean
vortex, ICW Dean vortex, and OCW Dean vortex.
in Fig. 9c and follows separation line S2 all the way to the end
of the curved section. Between positions P2 and P3 inset A of
Figs. 9a and 9c, the axial flow also separates. The overall flow
field in the corner region consists of the corner region Dean vortex
and a separation bubble, similar to that sketched in Fig. 7a.
Flow in the diverging section of the VACC, i.e., between the
angular positions 045 deg, is significantly different from that of
the CACC. The first noticeable effect of the variable c/s area is
near the OCW, just after the flow enters the curved section. The
inset of Fig. 9b shows that at the top half of the channel some of
the oil flow lines are upwards and backward. The upward flow is
indicative of conformal secondary flow at OCW, while the backward flow is due to axial flow separation. The separated axial flow
and the conformal secondary flow lead to combination C, which
leads to the formation of a base vortex, but of a different structure.
Figure 10a shows a sketch of this base vortex structure. Figure
9b shows that the axial flow separation initially starts near the
top and bottom walls, and gradually expands toward the center
angular position P4, along the separation line S3. From the angular position of P5 on the centerline, the axial flow starts to
reattach to the wall along attachment line A2. The reattachment is
triggered by separation of the axial flow from the ICW discussed
later, which pushes the fluid back toward the OCW.
091102-8 / Vol. 131, SEPTEMBER 2009
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Table 2 Summary of flow pattern along the 90 deg bend of two long curved channels, a CACC and b VACC.
Re= 2017 and De= 1008, AR = 2.3, CR = 5.6, Lav / Dh = 39.1, Lentry / Dh = 16.1, and Lexit / Dh = 16.1.
tex and the base vortex, within the boundaries of attachment line
A4 and separation line S7. The onset of the OCW vortex is also
delayed beyond 79 deg due to the existence of split base vortex,
unlike CACC where the OCW vortex appears just after the 45 deg
angular position. Along the outer curvature wall, since the flow
direction of OCW vortex is similar to that of split base vortex,
there is no additional attachment or separation line on OCW
other than S5 due to this Dean vortex. This is another noticeable
difference with CACC compare Fig. 9b with Fig. 6b.
The evolution of the complex flow structures in a VACC and its
difference with CACC are summarized in Table 2. Before closing
this discussion, it is important to re-emphasize that some features
of the complex flow structures is not yet well understood. Further
studies are much needed to shed light on those aspects.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(c)
(b)
(d)
(e)
Fig. 8 Formation of secondary flow at various angular positions along a 90 deg bend of a long diverging-converging VACC.
Re= 2017 and De= 1008 based on uniform flow velocity Vi and hydraulic diameter Dh at channel inlet, and properties at average
temperature Tav = Ti + Tw / 2. Lav / Dh = 39.1, Lentry / Dh = 16.1, and Lexit / Dh = 16.1. AR = 2.3, CR = 5.6. Views: a angular position 11 deg,
b 22.5 deg, c 56 deg, d 67.5 deg, and e 90 deg outlet from the curved portion of the channel or beginning of the straight
constant c/s area exit channel.
11b, the split in base vortex occurs very close to the centerline
heightwise. With increasing axial flow velocity, the location of
the split shifts away from the centerline and the split base vortex
covers a wider area of the channel. At De= 366 in Fig. 11c, an
early sign of formation of another pair of Dean vortex at the OCW
is also observed. The OCW Dean vortex grows in size with increasing De, as seen in Fig. 8e, De= 1008.
10
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(b)
(c)
Fig. 9 Oil flow lines on the curved section 90 deg bend of a long VACC. Re= 2017 and De= 1008, Lav / Dh = 39.1, Lentry / Dh = 16.1,
Lexit / Dh = 16.1. AR = 2.3, and CR = 5.6. a ICW or concave wall, b OCW or convex wall, and c TW. All views are looking into the
wall from front.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 10 Schematics of various 3D flow patterns in the curved section 90 deg bend of a
long VACC. Re= 2017 and De= 1008, Lav / Dh = 39.1, Lentry / Dh = 16.1, Lexit / Dh = 16.1. AR = 2.3, and
CR = 5.6. a Axial flow separation near OCW and the resulting base vortex structure at the
early part of the curved channel, b axial and secondary flow separation from ICW, and
formation of closed separation bubble, and c base vortex, split base vortex, and multiple
Dean vortices at the converging section of the channel.
axial flow separation from ICW in the diverging section, or secondary flow separation from ICW in the converging section. As a
result, Fig. 12 does not show an ICW Dean vortex and separation
bubble which were observed in Fig. 8.
The effect of flow velocity De on the formation of secondary
flow has also been investigated. The generic conclusions are very
similar to that of a long channel. At the lowest velocity condition
091102-12 / Vol. 131, SEPTEMBER 2009
of De= 92, only the base vortex appeared in the channel. The
onset of additional secondary vortices occurs at Decrit = 129, same
as that of the long channel. However, here the onset occurs in the
corner region.
To summarize, the secondary flow in the short, variable c/s
area, curved channel is predominantly characterized by the base
vortex and the split base vortex. The Dean vortices are smaller in
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(b)
(c)
Fig. 11 Effect of axial flow velocity Re and De on the growth of secondary
flow at the curvature outlet at the 90 deg bend of a long VACC. Lav / Dh
= 39.1, Lentry / Dh = 16.1, Lexit / Dh = 16.1. AR = 2.3, and CR = 5.6. a Re= 183, De
= 92; b Re= 366, De= 183; and c Re= 732, De= 366. Re and De are based on
uniform inlet flow velocity Vi, hydraulic diameter Dh at channel inlet, and
properties at average temperature Tav = Ti + Tw / 2.
size and often localized, i.e., they do not cover the entire length of
the curved section.
11
Pressure drop between the channel inlet and exit P is represented in a dimensionless form of friction factor f, defined as
f=
P Dh
1 2 4L
V
2 i
The effects of flow velocity Vi and channel geometry curvature and variable c/s area effect on f are investigated. For the
long channel P includes pressure drop over the entire length
Lentry, Lcurv, and Lexit, while for the short channel it is over Lcurv
only. For both the long and short channels, three different geometric configurations are investigated: the CASC, the CACC, and the
VACC. Traditionally for developing flow, the product of friction
factor and Reynolds number f Re is plotted against the axial
location x, expressed in dimensionless form as x+ = x / Dh Re.
Journal of Fluids Engineering
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
over the curved length Lcurv, the average flow area for VACC is
greater than that of the CACC. Therefore, the P for the former is
expected to be lower than its constant area counterpart, as observed in Fig. 13.
Comparison of the short length VACC with CASC in Fig. 13b
is quite striking. P or f Re for VACC is even less than that for
CASC. This is quite counter intuitive considering the presence of
secondary flow in the curved geometry. At the highest flow rate
condition of Re= 2017 corresponding L+ 0.004, f Re for
VACC is 25% lower than that of CASC. Evidently, the results
suggest that the area expansion effect overwhelms the secondary
flow effect. An order of magnitude, scaling analysis of P can
provide a reasonable explanation of this phenomenon. The inplane velocity vector plots in Fig. 12 show that even though the
secondary flow takes place over the entire channel c/s area, the
velocity boundary layer of the axial bulk flow is very thin. As a
result, the axial flow is uniform over most of the channel c/s area.
Short CACC and CASC also show a similar pattern. Therefore for
a short channel, P should scale with the pressure drop in the
core uniform flow area, or the difference in kinetic energy between the channel inlet and exit, i.e., P Vmax2 Vi2, where
Vmax is the maximum velocity at the channel exit. For the three
geometric cases shown in Fig. 13b, at Re= 2017, Vmax / Vi values
are 1.43, 1.39, and 1.30, respectively, for CACC, CASC, and
VACC. Figure 13b shows that at any Re, f Re follows the
decreasing trend in the same order. Based on the above Vmax valJournal of Fluids Engineering
12
Conclusion
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Acknowledgment
We gratefully acknowledge the support from Rockwell Collins
and the technical interaction with Mr. Matthew Wells and Mr.
Perry Baker. We also thank Mr. Travis Sidebotham and Mr. Benjamine Liu for their help in preparing the sketches and the data
analysis.
Nomenclature
Acs cross-sectional area of the flow channel at any
location variable in the curved section
AR channel aspect ratio= H / b
b channel width; variable, maximum at the 45
deg angular position from the start of the curvature, and for all calculations b value at the
channel inlet is used
CR channel curvature ratio= R / b
De Dean number= Re Dh / R
Decrit critical Dean number for the onset of Dean
vortices= Vi,critDh / Dh / R
Dh channel hydraulic diameter at the flow inlet
= 2bH / b + H
f friction factor= 2P / V2i Dh / 4L
H channel height
ho overall channel heat transfer coefficient= total
heat dissipated/total channel wall surface
areaTw Ti
Hmin magnitude of minimum helicity at a c/s area.
k f thermal conductivity of air temperature
dependent
Lav total channel length along the middle of the
cross section= Lentry + Lcurv + Lexit
L+ dimensionless channel length= Lav / Dh Re
Lcurv average length of the curved section, along the
middle of the cross section
091102-16 / Vol. 131, SEPTEMBER 2009
References
1 Garimella, S. V., 1996, Enhanced Air Cooling for Electronic Equipment, Air
Cooling Technology for Electronic Equipment, S. J. Kim and S. W. Lee, eds.,
CRC, Boca Raton, FL, Chap. 5, pp. 173202.
2 Carter, D. P., Crocker, M. T., Broili, B. M., Byquist, T. A., and Llapitan, D. J.,
2003, Electronics Assemblies With High Capacity Curved Fin Heat Sinks,
U.S. Patent No. 6,671,172.
3 McCormack, P. D., Welker, H., and Kelleher, M., 1970, Taylor Goertler Vortices and Their Effect on Heat Transfer, ASME J. Heat Transfer, 92, pp.
101112.
4 Dean, W. R., 1928, Fluid Motion in a Curved Channel, Proc. R. Soc. London, Ser. A, 121, pp. 402420.
5 Cheng, K. C., Nakayama, J., and Akiyama, M., 1977, Effect of Finite and
Infinite Aspect Ratios on Flow Patterns in Curved Rectangular Channels,
Proceedings of the International Symposium of Flow Visualization, Tokyo,
Japan, pp. 181186.
6 Sugiyama, S., Hayashi, T., and Yamazaki, K., 1983, Flow Characteristics in
the Curved Rectangular Channels Visualization of Secondary Flow, Bull.
JSME, 26216, pp. 961969.
7 Ligrani, P. M., and Niver, R. D., 1988, Flow Visualization of Dean Vortices in
a Curved Channel With 40 to 1 Aspect Ratio, Phys. Fluids, 3112, pp.
36053617.
8 Cheng, K. C., and Akiyama, M., 1970, Laminar Forced Convection Heat
Transfer in Curved Rectangular Channels, Int. J. Heat Mass Transfer, 13, pp.
471490.
9 Ghia, K. N., and Sohkey, J. S., 1977, Laminar Incompressible Viscous Flow
in Curved Ducts of Rectangular Cross-Section, ASME J. Fluids Eng., 99, pp.
640648.
10 Humphrey, J. A. C., Taylor, A. M. K., and Whitelaw, J. H., 1977, Laminar
Flow in a Square Duct of Strong Curvature, J. Fluid Mech., 833, pp. 509
527.
11 Hwang, G. J., and Chao, C.-H., 1991, Forced Laminar Convection in a
Curved Isothermal Square Duct, ASME J. Heat Transfer, 113, pp. 4855.
12 Chandratilleke, T. T., and Nursubyakto, 2003, Numerical Prediction of Secondary Flow and Convective Heat Transfer in Externally Heated Curved Rectangular Ducts, Int. J. Therm. Sci., 42, pp. 187198.
13 Beavers, G. S., Sparrow, E. M., and Magnuson, R. A., 1970, Experiments on
Hydrodynamically Developing Flow in Rectangular Ducts of Arbitrary Aspect
Ratio, Int. J. Heat Mass Transfer, 13, pp. 689702.
14 2006, FLUENT 6.3 Users Guide, Lebanon, NH.
15 Curr, R. M., Sharma, D., and Tatchell, D. G., 1972, Numerical Prediction of
Some Three-Dimensional Boundary Layers in Ducts, Comput. Methods Appl.
Mech. Eng., 1, pp. 143158.
16 Ebadin, M. A., and Dong, Z. F., 1998, Forced Convection, Internal Flow in
Ducts, Handbook of Heat Transfer, 3rd ed., W. M. Roshenow, J. P. Hartnett,
and Y. I. Cho, eds., McGraw-Hill, New York, Chap. 5.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
D. C. Langer
B. A. Fleck
e-mail: brian.fleck@ualberta.ca
D. J. Wilson
Department of Mechanical Engineering,
University of Alberta,
Edmonton, AB T6G 2G8, Canada
Introduction
This study examines the flow field, which occurs when a wall
jet impinges normally onto a forward facing step and is then deflected into a cross-flow. This geometry has not been examined in
the literature and later may be found to have applications for other
mixing or dispersion studies. The preliminary motivation for
studying this geometry was given by Wilson 1, who examined
the dispersion of gases released from pipeline ruptures. Underground pipelines can carry high concentrations of dangerous fluids
such as sour gas or natural gas. Occasionally these pipes fail,
leading to the release of significant quantities of toxic or flammable gasses. The rupture of a buried pipe produces a crater,
where the gas is released at high velocities and impinges onto the
crater wall. The interactions between the jet and the crater lead to
significant spreading and momentum losses in the jet, which can
drastically change the rate of dilution of the released gases in the
surrounding air. These ruptures can produce significant ill-effects
to the environment and danger for the populous surrounding the
rupture area. The ability to model the concentration distributions
of the pollutant within the atmosphere allows for safe placement
of toxic gas pipelines and proper emergency response in the event
of a poisonous gas release 1. Through better understanding the
shape and size of the jet entering into the cross-flowing fluid,
better predictions of the jets behavior can be made. The results of
this study could be used with Wilsons model 1 and be incorporated into hazard assessment models to more accurately determine
the dispersion of hazardous pollutants after a pipeline rupture.
Planar laser induced fluorescence PLIF measurements were
used for both flow visualization and to measure the concentration
profiles within the jet 2. The vertically deflected jet above the
step was studied to determine the shape of the jet entering into the
cross-flow. In preliminary experiments it was found that the jet
exiting the step was elliptical in nature, and the aspect ratio was
dependent on the jet velocity Vjet, cross-flow velocity V, step
height H, and the distance from the step to the jet L. This study
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received September 9, 2008; final manuscript received June 29, 2009; published online August 18, 2009. Assoc. Editor:
Juergen Kompenhans.
Experimental Method
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
pipe was mounted flush with the grid on the bottom of the channel, producing a wall jet, which was centered midway across the
water channel. Based on the logarithmic fit to the velocity profile
performed by Hilderman and Wilson 5, the boundary layer
thickness of the cross-flow at the jet exit location was calculated
to be approximately 15 mm 2d. A more detailed explanation of
the cross-flow is provided by Hilderman and Wilson 5. Dyed
water for the jet was stored in a 75 l tank, which was maintained
at 2.1 bars, and was connected to the pipe using flexible tubing.
Water from the jet and channel were of equal temperature, preventing any buoyancy effects. The velocity of the jet was controlled using a needle valve and a rotameter giving average jet
velocities of 0.47 m/s, 0.78 m/s, and 1.10 m/s. This corresponded
to pipe Reynolds numbers based on diameter of 4.1 103, 7.0
103, and 9.6 103. The step was constructed from a 1.2 m long
sheet of black acrylic, which covered the entire width of the water
channel. The height of the step was adjusted by adding strips of
1.27 cm thick acrylic which ensured that the step remained of
uniform height and was set at 2.54 cm, 3.81 cm, and 5.08 cm.
The front surface was perpendicular to the floor of the water channel and parallel to the face of the jet tube outlet. The jet outlet was
located 5d, 9d, and 15d from the step 43.75 mm, 78.75 mm, and
131.25 mm.
Measurements were taken to determine the shape of the jet
leaving the step with illumination in planes parallel to the water
channel floor. The system was composed of a 4 W Coherent
Innova 70 argon-ion laser, steering optics, a Powell lens, and a
12-bit Cooke SensiCam CCD camera. The laser was run in single
line mode, producing a single beam at 488 nm, with a rated power
of 2.1 W. The beam was passed through a focusing lens, decreasing the thickness of the laser sheet at the center of the water
channel to approximately 1 mm. It was then steered by two mirrors into a 30 deg Powell lens 8. The laser sheet was used to
illuminate fluorescein sodium salt C20H10Na2O5, which was premixed in the jet tank at a concentration of 0.20 0.02 mg/ l. The
absorption peak for this dye is approximately 488 nm, and the
emission peak is at 515 nm 2. A number 12 Kodak wratten
gelatin filter was placed in front of the camera to effectively attenuate the diffracted and reflected laser light from the recorded
images. The filter was rated to have zero transmittance at 488 nm
and approximately 25% transmittance at 515 nm the wavelength
of peak fluorescence emission. The camera was oriented effectively perpendicular to the laser sheet and focused using a 75 mm
Cosmicar TV zoom lens. The resolution of the images was found
to be 0.25 mm/pixel. The experimental signal to noise ratio was
091103-2 / Vol. 131, SEPTEMBER 2009
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(b)
Fig. 3 Single experimental image of the concentration profile within the jet.
A linear concentration scale is shown in a and a sawtooth scale used to
emphasize the structures apparent in the jet is shown in b. The cross-flow
moves from left to right, with the jet fluid moving out of the page. The step
is located at 0 mm.
P = 2a2 + b2
S=
a
b
a b2
2
1/2
1
2
CFD Constraints
15 diameters were chosen to be compared directly with the experimental measurements at the same locations. The intermediate
length was chosen to be ten diameters for the CFD results to
provide a central value for curve fitting. The experimental length
of nine diameters differs from the ten diameters used for the CFD
only because of experimental limitations.
3.1.1 Boundary Conditions. The locations of the boundary
conditions are shown in Fig. 6b. The inlet condition for the
cross-flow had three initial normal speeds: 0.041 m/s, 0.061 m/s,
and 0.081 m/s, all with an initial turbulence intensity of 1%. The
inlet was defined as a subsonic flow with a scalar concentration of
0 mg/l. The jet inlet was also a subsonic flow, with an initial
concentration of 0.7 kg/ m3. The jet velocity was set at three different normal speeds: 0.47 m/s, 0.78 m/s, and 1.10 m/s, with a
turbulence intensity of 5%. An average static pressure condition
was used at the outlet, which was set at 0 Pa relative to the flow.
The bottom of the channel, the step, and the side wall were all set
as smooth walls with the no-slip condition. The top of the channel
was set as a wall with free slip to best simulate the free surface at
the top of the channel.
3.1.2 Meshing. Three computational grids were used to determine the effect of grid dependence: the grids were designated as
coarse, medium, and fine. The grid used for all of the simulations
except the grid convergence study is shown in Fig. 6c. Due to
computer limitations, the coarse grid was used for the parametric
study, and only one simulation was run using each the medium
and fine grids. Since the purpose of the CFD study was to determine the feasibility of using computer simulations to predict the
perimeter of the jet, the usage of the coarse grid was justified. It
was further justified by the parametric study that was analyzed,
which required several different simulations, limiting the available
computer time. An analysis of the effects of using the coarse grid
on the accuracy due to grid dependence is discussed in Sec. 3.1.4.
The three grids used had spacings of 2.5 cm, 2 cm, and 1.5 cm
between nodes in the free stream region. Along the floor of the
channel and on the step, inflation was used to increase the number
of elements near to the walls. In all cases the inflated layer had a
SEPTEMBER 2009, Vol. 131 / 091103-3
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(b)
(c)
(d)
Fig. 4 Experimental a and b and CFD c and d concentration profiles. Linear and sawtooth concentration scales are shown so flow structures can be compared. The cross-flow moves from left to right with the jet
fluid moving out of the page. The step is located at 0 mm.
The advection scheme selected in the CFX software was the high
resolution advection scheme. For this, CFX defines a blending
factor, which will be referred to as but is in the CFX manual,
which represents the order of the numerical scheme, with = 0
representing an upwind differencing scheme and = 1 representing a second order accurate scheme. The high resolution advection
scheme set as close to 1 as possible without introducing oscillations into the flow. This resulted in being variable throughout
the flow with a solution, which was greater than first order accurate throughout the domain though not fully second order. Values
of for a line along the jet centerline for the vertical velocity
component is shown as an example in Fig. 7.
The standard k- model with a scalable wall function was used
for the turbulence within the flow. The coefficients for this model
were C = 0.09, C1 = 1.44, C2 = 1.92, k = 1.0, and = 1.3. The
scalable wall function 10 limits the value of y to be greater
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(b)
Fig. 5 Measured concentrations along the maximum depth line a and the maximum width line b. The
horizontal line shows the 70% value where the depth and width are measured. A threshold value of 60% of the
maximum concentration was used to determine the Gaussian fit.
than 11.06 to avoid placing mesh points within the viscous sublayer which produces the inconsistencies outlined by Grotjans
and Menter 11. Although the k- method has been shown to
overpredict turbulent spreading in round jets and can fail profoundly for three-dimensional flows 12, it provided fast convergence and was sufficiently robust for all of the flow conditions.
Other methods, such as the shear stress transport model and two
Reynolds stress models Baseline Wilcox Model and SpezialeSarkar-Gatski Model, were tested but produced oscillatory results
for the conditions used and would not converge.
A stationary solution was found for the simulation in order to
minimize its computational cost. Since the experiments were performed to determine averaged concentration intensities, a steady
state approximation was all that was required. The timescale and
length scale used in the solution used the default settings: automatic and conservative, respectively. The required convergence
level was set to an average residual of 106. The effect of the
residual on the jet perimeter is shown in Table 1. This shows that
the choice of residual causes less than a 1% variation in the perimeter. The passive scalar transport equation was used to determine the movement of the scalar throughout the flow. It was assumed the concentration variable was a volumetric scalar. The
remainder of the tuning parameters present in the software remained at their default settings.
(a)
(b)
(c)
Fig. 6 Physical definitions in the CFD code. Dimensions were
chosen to match the water channel facility used for the experiments a. The type and locations of the boundary conditions
are given in b. The coarse mesh, which was used for the
simulations with the exception of the grid refinement study, is
given viewed along the symmetry plane in c.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Table 1 Effect of the average maximum residual used for convergence on perimeter measurement. To save computational
time, a residual of 106 was used for the simulations.
Residual
Perimeter
mm
1 1005
1 1006
1 1007
1 1008
220.6
222.9
220.9
220.9
where Fs is a safety factor, for which a value of 3.0 is recommended for conservative estimates of the errors 13. E2 is defined
in Eq. 4
Ecoarse
=
2
r p
1 rp
p = ln
f3 f2
/lnr
f2 f1
Experimental Results
Experiments were done to determine the effects of three parameters: L, H, and R. The results of this parametric study are shown
in Fig. 8. The data plotted here show the average of all of the tests
taken at each flow condition typically two runs. The error bars
were determined by taking two times the standard deviation of
one condition where five tests were run. Figure 8a shows the
relation between perimeter and length is nearly linear for all of the
step heights. Figure 8b shows the effect of the step height, which
is also possesses linearity. Figure 8c shows the perimeter compared with the logarithm of the velocity ratio for a constant step
height H = 38.1 mm and a constant length L = 9d. It can be
seen in these figures that the data is dependent on L, H, and R.
These data were used to create the universal scaling law for this
geometry proposed in Sec. 5.
Discussion
5.1 Empirical Perimeter Correlation. An empirical correlation predicting the perimeter of the jet based on R, L, and H was
determined using the data presented above. In order to better fit
the model, the length L was altered to the distance from a point
momentum source to the step L, located five diameters behind
the jet outlet. This produced the empirical correlation given in Eq.
6, and shown in Fig. 9. From the fit, it can be seen that this
model predicts the perimeter within 13% error for all of the data
points. This correlation provides allows for the perimeter of the jet
entering the cross-flow to be predicted for any H, L, and R within
the range tested
P
HL
= 1.84R0.2
d
d2
1/2
5.79
Both Eq. 6 and Fig. 9 normalize the perimeter by the inlet jet
diameter. This scaling is not fully justified by the experiments, as
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
(a)
(b)
(c)
Fig. 8 Measured values of perimeter compared with the different inlet conditions, normalized with the jet diameter d. The effect of the step position L with different velocity ratios
and a constant step height is shown in a. The effect of the step height H at different
velocity ratios is shown in b. The effect of velocity ratio R for a single step height and
length is shown in c.
S = 8.38
H
L
0.3
R0.4
5.3 Comparison of CFD to Experimental Data. A comparison of the concentration profiles along the measurement plane is
shown in Fig. 4. Parts a and c show the concentration profile
on a linear scale, which appears to have similar shapes for both
the actual and modeled jets. Parts b and d used a sawtooth
concentration profile, which emphasized the structures in the flow.
From this, it can be seen that the simulation accurately predicted
SEPTEMBER 2009, Vol. 131 / 091103-7
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 10 Experimental correlation for the jet aspect ratio S defined in Eq. 7 based on the step height H, effective distance
to the step L, and velocity ratio R
the elliptical shape of the vertical jet, though with slightly larger
depth than the experimental case. The location of the counterrotating vortex pair downstream of the jet is also accurately predicted by the simulation. The comparison of the flow visualization
images shows that the CFD code can clearly predict mean characteristics of this flow.
The perimeters of the ellipses formed by the 70% contours for
all of the experiments and simulations are shown in Fig. 11. A
comparison of the aspect ratios is shown in Fig. 12. This figure
shows significant overlap between the CFD and experimental results. It can be seen that the perimeters predicted by the simulation
do not match those determined experimentally within error. As
can be seen in Fig. 11, the jet perimeter predicted through the
numerical model is significantly larger than the measured values.
According to Pope 12, a well-known deficiency of the
k-epsilon model is that it significantly overpredicts the rate of
spreading for the round jet. This increased spreading would lead
to an increase in the jets perimeter. Since the overprediction of
the spreading is equal in all directions, it has little effect on the
Conclusions
Acknowledgment
The work presented here was made possible by NSERC research grants to B.A.F. and D.J.W., and an NSERC graduate
scholarship to D.C.L.
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Nomenclature
a
b
d
Fs
H
L
L
p
P
r
R
S
V
Vjet
z
References
1 Wilson, D. J., 1979, The Release and Dispersion of Gas From Pipeline Ruptures, Alberta Environment, Contract No. 790686.
2 Walker, D. A., 1987, A Fluorescence Technique for Measurement of Concentration in Mixing Liquids, J. Phys. E, 20, pp. 217224.
3 Cusworth, R. A., and Sislian, J. P., 1987, Computation of Turbulent Free
Isothermal Swirling Jets, University of Toronto Institute for Aerospace Studies, Technical Report No. 315.
4 Behnia, M., Parneix, S., and Durbin, P. A., 1998, Prediction of Heat Transfer
6
7
8
9
10
11
12
13
14
15
16
17
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Sudipto Sarkar
Department of Mechanical Engineering,
Indian Institute of Technology Kanpur,
Kanpur-208016, India
Large-eddy simulations (LESs) of flow past a circular cylinder in the vicinity of a flat
plate have been carried out for three different gap-to-diameter G / D ratios of 0.25, 0.5,
and 1.0 (where G signifies the gap between the flat plate and the cylinder, and D signifies
the cylinder diameter) following the experiment of Price et al. (2002, Flow Visualization
Around a Circular Cylinder Near to a Plane Wall, J. Fluids Struct., 16, pp. 175191).
The flow visualization along with turbulent statistics are presented for a Reynolds number
of Re 1440 (based on D and the inlet free-stream velocity U). The three-dimensional
time-dependent, incompressible NavierStokes equations are solved using a symmetrypreserving finite-difference scheme of second-order spatial and temporal accuracy. The
immersed-boundary method is employed to impose the no-slip boundary condition at the
cylinder surface. An attempt is made to understand the physics of flow involving interactions of shear layers shed from the cylinder and the wall boundary layer. Present LES
reveals the shear layer instability and formation of small-scale eddies apart from their
mutual interactions with the boundary layer. It has been observed that G / D ratio has a
large influence on the modification of wake dynamics and evolution of the wall boundary
layer. For a low gap ratio, it is difficult to identify the boundary layer because of its
strong interactions with the shear layers; however, a rapid transition to turbulence of the
boundary layer, which is similar to bypass transition, is observed for a large gap ratio.
DOI: 10.1115/1.3176982
Keywords: wake-boundary layer interactions, shear layer transition, LES, IB method
Introduction
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
absence of spanwise direction in 2D simulation failed to redistribute the energy in the third direction developing unrealistic large
structures 15.
There are number of experiments illustrating flow around a circular cylinder in the vicinity of a wall, but the discussions were
largely focused on vortex shedding process and modification of
aerodynamic forces on the cylinder. Most of the numerical simulations made so far on this flow configuration are in twodimensional framework and elucidate the flow characteristics with
limited success particularly for high Reynolds numbers. Thus, the
objective of the present paper is to characterize the vortex shedding, either qualitatively or quantitatively, illustrating the flow
physics in terms of vortex stretching, breakdown and turbulence
generation. Here, we try to explain the mutual interactions, i.e.,
modification of wake dynamics due to the presence of the boundary layer and the excitation of boundary layer because of the
wake. It should be noted that the eddy motions and their interactions are expected to be resolved to understand the dynamical
process in turbulent flows for the concerned problem. A direct
numerical simulation DNS would have been ideal for the problem, but it is extremely costly and Reynolds-averaged Navier
Stokes calculations will not provide the answer. Hence, LES can
be pursued as an alternative means with the benefit being a considerably lower computational effort. In the present study, a 3-D
LES with a dynamic subgrid-scale model has been used to investigate the flow around a cylinder in proximity to a wall boundary
layer at Re= 1440 and for G / D of 0.25, 0.5 and 1.0. The Reynolds
number considered here is in the range of shear layer transition
3.
Numerical Methods
u j
=0
xj
1 P 1 2
ui
ij
jui =
ui
u
+
+ fi
+
t xj
xi Re
xj
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
v = 0, u / y = w / y = 0
Computational domain
Blasius
profile
v=w=0
/ t + U / x = 0
Flat plate
5D
u=v=w=0
Boundary
layer
5D
X
25 D
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
G/D = 0.25
4
x/D =
10
15
y/D 2
1
0
3
__
(a)
U /U
Fig. 3 Instantaneous spanwise vorticity z: a and b experiment at Re= 1900 1; c and d present LES at Re= 1440.
a and c G / D = 0.25 and b and d G / D = 0.5.
G/D = 0.25
4
x/D = 2
10
15
y/D 2
1
0
0.1
0.2
(b)
0.3
0.4
0.5
TKE
Fig. 2 Grid independent test for G = 0.25D and Re= 1440 con / U and b turbulent
sidering a mean streamline velocity U
Dimension of Box Lx Ly Lz
0.25
30D 8.5D 3D
0.5
1.0
30D 8.5D 3D
30D 9D 3D
Grid nx ny nz
288 192 32 grid
384 160 32 grid
384 192 32 grid
384 192 64 grid
384 192 32
384 224 32
1
2
3
4
x+
y +
z+
Comments
7.96
3.82
4.35
4.23
4.07
4.87
0.68
0.86
0.71
0.74
0.68
0.8
4.11
3.85
4.41
2.25
4.12
4.93
Test case
Test case
Useful
Test case
Useful
Useful
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
t/T = 0.0
y/D
1
0
3
10
10
t/T = 0.2
y/D
1
0
3
Separation
bubble
0
t/T = 0.4
y/D
1
10
t/T = 0.6
y/D
1
0
3
t/T = 0.0
Inner shear
2
1
0
0
10
10
t/T = 0.8
y/D
y/D
0
Negative roll
0
t/T = 0.2
10
A+B
0
x/D
B
Fig. 6 Span averaged vorticity contours for G / D = 0.5. For details refer Fig. 5.
y/D
1
0
3
t/T = 0.4
Stretching of
inner vortex sheet (C)
10
10
10
y/D
1
0
3
2
Separation
t/T = 0.6
y/D
1
0
3
2
y/D
1
0
10
x/D
Fig. 5 Span averaged vorticity contours for G / D = 0.25. A total
of 20 nondimensional contours are considered in between 5
and +5. The negative vorticity is represented by dotted line.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
t/T = 0.0
Wall vortex
G/D = 0.25
Inner vortex
y/D
Outer vortex
1
y/D
0
3
10
t/T = 0.2
y/D
B
1
0
3
(a)
101
10
x/D
G/D = 0.25
10
t/T = 0.4
100
| z|max
y/D
1
Wall vortex
10-1
Inner vortex
0
3
t/T = 0.6
Interaction of vortex C
and boundary layer
10
10
Roll-up of
boundary layer
t/T = 0.8
y/D
1
0
10
x/D
10-2
(b)
y/D
0
Outer vortex
10
x/D
Fig. 7 Span averaged vorticity contours for G / D = 1.0. For details refer Fig. 5.
= 0.25, vortices from the inner-shear layer and the wall boundary
layer travel in a parallel trajectory owing to the strong pairing
between them Fig. 8a and move in the wall-normal direction
immediately behind the cylinder x / D 3.5. The outer shear
layer vortices move almost parallel to the wall. Near x / D = 4, a
strong mutual interaction between vortices is reflected by the
cross over of trajectories of the shed vortices emerging a single
negative vortex sheet that travels downstream and remains almost
parallel to the wall. It is difficult to identify whether this negative
vortex is produced by the upper shear layer or by the wall vortices
as explained before. Figure 8b indicates that the peak values of
all vortices decay in same pace with increasing streamwise coordinate; they decay rapidly in the vicinity of the cylinder and at a
slower rate downstream of x / D = 4. This is attributed to the break5
Wall vortex
Inner vortex
Outer vortex
G/D = 1.0
4
3
y/D
2
1
0
(a)
10
x/D
101
10
G/D = 1.0
Wall vortex
Inner vortex
Outer vortex
| z|
max
10-1
10-2
(b)
x/D
10
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Upstream
bubble
Secondary
bubble
2.74 D
Primary
bubble
G/D=0.25
3
2
y/D
1
0
-3
-2
-1
1
1.8 D
10
11
10
11
10
11
G/D=0.5
3
2
y/D
1
0
-3
-2
-1
1.1 D
G/D=1
3
2
y/D
1
0
-3
-2
-1
x/D
Fig. 10 Mean streamlines for G / D = 0.25, 0.5, and 1.0
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
G/D = 0.25
3
2
0.2
y/D
0.3
0.4
-0.1
0.1
y/D
-2
3
2
0.2
y/D
0.3
0.4
-0.1
y/D
y/D
1
0.6
-2
-2
0.3
2
-0.4
0.8
0.6
1
0.8
0.6
-0.6
-2
0.2
-0.2
0
-2
0.2
y/D
0.2
0.4
0.6
0.8
-0.6
-0.4
-0.8
4
0
-0.2
-0.4
0.4
0.6
-0.8
-0.4
G/D = 1.0
0.1
-0.4
0.1
y/D
-0.6
-1
-0.6
0
-0.2
-0.4
-0.6
G/D = 1.0
-0.8
G/D = 0.5
Cp
1
0
-0.6
-0.6
0
0.2
-0.2
-0.6
1
0.8
G/D = 0.5
(a)
0.6
0.1
-0.3
0
-0.4
1
0.1
G/D = 0.25
Cp
-2
1.2
x/D
(b)
x/D
Fig. 11 The mean a streamwise and wall-normal velocity contours and b coefficient of pressure
distributions for different G / D. Streamwise velocity contours are drawn by flood and wall-normal
velocity contours are drawn by line black color diagrams. The thick white line indicates zero contour level for mean streamwise velocity.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Primary bubble
Secondary bubble
G/D
Start
End
Length
Start
End
Length
Start
End
Length
0.25
0.5
1.0
2.58
-
0.71
-
1.87
-
1.69
1.06
1.34
10.06
4.57
3.04
8.37
3.51
2.7
2.89
1.95
-
3.64
2.52
-
0.75
0.57
-
(a)
3
y/D = 0.25
z/D
1
0
(b)
-2
10
x/D
Fig. 14 Streamwise velocity contour and vector plot of velocity fluctuations u , v for a front view x-y plane along with a
zoomed section, b top view x-z plane considering a gap ratio of G / D = 1.0. For better visualization, two vectors are
skipped in the x-direction and four vectors in the y-direction,
respectively.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
x/D = 7
y/D = 0.18
R uu
y/D = 0.25
R uu R vv
,
1.5
y/D = 0.32
1
0.5
R vv
0
-0.5
0.3
0.6
0.9
1.2
1.5
(a)
z/D
Fig. 15 Spanwise correlation functions of the streamwise and
the wall-normal velocity fluctuations
mum mean velocity in between the plate and the wake. The experimental data of Matsubara and Alfredsson 38 are superimposed, which were obtained from a transitional boundary layer
subjected to moderate free-stream turbulence. The computed profiles of urms exhibit an approximate self-similarity near the wall
and within x / D = 3.5 to 5.0. High values attributing to the wake
flow are observed in the outer region, Fig. 16a. The agreement
with the measured data 38 appears good illustrating that the
mode of transition and formation of streaks in this case are similar
to those of boundary layer perturbed by strong levels of fst. The
location of urms peak shifts toward the wall depicting the end of
transition downstream of x / D = 7.5, Fig. 16b. In the region of
x / D = 5 7, a large value of urms is noticed, which is attributed to
interactions of large-scale eddies of the shear layers with the
boundary layer.
For G / D = 0.5, some important features of shear layer transition
under the influence of boundary layer can be illustrated by instantaneous velocity fluctuations along with streamwise velocity u
091201-10 / Vol. 131, SEPTEMBER 2009
(b)
Fig. 16 Time-averaged nondimensional urms near the wall: a
between x / D = 3.5 5.0 and b after x / D = 7.5
contours in Fig. 17a: a sectional view along the x-y plane. The
wake flow seems to be turbulent downstream of x / D = 3.0 with
appearance of large- and small-scale structures. The coherent vortices formed due to mutual interactions of shear layers suffer from
(a)
x/D = 1.5
x/D = 5
y/D
x/D = 10
x/D = 15
y/D
(b)
z/D
z/D
Fig. 17 Streamwise velocity contour and vector plot of velocity fluctuations u , v for a front view x-y plane and b side
views y-z plane considering a gap ratio of G / D = 0.5. For details refer Fig. 14.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
10 3
10 1
10 2
10
10
10 -1
10 0
Euu
____
-1
10
-2
10
-3
10
-4
Eww
____
2
U D
-2
10 -3
10
10
(-5/3)
2
U D
(-5/3)
10
G/D = 0.5
-4
G/D = 0.5
10 -5
x/D = 3.3
x/D = 3.3
10 -6
x/D = 5.0
x/D = 5.0
10
-5
10
-6
10
x/D = 10.0
-1
(a)
10
10
10
St
10
-7
10
-8
10
x/D = 10.0
-1
10
(b)
10
10
St
Fig. 18 The power spectra at Re= 1440 for G / D = 0.5: a streamwise velocity and b spanwise
velocity
rms
layer. The presence of the velocity deficit is visible far downstream even up to x / D = 10 and is almost parallel to the wall. As
explained earlier, the strong coupling between the inner-shear
layer and the boundary layer in this case causes the deflection of
urms / U
0.2
0.4
0.6
0.8
1.2
y/D 2
1
0
3
__
U /U
0.8
1.2
urms / U
4
8
(a)
0.2
0.4
0.6
y/D 2
1
0
(b)
3
__
U /U
0.8
1.2
urms / U
4
8
0.2
0.4
0.6
y/D 2
(c)
3
__
U /U
rms
solution
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
x/D
7.5
10
15
20
Kim et al. (1968)
Spalart (1988)
0.2
Production
(a)
0.04
0.08
0.12
0.16
0.2
Zone of
high k
3
y/D 2
3
2
G/D = 0.25
1
0
0.1
(a)
0.1
0.2
0.3
0.4
0.5
0.16
0.2
TKE
__
+ + +
<u/v/> dU / dy
0.3
Production
(b)
20
40
60
80
y+
Fig. 20 Time-averaged nondimensional turbulent kinetic energy productions near the wall and after x / D = 7.5
0.04
0.08
0.12
Zone of
high k
G/D = 0.5
4
3
y/D 2
(b)
0.1
0.2
0.3
0.4
0.5
0.16
0.2
TKE
Production
(c)
0.04
0.08
0.12
G/D = 1.0
Zone of
high k
y/D 2
1
0
1
0
0.1
(c)
0.2
0.3
0.4
0.5
TKE
Conclusions
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
At a moderate gap ratio when G , the inner- and the outerlayer shed from the cylinder tends to curl up in an alternate fashion, although lack of symmetry is strongly evident. Here also the
boundary layer is stretched because of the attraction of the innershear layer evolving an unsteady bubble, the periodicity of which
appears to be same of vortex shedding. The mutual interactions of
the wall vortices and that with the shed vortices create a remarkable difference in the wake dynamics: trajectories of shed vortices
cross each other and an inversion on the position of vortices occurs with respect to the unbounded case. The presence of a wall
causes movements of stagnation and separation points changing
the wake size and base pressure and thus an upward lift on the
cylinder is generated.
When the cylinder is far apart from the wall G = 2.4, the
vortex shedding is similar to the Krmn sheet and the Strouhal
frequency becomes close to 0.2, but asymmetry in wake shape is
still evident. The Krmn vortices excite the boundary layer in an
alternate fashion that triggers its transition.
An attempt is made to describe the transition mechanism of the
shear layers and that of the boundary layer with respect to the gap
ratio. For a relatively low gap ratio G = 0.6 1.2, the shear
layer instability sets in with appearance of a spanwise wavelength
that becomes three-dimensional and breakdowns to turbulent flow
downstream. Here, the coherent vortices of the shear layers suffer
from the intrusion of the deflected boundary layer. The boundary
layer losses its identity downstream because of its strong coupling
with the shear layers. Thus, the generation of turbulence and production in this case is mainly due to shear in the wake fluid illustrating a high active outer region along the wake trajectory and a
relatively calm near-wall. On the other hand, when the cylinder is
well outside the boundary layer G = 2.4, the shear layers and
the boundary layer can be separately identified. The breakdown of
both the outer- and inner-shear layer occurs relatively closer to the
cylinder and the wake tends to retain its coherent structures with
appearance of small-scale eddies. These coherent structures inject
momentum to the boundary layer during their convection that in
turn evolves streaks leading to transition to turbulence of the
boundary layer, an example of receptivity to external disturbances.
Thus, the transition mechanism of the boundary layer for a large
gap ratio is somewhat similar to the bypass transition by fst, although here the excitation is caused by the Krmn vortices.
Downstream of x / D = 7.5, the boundary layer approaches a canonical turbulent layer. This is noteworthy in terms of differences
in response of the boundary layer close to a migrating wake as gap
ratio changes.
Acknowledgment
The authors wish to express their gratitude to the reviewers for
their valuable comments and suggestions.
References
1 Price, S. J., Sumner, D., Smith, J. G., Leong, K., and Paidoussis, M. P., 2002,
Flow Visualization Around a Circular Cylinder Near to a Plane Wall, J.
Fluids Struct., 16, pp. 175191.
2 Norberg, C., 1994, An Experimental Investigation of the Flow Around a
Circular Cylinder: Influence of Aspect Ratio, J. Fluid Mech., 258, pp. 287
316.
3 Williamson, C. H. K., 1996, Vortex Dynamics in the Cylinder Wake, Annu.
Rev. Fluid Mech., 28, pp. 477539.
4 Taneda, S., 1965, Experimental Investigation of Vortex Streets, J. Phys. Soc.
Jpn., 20, pp. 17141721.
5 Bearman, P. W., and Zdravkovich, M. M., 1978, Flow Around a Circular
Cylinder Near a Plane Boundary, J. Fluid Mech., 89, pp. 3347.
6 Grass, A. J., Raven, P. W. J., Stuart, R. J., and Bray, J. A., 1984, The Influence of Boundary Layer Velocity Gradients and Bed Proximity on Vortex
Shedding From Free Spanning Pipelines, ASME J. Energy Resour. Technol.,
106, pp. 7078.
7 Taniguchi, S., and Miyakoshi, K., 1990, Fluctuating Fluid Forces Acting on a
Circular Cylinder and Interference With a Plane Wall, Exp. Fluids, 9, pp.
197204.
8 Lei, C., Cheng, L., and Kavanagh, K., 1999, Re-Examination of the Effect of
a Plane Boundary on Force and Vortex Shedding of a Circular Cylinder, J.
Wind. Eng. Ind. Aerodyn., 80, pp. 263286.
9 Buresti, G., and Lanciotti, A., 1979, Vortex Shedding From Smooth and
Roughened Cylinders in Cross-Flow Near a Plane Surface, Aeronaut. Q., 30,
pp. 305321.
10 Angrilli, F., Bergamaschi, S., and Cossalter, V., 1982, Investigation of Wall
Induced Modifications to Vortex Shedding From a Circular Cylinder, ASME
J. Fluids Eng., 104, pp. 518522.
11 Zdravkovich, M. M., 1985, Forces on a Circular Cylinder Near a Plane Wall,
Appl. Ocean Res., 7, pp. 197201.
12 Vada, T., Nestegard, A., and Skomedal, N., 1989, Simulation of Viscous Flow
Around a Circular Cylinder in the Boundary Layer Near a Wall, J. Fluids
Struct., 3, pp. 579594.
13 Liou, T. M., Chen, S. H., and Hwang, P. W., 2002, Large Eddy Simulation of
Turbulent Wake Behind a Square Cylinder With a Nearby Wall, ASME J.
Fluids Eng., 124, pp. 8190.
14 Dipankar, A., and Sengupta, T. K., 2005, Flow Past a Circular Cylinder in the
Vicinity of a Plane Wall, J. Fluids Struct., 20, pp. 403423.
15 Mittal, R., and Balachandar, S., 1995, Effect of Three-Dimensionality on the
Lift and Drag of Nominally Two-Dimensional Cylinders, Phys. Fluids, 7, pp.
18411865.
16 Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., 1991, A Dynamic
Subgrid-Scale Eddy Viscosity Model, Phys. Fluids A, 3, pp. 17601765.
17 Lilly, D. K., 1992, A Proposed Modification of the Germano Subgrid-Scale
Closure Method, Phys. Fluids A, 4, pp. 633635.
18 Chorin, A. J., 1968, Numerical Solution of the NavierStokes Equations,
Math. Comput., 22, pp. 745762.
19 Zhang, S. L., 1997, GPBi-CG: Generalized Product-Type Methods Based on
Bi-CG for Solving Nonsymmetric Linear Systems, SIAM J. Sci. Comput.
USA, 18, pp. 537551.
20 Mittal, R., and Moin, P., 1997, Suitability of Upwind-Biased FiniteDifference Schemes for Large-Eddy Simulation of Turbulent Flows, AIAA J.,
35, pp. 14151417.
21 Morinishi, Y., Lund, T. S., Vasilyev, O. V., and Moin, P., 1998, Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow,
J. Comput. Phys., 143, pp. 90124.
22 Beaudan, P., and Moin, P., 1994, Numerical Experiments on the Flow Past a
Circular Cylinder at Sub-Critical Reynolds Number, Department Mechanical
Engineering, Stanford University, Report No. TF-62.
23 Kravchenko, A. G., and Moin, P., 2000, Numerical Studies of Flow Over a
Circular Cylinder at ReD = 3900 , Phys. Fluids, 12, pp. 403417.
24 Orlanski, I., 1976, Simple Boundary Condition for Unbounded Hyperbolic
Flows, J. Comput. Phys., 21, pp. 251269.
25 Spalart, P. R., 1986, Numerical Study of Sink Flow Boundary Layers, J.
Fluid Mech., 172, pp. 307328.
26 Fadlun, E. A., Verzicco, R., Orlandi, P., and Mohd.-Yusof, J., 2000, Combined Immersed Boundary Finite Difference Methods for Three Dimensional
Complex Flow Simulations, J. Comput. Phys., 161, pp. 3560.
27 Muldoon, F., and Acharya, S., 2005, Mass Conservation in Immersed Boundary Method, Proceedings of the FEDSM, pp. 19, Paper No. FEDSM-77301.
28 Sarkar, S., 2008, Identification of Flow Structures on a LP Turbine Blade Due
to Periodic Passing Wakes, ASME J. Fluids Eng., 130, p. 061103.
29 Sarkar, S., 2007, The Effects of Passing Wakes on a Separating Boundary
Layer Along a Low-Pressure Turbine Blade Through Large-Eddy Simulation,
Proc. Inst. Mech. Eng., Part A, 221, pp. 551564.
30 Sarkar, S., and Voke, P. R., 2006, Large-Eddy Simulation of Unsteady Surface Pressure on a LP Turbine Blade Due to Interactions of Passing Wakes and
Inflexional Boundary Layer, ASME J. Turbomach., 128, pp. 221231.
31 Sarkar, S., and Sarkar, S., 2007, Large-Eddy Simulations of Cylinder Boundary Layer Interactions, Proceedings of the ACFD7, pp. 13491360.
32 Sarkar, S., and Sarkar, S., 2007, Immersed Boundary Method for Simulating
Complex Flows, Proceedings of the ICFD07, University of Reading, UK.
33 Jacobs, R. G., and Durbin, P. A., 2001, Simulations of Bypass Transition, J.
Fluid Mech., 428, pp. 185212.
34 Ovchinnikov, V., Piomelli, U., and Choudhari, M. M., 2006, Numerical
Simulations of Boundary-Layer Transition Induced by a Cylinder Wake, J.
Fluid Mech., 547, pp. 413441.
35 Zhou, M. D., and Squire, L. C., 1985, The Interaction of a Wake With a
Turbulent Boundary-Layer, J. Aeronaut. Sci., 89, pp. 7281.
36 Zovatto, L., and Pedrizzetti, G., 2001, Flow About a Circular Cylinder Between Parallel Walls, J. Fluid Mech., 440, pp. 125.
37 Wu, X., Jacobs, R. G., Hunt, J. C. R., and Durbin, P. A., 1999, Simulation of
Boundary Layer Transition Induced by Periodically Passing Wakes, J. Fluid
Mech., 398, pp. 109153.
38 Matsubara, M., and Alfredsson, H., 2001, Disturbance Growth in Boundary
Layers Subjected to Free-Stream Turbulence, J. Fluid Mech., 430, pp. 149
168.
39 Spalart, P. R., 1988, Direct Simulation of Turbulent Boundary Layer Up to
R = 1410 , J. Fluid Mech., 187, pp. 6198.
40 Kim, H. T., Kline, S. J., and Reynolds, W. C., 1968, An Experimental Study
of Turbulence Production Near a Smooth Wall in a Turbulent Boundary Layer
With Zero Pressure Gradient, Stanford University, CA, Report No. MD-20.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Paul R. Chiarot
Sergey I. Gubarenko
Ridha Ben Mrad
Pierre E. Sullivan
e-mail: sullivan@mie.utoronto.ca
Department of Mechanical and Industrial
Engineering,
University of Toronto,
5 Kings College Road,
Toronto, ON M5S 3G8, Canada
Introduction
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
duced to a critical curve. The model reveals that the angle at the
center of the interface, the apex angle , dramatically affects the
location of the critical surface in the operational space.
Equilibrium means all forces are balanced on the interface,
the apex angle is 0 deg, and the interface does not emit a spray. In
quasi-equilibrium, all forces are balanced on the interface; however, the apex angle is nonzero and the interface emits a spray see
Fig. 1 for examples of an interface in quasi-equilibrium. Fluids of
high conductivity have a small jet formation region located almost
entirely at the apex of the cone 8. This means that an interface in
quasi-equilibrium can accurately be described as a continuous
function after specifying an apex angle i.e., the jet formation
region can be omitted in the model. Highly conductive fluids are
the focus of both this study and the analytical model 17.
An electrified interface has an operating point that is a function
of the applied voltage, electrode separation, and pressure difference residing inside the operational space. The location of the
operating point defines the important characteristics of the interface. By using the equilibrium model, properties of the interface
can be predicted or selected, including the equilibrium status of
the interface, the presence of an electrospray, and the shape of the
meniscus. The operating point can also move within the operational space and transfer between equilibrium subdomains.
To use the equilibrium model, it is necessary to specify the
value of the apex angle and calculate the pressure difference
across the interface. The sine of the apex angle defines the location of the critical surface, and the pressure difference dictates
along with the applied voltage and electrode separation distance
the location of the interface operating point in the operational
domain. The apex angle term is calculated from the experimental
data by fitting a linear function to several points on the order of
1020 to a portion of the interface in the vicinity of the apex. It
is important to note that the apex angle used here and the classic
Taylor angle T 4 are related as + T = 90 deg.
For an interface under pressure, surface tension, and electric
stress, the pressure term is found by solving a stress balance for
the interface, knowing the curvature of the interface and surface
tension coefficient giving the stress from the surface tension, and
the applied voltage and counterelectrode separation giving the
electric stress. Note that the pressure term is the excess internal
pressure at the interface that balances the surface tension and electric stress. The position of the interface as a function of pressure
difference can be expressed using the stress balance, and compared with the position of the interface at experimentally found
points. By minimizing the difference between the measured and
calculated points, the pressure difference across the interface can
be found 24. The pressure difference for interfaces not under
electric stress can be determined by solving the YoungLaplace
equation. Pressure values are gauge pressures.
The equilibrium model accounts for the observed transient phenomena of an electrified interface. The phenomena considered in
this study are intermittent or pulsed cone-jet mode and smooth
and abrupt transitions of the interface in response to a step voltage. The equilibrium model accounts for these phenomena as
movement of the operating point of the interface within the operational domain and the transition between equilibrium and
quasi-equilibrium states. This approach simplifies the treatment of
these transient phenomena while still accurately describing the
observed behavior.
091202-2 / Vol. 131, SEPTEMBER 2009
3.1 Materials and Experimental Setup. Testing was performed using metal tubing with a radius of 150 m. The electrified interface was located at the edge of the metal tubing that was
positioned directly in front of the counterelectrode. For the imaging of the interface, the bulk fluid used in this study was a
100 M solution of Rhodamine B in 70:30 MeOH: H2O. The
surface tension values were determined using Ref. 25.
Rhodamine B is a fluorescent dye and was used to improve the
recording, analysis, and display of the interface. For the experiments that involved measuring the emitted ion current, the bulk
fluid is a 100 M solution of sodium iodide NaI in 90:10
MeOH: H2O. For some experiments, the bulk fluid is modified by
adding 1% AcOH.
The flow rate was controlled using a syringe pump Cole
Parmer, Montreal, QC, Canada, and the interface was observed
and recorded using an inverted fluorescent microscope Leica,
Wetzlar, Germany and charge-coupled device CCD camera
Sony, Tokyo, Japan. Using optical filters and a dichromatic mirror a filter cube, the microscope can illuminate the interface with
light at a wavelength of 515560 nm, and pass light at a wavelength exceeding 590 nm to the CCD camera. These wavelengths
are ideal for the excitation/emission spectrum of Rhodamine B.
Voltage is applied to the metal tubing using a high voltage source
Labsmith HVS448 with an operating range 03000 V. The radius of the interface in all of the captured images is 150 m,
which is the same as the outer radius of the metal tubing. Images
of the interface are captured at a resolution of 640 480 in PNG
format and processed using MATLAB 26.
The use of Rhodamine B allows for the interface to be more
easily visualized and processed. However, it was typically found
that the emitted jet could not be visualized when illuminated by
filtered light. This is likely because of the narrow bandwidth applied and the small volume of the jet. The jet was visualized when
the full spectrum bandwidth of light is applied. All frames are
aligned so that the bottom of the image is aligned with the end of
the metal tubing. This was confirmed visually and simplifies image processing.
Emitted ion current measurements are made using a current
amplifier Keithley, Cleveland, OH, and the output is monitored
using a digital oscilloscope Agilent, Santa Clara, CA. Ions are
emitted from the apex of the Taylor cone and neutralized on the
counterelectrode. This current was measured, converted to a voltage, and then amplified for measurement.
3.2 Intermittent or Pulsed Cone-Jet Mode. The intermittent
or pulsed cone-jet mode was identified by Cloupeau and PrunetFoch 8 who reported that this occurs when the voltage is
slightly lower than for which a single permanent jet is obtained,
and this mode is one focus of this study. Marginean et al. 19
studied transitions between the operating regimes of an electrospray in experimental studies, and Li 27 investigated meniscus
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 3 Operational domain for an electrode separation distance of 7 mm and surface tension coefficient of 27.48 mN/m.
The operating point moves left and right in the domain across
the critical curves between equilibrium and quasi-equilibrium.
librium state. Due to the low applied flow rate, the interface collapses and returns to the equilibrium state shown by the broken
arrow. Once back in equilibrium, the pressure builds again and
the process repeats itself. This cyclic process is observed to occur
at a very regular frequency of approximately 6 Hz.
The question remains as to why the interface in Fig. 2b, after
moving out of the quasi-equilibrium region, collapses to an equilibrium state rather than moving to a second quasi-equilibrium
state as seen for the smooth transition below. This is because of
the supplied flow rate. For similar operating conditions at higher
flow rates, a stable cone-jet can form 17. Experimentally, it was
found that higher flow rates create larger pressure differences in
the negative direction across the interface. A similar finding was
reported in other electrospray studies 28. At low flow rates, the
pressure difference across the interface remains small, and the
operating point of the interface is limited by the amount it can
move in the negative pressure direction to the left in Fig. 3. As a
result, the operating pointafter crossing into nonequilibrium
must move to the right, back to the equilibrium region. Thus, the
contention of Cloupeau and Prunet-Foch 8 that pulsed cone-jet
mode occurs at voltages slightly below those required for a stable
cone-jet should also include the importance of sufficient pressure
related to flow rate as well. Fernandez de la Mora and Loscertales 10 and Chen and Pui 29 also reported a required minimum flow rate to achieve a stable electrospray. In their work,
scaling laws for the minimum flow rate Qmin are expressed as a
function of the bulk fluid parameters, but they do not specifically
address the relationship of flow rate and pulsed mode.
Ion current measurements confirm the regularity of pulsed
cone-jet mode and the impact of flow rate pressure on the pulse
frequency. Figure 4 shows the measured ion current emitted by an
interface in pulse cone-jet mode for two different flow rate conditions. Figure 4a shows a flow rate of 2 l / min and Fig. 4b
shows a flow rate of 0.5 l / min. In both cases, the bulk fluid is a
SEPTEMBER 2009, Vol. 131 / 091202-3
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
100 M solution of NaI in 90:10 MeOH: H2O. The current between pulses is primarily attributed to external electrical noise. In
Fig. 4, the current peaks correspond to a quasi-equilibrium state of
the interface i.e., Fig. 2b, and the zero current condition corresponds to the equilibrium state i.e., Fig. 2a. The pulsed mode
seen in Fig. 4 is the low frequency component of what Juraschek
and Rollgen 28 identified as axial spray mode I. The low
frequency oscillation of the spraying jet corresponds to the jumping i.e., movement of the operating point in the parameter space
between equilibrium and quasi-equilibrium subdomains. This
measurement shows that the magnitude of the emitted current remains constant despite changes in the pulse frequency, and that
the pulse frequency the emitted ions can be controlled via the
pressure difference across the interface i.e., the supplied flow
rate.
3.3 Smooth and Abrupt Transitions in Response to a Step
Voltage. Experimental evidence reveals that it is possible to have
both smooth and abrupt transitions from equilibrium to final
quasi-equilibrium state. The type of transition for fixed electrode
separation depends on how smoothly the pressure difference and
voltage are applied to the interface. Smoothly applying an increasing pressure difference across the interface is done by operating
the syringe pump at a constant flow rate. Smoothly applying the
voltage is more difficult, as most high voltage supplies including
the unit used in the study do not allow gradual transitions over a
wide operating range, instead, allowing only finite steps. Nevertheless, the impact of varying the spray voltage, as well as the
resultant response of the interface, was previously addressed by
Marginean et al. 19.
Examples of transitions are shown in Figs. 5 and 6. Figure 5
represents an abrupt transition where pressure is smoothly applied, and the voltage is abruptly applied a voltage step from 0 V
up to 3000 V. Figure 6 represents a smooth transition where the
voltage is applied before any fluid is present at 3000 V, and then
pressure is smoothly applied. The time delay between each frame
is shown in each figure.
The transition of the interface in Figs. 5 and 6 can be tracked
using the operational domains shown in Figs. 7 and 8, respectively. The solid and dashed arrows represent smooth transitions
Fig. 7 Operational domain for an electrode separation distance of 10 mm and surface tension coefficient of 27.48 mN/m.
The solid arrow represents the smoothly increasing pressure
difference, and the dotted arrow represents a discontinuous or
abrupt increase in voltage. The application of voltage is discontinuous and not smooth. The pressures at a, b, and d are
known exactly. There is no pressure to report for c since the
interface has ruptured.
Fig. 8 Operational domain for an electrode separation distance of 10 mm and surface tension coefficient of 27.48 mN/m.
The solid arrow represents smoothly increasing pressure difference, and the dashed arrow represents smooth relocations
of the operating points. The transition is smooth; therefore
there is no abrupt loss of mass, as seen in Fig. 5.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
tial mass of fluid is ejected. Since the interface loses mass, the
pressure difference across the interface moves in the negative direction. The operating point after transition Fig. 5d is to the left
of the critical curve for = 50 deg in a quasi-equilibrium state,
since the interface emits a spray after the voltage is applied.
Reznik et al. 16 also found a total rupturing of the interface,
as in Fig. 5c they defined this as droplets jumping off, but
only for large contact angles in excess of 0.8. Here, rupturing
occurs at much smaller angles. Reznik et al. 16 dealt with isolated dropletswhile here droplets are infused with mass from
the syringe pump.
Qualitatively, it has been observed that when the pressure difference across the interface is large before an abrupt transition, the
sudden loss of mass from the interface in transition is large. As
seen in Fig. 5, when the applied pressure before transition is near
the equilibrium limit, the mass lost when the voltage is applied
was large. For smaller pressure differences before abrupt transition, the amount of mass lost from the interface was observed to
be comparatively less.
For the smoothly transitioning interface in Fig. 6, the operating
point starts at zero pressure difference but at an applied voltage of
3000 V the step voltage is applied before any fluid is present.
Mass is continually infused with the syringe pump, and the pressure difference smoothly increases in the positive direction, causing the operating point on the domain in Fig. 8 to start at S and
move to the right. After crossing the critical curve for = 0 deg,
the operating point immediately relocates to the quasi-equilibrium
region, where the pressure increases again. Here the operating
point moves to the right, crossing into nonequilibrium, and adjusts
to find equilibrium. The flow rate in this case is sufficiently high;
therefore the operating point can move to the left for decreasing
pressure, as also noted in Ref. 28, and the interface moves into
a second quasi-equilibrium state. This process repeats until the
interface reaches a final steady operating location.
For the operating domain in Fig. 8, only three discrete transitions into three different quasi-equilibrium regions are shown and
correspond to the three time delays of the interface shown in Fig.
6. It is possible to take more intermediate images of the interface
and plot additional operating points. This would show a series of
transitions between quasi-equilibrium states until reaching the
steady configuration shown in Fig. 6d, and would mean that the
transitional operating points and critical curves would get increasingly closer together.
Overall, the interface in Fig. 6 smoothly transitions to a final
quasi-equilibrium state: cone-jet mode. On the operating domain,
the operating point moves in a smooth manner. Most notably,
there is no dramatic rupture of the interface that is found for the
abrupt transition found in Fig. 5, just a smooth transition to an
interface in cone-jet mode.
The smooth transition can be compared with the pulsed cone-jet
configuration considered in Sec. 3.2. In this case, the interface
does not collapse and is able to return to a quasi-equilibrium state
after a critical curve is crossed. The most notable difference in
operating conditions is the supplied flow rate, where the smooth
transition in Fig. 6 is an order of magnitude larger than the conditions in Fig. 2. With adequate flow rate, the pressure difference
can decrease sufficiently, such that the operating point can move
to the left. The operating point will relocate in a smooth manner
between quasi-equilibrium regions, until a steady value for pressure difference is reached.
Summary
Acknowledgment
The authors gratefully acknowledge the financial support of the
Natural Sciences and Engineering Research Council of Canada
NSERC and Engineering Services Inc. ESI.
References
1 Fenn, J. B., Mann, M., Meng, C. K., Wong, S. F., and Whitehouse, C. M.,
1989, Electrospray Ionization for Mass Spectrometry of Large Biomolecules, Science, 2464926, pp. 6471.
2 Satterley, C., Perdigao, L., Saywell, A., Magnano, G., Rienzo, A., Mayor, L.,
Dhanak, V., Beton, P., and, OShea, J., 2007, Electrospray Deposition of
Fullerenes in Ultra-High Vacuum: In Situ Scanning Tunneling Microscopy and
Photoemission Spectroscopy, Nanotechnology, 18, p. 455304.
3 Barrero, A., and Loscertales, I. G., 2007, Micro- and Nanoparticles Via Capillary Flows, Annu. Rev. Fluid Mech., 39, pp. 89106.
4 Taylor, G. I., 1964, Disintegration of Water Drops in Electric Field, Proc. R.
Soc. London, 2801382, pp. 383397.
5 Basaran, O., and Scriven, L. E., 1990, Axisymmetric Shapes and Stability of
Pendant and Sessile Drops in an Electric Field, J. Colloid Interface Sci.,
1401, pp. 1030.
6 Pantano, C., Ganan-Calvo, A. M., and Barrero, A., 1994, Zeroth-Order, Electrohydrostatic Solution for Electrospraying in Cone-Jet Mode, J. Aerosol Sci.,
25, pp. 10651077.
7 Stone, H., Lister, J., and Brenner, M., 1999, Drops With Conical Ends in
Electric and Magnetic Fields, Proc. R. Soc. London, Ser. A, 455, pp. 329
347.
8 Cloupeau, M., and Prunet-Foch, B., 1994, Electrohydrodynamic Spaying
Functioning Modes: A Critical Review, J. Aerosol Sci., 256, pp. 1021
1036.
9 Fernndez de la Mora, J., 1992, The Effect of Charge Emission From Electrified Liquid Cones, J. Fluid Mech., 243, pp. 561574.
10 Fernndez de la Mora, J., and Loscertales, I. G., 1994, The Current Emitted
by Highly Conducting Taylor Cones, J. Fluid Mech., 260, pp. 155184.
11 Ganan-Calvo, A. M., Davila, J., and Barrero, A., 1997, Current and Droplet
Size in the Electrospraying of Liquids. Scaling Laws, J. Aerosol Sci., 282,
pp. 249275.
12 Cherney, L., 1999, Structure of Taylor Cone-Jets: Limit of Low Flow Rates,
J. Fluid Mech., 378, pp. 167196.
13 Fernandez de la Mora, J., 2007, The Fluid Dynamics of Taylor Cones, Annu.
Rev. Fluid Mech., 39, pp. 217243.
14 Collins, R. T., Jones, J. J., Harris, M. T., and Basaran, O., 2008, Electrohydrodynamic Tip Streaming and Emission of Charged Drops From Liquid
Cones, Nat. Phys., 4, pp. 149154.
15 Notz, P. K., and Basaran, O., 1999, Dynamics of Drop Formation in an
Electric Field, J. Colloid Interface Sci., 213, pp. 218237.
16 Reznik, S., Yarin, A., Theron, A., and Zussman, E., 2004, Transient and
Steady Shapes of Droplets Attached to a Surface in a Strong Electric Field, J.
Fluid Mech., 516, pp. 349377.
17 Gubarenko, S., Chiarot, P., Ben Mrad, R., and Sullivan, P., 2008, Plane Model
of Fluid Interface Rupture in an Electric Field, Phys. Fluids, 204, p.
043601.
18 Zhang, X., and Basaran, O., 1996, Dynamics of Drop Formation From a
Capillary in the Presence of an Electric Field, J. Fluid Mech., 326, pp. 239
263.
19 Marginean, I., Nemes, P., and Vertes, A., 2007, A Stable Regime in Electrosprays, Phys. Rev. E, 76, p. 026320.
20 Gomez, A., and Tang, K., 1994, Charge and Fission of Droplets in Electrostatic Sprays, Phys. Fluids, 61, pp. 404414.
21 Venter, A., Sojka, P., and Cooks, R., 2006, Droplet Dynamics and Ionization
Mechanisms in Desorption Electrospray Ionization Mass Spectrometry, Anal.
Chem., 78, pp. 85498555.
22 Choi, H., Park, J., Park, O., Ferreira, P., Georgiadis, J., and Rogers, J., 2008,
Scaling Laws for Jet Pulsations Associated With High-Resolution Electrohydrodynamic Printing, Appl. Phys. Lett., 92, p. 123109.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
23 Paine, M., Alexander, M., Smith, K., Wang, M., and Stark, J., 2007, Controlled Electrospray Pulsation for Deposition of Femtoliter Fluid Droplets onto
Surfaces, J. Aerosol Sci., 38, pp. 315324.
24 Chiarot, P., Gubarenko, S., Ben Mrad, R., and Sullivan, P., 2008, Application
of an Equilibrium Model for an Electrified Fluid InterfaceElectrospray Using a PDMS Microfluidic Device, J. Microelectromech. Syst., 176, pp.
13621375.
25 Lide, D., ed., 2006, CRC Handbook of Chemistry and Physics, 87th ed., Taylor
and Francis, New York.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Stephen P. Decent
School of Mathematics Edgbaston,
University of Birmingham,
Birmingham, B15 2TT, United Kingdom
Applications of the breakup of a liquid jet into droplets are common in a variety of
different industrial and engineering processes. One such process is industrial prilling,
where small spherical pellets and beads are generated from the rupture of a liquid thread.
In such a process, curved liquid jets produced by rotating a perforated cylindrical drum
are utilized to control drop sizes and breakup lengths. In general, smaller droplets are
observed as the rotation rate is increased. The addition of surfactants along the free
surface of the liquid jet as it emerges from the orifice provides a possibility of further
manipulating breakup lengths and droplet sizes. In this paper, we build on the work of
Uddin et al. (2006, The Instability of Shear Thinning and Shear Thickening Liquid Jets:
Linear Theory, ASME J. Fluids Eng., 128, pp. 968975) and investigate the instability
of a rotating liquid jet (having a power law rheology) with a layer of surfactants along its
free surface. Using a long wavelength approximation we reduce the governing equations
into a set of one-dimensional equations. We use an asymptotic theory to find steady
solutions and then carry out a linear instability analysis on these solutions.
DOI: 10.1115/1.3203202
Introduction
Problem Formulation
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
u=0
ut + u u =
1
1
p + 2 u r
where and p are the density and pressure, t is time, and is the
apparent viscosity such that
=m
:
2
1/2
,v,w
=
u
p =
t + s us + s nu n = S,Bs + Dss2
= ,
m
= 1 n,R,
,R
n
a
and
w = n2w2s, ,t +
p = p0s, ,t + np1s, ,t +
and
,Z
,t, = 1 s,X,Z,Ut,a
,X
s
s0
where U is the exit speed of the jet in the rotating frame, is the
aspect ratio of the jet, and u, v, and w are the tangential, radial,
and azimuthal velocity components relative to the centerline of
the jet, respectively.
We exploit the fact that the radius of the orifice is small when
compared to the radius of the cylinder by expanding u, v, w, and
p in Taylor series in n see Refs. 26,27, and R, X, and Z in
asymptotic series in . We assume that the leading order axial
component of the velocity is independent of and that the centerline of the jet is unaffected by small perturbations to leading
order. Thus we have
n Rs,t,
n Rs,t,
n n =
p
,
U2
1
u, v,w,
U
R = R0s,t + R1s, ,t +
X = X0s + X1s,t + Z = Z0s + Z1s,t +
In addition we have the following expansions for the surfactant
concentration and variable surface tension as
= 0s,t + 1s,t + O2
= 0s,t + 1s,t + O2
n = Rs, ,t
where S = X0sZ0ss X0ssZ0s. The second tangential stress condition is automatically satisfied at leading order and at order 2 we
have
3R20R1v1 + R40w2 + v2 = 0
10
where both these quantities are assumed to vary on the free surface of the jet and only along the jet.
The resulting equations contain a group of nondimensional parameters, which include the Reynolds number Re
= / ms0 U2, the Rossby number Rb= U / s0, and the Weber
number We= U2a / .
The continuity equation at leading order gives v1 = u0s / 2 and
at order n we obtain
w 2 + v 2 = 0
for
u1 = u0S cos
12
13
14
Therefore,
v2 3v2 = u0Ss
u0s
S cos
2
15
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
v2 =
w2 =
16
17
1 u0s
S u0Ss cos
4 2
1 u0s
S u0Ss sin
4 2
= + On
2
su0s
Re
18
3 R0s u0ss
Re 0s
+
u2 = u0s
+
2
R0
4
We 2R0
6
u0sS + u0Ss cos
2
u0s
Re
0
,
R0We
p1 =
R1 + R1
0
+ S cos
R0We
R 02
4
1
v2 +
Re
R0We
R0t +
21
0t + u00s +
24
25
26
u0s
=0
2
27
We see that we have now five equations 22 and 2427 for the
five unknowns u0, R0, 0, X, and Z. The dimensionless equation of
state relating the surfactant concentration to the surface tension
of the liquid-gas interface is given by the nonlinear Szyskowsky
equation
= 1 + log1
28
Steady-State Solutions
We search for steady-state solutions by considering all the variables to be functions of s. In this case we have
uu0s =
+
22
This equation is valid only if the leading order terms in the expansion of X and Z are independent of t. If we had retained leading order translational velocity terms v0 and w0 in our expansions
for v and w, and had X0t 0 and Z0t 0, then the right hand side
of Eq. 22 would contain some additional unsteady terms in E
= ZsXt ZtXs. However E 0 has already been found in Ref. 6 to
be a very accurate approximation between the orifice and the
break up point of the jet. In particular Pru et al. 7 considered
an expansion of the form Xs , t = X0s + Xs , t and Zs , t
= Z0s + Zs , t and then solved the linearized version of the resulting equations. Perturbations of the steady centerline were
found to be within a maximum deviation of order 102 and thus
are relatively small when compared to O1 values of X0s and
Z0s effectively making the trajectory steady. This is backed up
by the experimental observations of Wong et al. 8, which show
that the centerline of the jet is steady in the rotating frame,
giving Xt 0, Zt 0 and E 0 for experimental observations.
Journal of Fluids Engineering
u0s
R0 + u0R0s = 0
2
Xs2 + Zs2 = 1
7 u0s
1
2
u0
+
u20 +
Rb
2 Re
R0We
X + 1Zs ZXs 1
+
u0XsZsss XsssZs
Rb2
Re
20
X0 + 1X0s + Z0Z0s
Rb2
And finally the last equation to be found is the convectiondiffusion equation, which at leading order is
3
20s
R0s
u0ss + 2u0s
+
R
Re
WeR0
19
2
X + 1Zs ZXs
u0 +
cos
p1 = u20S
Rb
Rb2
Re
1 0
We R0
u0t + uu0s =
23
X + 1Xs + ZZs
+
u0ss + 4u2 + u2
2
Rb
Re
1 0
We R0
+
s
20s
WeR0
29
X0sZ0ss X0ssZ0s
=
X0 + 1X0s + Z0Z0s 3
R0s
+
u0ss + 2u0s
2
R
Rb
Re
7 u0s
1
2
u0
+
u20 +
Rb
2 Re
R0We
X0 + 1Z0s Z0X0s 1
+
u0X0sZ0sss X0sssZ0s
Rb2
Re
30
u0s
R0 + u0R0s = 0
2
31
u0s
=0
2
32
and
u00s +
From Eqs. 31 and 32 we observe that R20u0 and 20u0 are constant and by using R00 = 1, 00 = where is the initial surfactant concentration at the nozzle and u00 = 1, we have R20u0
= 1 and 20u0 = 2 so that we can substitute R0 and u0 in the previous equations, which now become
SEPTEMBER 2009, Vol. 131 / 091203-3
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
3
2.5
2
1.5
1
0.5
0.5
4
Z
1 + log1 u1/2
u0s
0 u0s
+
We
2Weu
1
u1/2
2 u0
0
0
u2
X + 1Xs + ZZs 3
+
u0ss 0s
2
Rb
Re
u0
u1/2
0
u0
33
34
35
The unknowns in Eqs. 3335 are X0, Z0, and u0. This nonlinear
system of ordinary differential equations could be solved using a
centered finite difference method as in Ref. 7 for the case where
= 1, however, since it can be shown that viscosity has a negligible influence see Refs. 7,9,28 on the steady-state solutions
for X, Z, and u0 except in the case of highly viscous fluids which
we do not consider here, we will solve Eqs. 3335 in the
simpler inviscid limit i.e., as Re . In this inviscid limit the
resulting nonlinear equations are similar to those found in Ref. 9
except here there are some extra terms due to variable surface
tension. We solve this system using a RungeKutta method with
initial conditions R0 = u0 = X0s = 1 and X0 = Z0 = Z0s = 0 for s = 0. The
effect of changing the initial surfactant concentration on the
steady centerline of the jet is shown in Fig. 1. The trajectory is
seen to curve less when surfactants are present.
= 1 + log1 + = 1 + log1
= e E
1
1
1
2R0ss
0
2 1/2
2 3/2
We
R01 + R0s
1 + R0s
20s
WeR0
38
3
Re
k2 + + iku + iku2
k 2 R 1
k2 + iku
2We R2
Ek2
+ iku = 0
2RWe
39
3
Re
k 2 r
3
2 Re
k2 +
k
2
1
k2
E
R 2 k2
=0
2We
R
R
3 k
Re
2 1
2E
k2 R
We R2
RWe
40
1/2
41
Differentiating the above equation to find the wavenumber for
which r is a maximum which we shall henceforth refer to as k
we find
Temporal Instability
u0t + u0u0s =
and finally
37
X0 + 1X0s + Z0Z0s 3
R0s
+
+
u0ss + 2u0s
2
R0
Rb
Re
=
+
Rb
Rb2
2
2
+ Z0s
=1
X0s
R eikset
7 u0s
1
1 + log1
We
2 Re
1
+
u0X0sZ0sss X0sssZ0s
Re
u0
R = R0 +
0
u0u0s =
36
k =
1
E
2R31/4 2R + 3Oh
1/2
42
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
0.12
85
=0.0
80
breakup length / orifice radius
0.1
= 0.5
0.08
r
= 0.5
0.06
0.04
0.02
75
70
= 0.25,
=0.5
65
60
55
50
= 0.0,
=0.0
45
0
0.2
0.4
0.6
k
0.8
40
0.6
1.2
Fig. 2 The growth rate of disturbances of a rotating shear thinning liquid jet subject to arbitrary disturbance wavenumbers k
at a short distance away from the orifice. The effect of increasing the initial surfactant concentration is seen to lower
growth rates and also to lower the most unstable wavenumber
where r attains a minimum. Here the parameters are We
= 10, Rb= 1.0, Re= 24, = 0.2, and = 0.8.
1
2R301/4
43
3Oh + 2R0
As with previous works see Ref. 9 we find that due to variations in the steady-state variables R and u, both the most unstable
wavenumber and its associated growth rate vary along the jet. In
Fig. 2 we plot the growth rate of disturbances r as given by Eq.
41 at a distance s = 0.3 from the orifice for a shear thinning
rotating liquid jet as the initial surfactant concentration is altered. We find that the effect of adding surfactants is to reduce the
growth rate of disturbances to all unstable wavenumbers, as well
as reducing the wavenumber k at which the growth rates attain a
maximum. Furthermore, the effect of increasing the initial surfactant concentration is to reduce the range of instability. In Fig. 3 we
plot the growth rate of disturbances to arbitrary k for a surfactant
laden shear thickening and shear thinning rotating liquid jet. For
these parameter values we find the shear thinning fluid to have
0.08
= 1.3
= 0.7
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
0.2
0.4
0.6
0.8
0.8
1
1.2
flow index number,
1.4
is shown in
most unstable wavenumber at the breakup point k = kbr
Fig. 7 for a rotating liquid jet at moderately high rotation rates
when the initial surfactant concentration on the free surface is
and
of the value of shear thinning jets have larger values of kbr
that the difference between them decreases as the initial surfactant
concentration is increased. When the jet eventually ruptures it is to
be expected in the context of linear theory that the volume of
SEPTEMBER 2009, Vol. 131 / 091203-5
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
0.8
=
1.0
0.6
=
1.5
=
0.5
0.5
X
1.9
Droplet radius / radius of orifice
0.7
0.4
0.3
0.2
0.1
0
0.2
0.4
0.6
0.8
1.8
1.7
1.6
1.5
1.4
= 0.0
= 0.5
1.3
0.8
0.7
Rb = 1.0
Rb = 4
0.6
0.5
1.5
2
2.5
3
Rossby number, Rb
3.5
duced as the rotation rate is decreased in agreement with Newtonian results obtained by Pru et al. 7. Moreover, it can be
seen that the there exists a nonmonotonic behavior with regard to
the size of droplets as Rb is changed; in particular, we have that
surfactant free jets produce larger droplets for large rotation rates
as compared with jets with surfactants, however, for smaller rotation rates i.e., larger values of Rb this result is vice versa. A
similar qualitative picture is obtained in Fig. 9 for a shear thinning
liquid jet. Finally, in Fig. 10 we show the effect of changing the
Weber number on droplet sizes. Droplets become smaller as the
Weber number is increased i.e., the effects of surface tension are
reduced.
0.5
X
Rb = 0.5
0.4
0.3
0.2
0.1
0
0.2
0.4
Z
0.6
0.8
Conclusions
0.8
= 0.6
= 1.4
1.9
Droplet radius / radius of orifice
0.75
k br
0.7
0.65
0.6
1.8
1.7
1.6
1.5
1.4
= 0.0
= 0.5
1.3
0.55
0.1
0.2
0.3
0.4
0.5
0.6
0.5
1.5
2
2.5
3
Rossby number, Rb
3.5
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
1.9
1.88
1.86
1.84
1.82
1.8
1.78
1.76
1.74
1.72
10
15
20
25
Weber number, We
30
35
Acknowledgment
J.U. would like to thank EPSRC-GB for their financial support.
References
1 Andersen, K. G. and Yttri, G., 1997, Porso vert, Universitets for laget Oslo.
2 Ornek, D., Gurkan, T., and Oztin, C., 2000, Prilling of Aluminium Sulphate,
J. Chem. Technol. Biotechnol., 758, pp. 689694.
3 Wallwork, I. M., Decent, S. P., King, A. C., and Schulkes, R. M. S., 2002,
The Trajectory and Stability of a Spiralling Liquid Jet. Part 1, Inviscid
Theory, J. Fluid Mech., 459, pp. 4365.
4 Wallwork, I. M., 2002, The Trajectory and Stability of a Spiralling Liquid
Jet, Ph.D. thesis, University of Birmingham, Birmingham.
5 Decent, S. P., King, A. C., and Walwork, I. M., 2002, Free Jets Spun From a
Prilling Tower, J. Eng. Math., 42, pp. 265282.
6 Pru, E. I., Decent, S. P., King, A. C., Simmons, M. J. H., and Wong, D. C.
Y., 2006, Nonlinear Travelling Waves on a Spiralling Liquid Jet, Wave Motion, 43, pp. 599618.
7 Pru, E. I., Decent, S. P., Simmons, M. J. H., Wong, D. C. Y., and King, A. C.,
2007, Nonlinear Viscous Liquid Jets From a Rotating Orifice, J. Eng. Math.,
57, pp. 159179.
8 Wong, D. C. Y., Simmons, M. J. H., Decent, S. P., Pru, E. I., and King, A. C.,
2004, Break-Up Dynamics and Drop Sizes Distributions Created From Spiralling Liquid Jets, Int. J. Multiphase Flow, 305, pp. 499520.
9 Uddin, J., Decent, S. P., and Simmons, M. H., 2006, The Instability of Shear
Thinning and Shear Thickening Liquid Jets: Linear Theory, ASME J. Fluids
Eng., 128, pp. 968975.
10 Uddin, J., Decent, S. P., and Simmons, M. H., 2008, Nonlinear Waves Along
Rotating Non-Newtonian Liquid Jets, Int. J. Eng. Sci., 46, pp. 12531265.
11 Uddin, J., Decent, S. P., and Simmons, M. H., 2008, The Effect of Surfactants
on the Instability of a Rotating Liquid Jet, Fluid Dyn. Res., 40, pp. 827851.
12 Hawkins, V. L., Simmons, M. J. H., Brisbane, C. J., Uddin, J., and Decent, S.
P., 2007, Break-Up of Spiralling Non-Newtonian Liquid Jets, Sixth International Conference on Multiphase Flow, Leipzig, Germany.
13 Whitaker, S., 1976, Studies of Drop-Weight Method for Surfactant Solutions
III, J. Colloid Interface Sci., 54, pp. 231248.
14 Hansen, S., Peters, G. W. M., and Meijer, H. E. H., 1999, The Effect of
Surfactant on the Stability of a Fluid Filament Embedded in a Viscous Fluid,
J. Fluid Mech., 382, pp. 331349.
15 Timmermans, M.-L. E., and Lister, J. R., 2002, The Effect of Surfactant on
the Stability of a Liquid Thread, J. Fluid Mech., 459, pp. 289306.
16 Kwak, S., and Pozrikidis, C., 2001, Effect of Surfactants on the Instability of
a Liquid Thread or Annular Layer. Part I: Quiescent Fluids, Int. J. Multiphase
Flow, 17, pp. 137.
17 Blyth, M. G., Luo, H., and Pozrikidis, C., 2006, Stability of Axisymmetric
Core-Annular Flow in the Presence of an Insoluble Surfactant, J. Fluid
Mech., 548, pp. 207235.
18 Ambravaneswaran, B., and Basaran, O. A., 1999, Effects of Insoluble Surfactants on the Nonlinear Deformation and Breakup of Stretching Liquid
Bridges, Phys. Fluids, 115, pp. 9971015.
19 Stone, H. A., and Leal, L. G., 1990, The Effect of Surfactant on Drop Formation and Breakup, J. Fluid Mech., 220, pp. 161186.
20 Blyth, M. G., and Pozrikidis, C., 2004b, The Stability of Two-Layer Channel
Flow With Surfactants, J. Fluid Mech., 505, pp. 5986.
21 Blyth, M. G., and Pozrikidis, C., 2004c, The Effect of Surfactant on the
Stability of Film Flow Down an Inclined Plane, J. Fluid Mech., 521, pp.
241250.
22 Blyth, M. G., and Pozrikidis, C., 2004a, Evolution Equations for the Surface
Concentration of an Insoluble Surfactant; Applications to the Stability of an
Elongating Thread and a Stretched Interface, Theor. Comput. Fluid Dyn., 17,
pp. 147164.
23 Xue, Z., Corvalan, C. M., Dravid, V., and Sojka, P. E., 2008, Breakup of
Shear Thinning Liquid Jets With Surfactants, Chem. Eng. Sci., 63, pp. 1842
1849.
24 Ribe, N. M., 2004, Coiling of Viscous Jets, Proc. R. Soc. London, Ser. A,
460, pp. 32233239.
25 Entov, V. M., and Yarin, A. L., 1984, The Dynamics of Thin Liquid Jets in
Air, J. Fluid Mech., 140, pp. 91113.
26 Eggers, J., 1997, Nonlinear Dynamics and Breakup of Free Surface Flows,
Rev. Mod. Phys., 693, pp. 865929.
27 Hohman, M. M., Shin, M., Rutlidge, G., and Brenner, M. P., 2001, Electrospinning and Electrically Forced Jets. II. Applications, Phys. Fluids, 138,
pp. 22212236.
28 Uddin, J., 2007, An Investigation Into Methods to Control Breakup and Droplet Formation in Single and Compound Liquid Jets, Ph.D. thesis, University
of Birmingham, Birmingham, UK.
29 Rayleigh, L., 1878, On the Instability of Jets, Proc. London Math. Soc.,
s1-10, pp. 413.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Kakani Katija
John O. Dabiri
Graduate Aeronautical Laboratories and
Bioengineering,
California Institute of Technology,
Pasadena, CA 91125
Vortex rings are one of the fundamental flow structures in nature. In this paper, the
generation of circulation and vortex rings by a vortex generator with a static converging
conic nozzle exit is studied numerically. Conic nozzles can manipulate circulation and
other flow invariants by accelerating the flow, increasing the Reynolds number, and by
establishing a two-dimensional flow at the exit. The increase in the circulation efflux is
accompanied by an increase in the vortex circulation. A novel normalization method is
suggested to differentiate between two contributions to the circulation generation: a
one-dimensional slug-type flow contribution and an inherently two-dimensional flow contribution. The one-dimensional contribution to the circulation increases with the square
of the centerline exit velocity, while the two-dimensional contribution increases linearly
with the decrease in the exit diameter. The two-dimensional flow contribution to the
circulation production is not limited to the impulsive initiation of the flow only (as in
straight tube vortex generators), but it persists during the entire ejection. The twodimensional contribution can reach as much as 44% of the total circulation (in the case
of an orifice). The present study offers evidences on the importance of the vortex generator geometry, and in particular, the exit configuration on the emerging flow, circulation
generation, and vortex ring formation. It is shown that both total and vortex ring circulations can be controlled to some extent by the shape of the exit nozzle.
DOI: 10.1115/1.3203207
Keywords: vortex ring, circulation, nozzle, laminar flow, CFD
Introduction
The present paper studies the generation of circulation and vortex rings using conic converging static nozzles in an attempt to
control or manipulate the emerging flow field. The kinematics and
dynamics of vortices formed by starting the flow through parallelwalled tubes were studied extensively in experiments over the
past 70 years 18. Theories and numerical simulations were
implemented to reproduce many of the salient features of the vortex formation process, such as initial vortex sheet rollup 9 and
vortex pinch-off 1014.
A common way of generating vortex rings is by pushing a
piston in a tube 1,35. Two main vortex generator configurations
have been considered: either a tube with an exit protruding into
the flow field, referred to as a nozzle, or a tube with an exit flushed
with a vertical plate, referred to as an orifice. Rosenfeld et al. 10
did not observe significant differences in the maximal total circulation between the two configurations 2%, but the formation
number see definition in Gharib et al. 5 was slightly lower for
the orifice 3.83 versus 3.97 due to vorticity cancellation of the
forming vortex with the side walls. Gharib et al. 5 found experimentally larger differences in the formation number 3.6 and 4.2
for the orifice and nozzle cases, respectively. Mohseni et al. 13
generated numerically vortex rings by applying a nonconservative
force rather than conventional piston/tube arrangements. Nevertheless, the properties of the formed vortex rings were similar to
those generated by a piston in a tube, and they were also affected
by the generation characteristics, such as the amplitude and duration of the forcing.
In addition to these vortex generator types, one can envision
employing nozzles with an axially varying cross section to decrease the exit area, and thus, to control the exit flow conditions,
such as the mean velocity and the velocity profile. Such exit conContributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received April 14, 2009; final manuscript received July 8, 2009; published online August 18, 2009. Assoc. Editor: Mark
Stremler.
figurations can be found in numerous biological systems. Intracardiac blood flow and aquatic propulsion are two areas of active
research in which the dynamics of the starting flow through coniclike nozzle geometries is known to govern the overall system
performance 15,16. The shape of the vortex generator, and in
particular, the exit geometry, may have a significant influence on
the vortex ring formation and pinch-off since the rate of generation of the motion invariants circulation, impulse, and kinetic
energy depend on the details of the flow in the vortex generator.
Indeed, Allen and Naitoh 17 and Dabiri and Gharib 18 manipulated the circulation production rate, vortex circulation, and
energy by dynamically closing an orifice or opening/closing a
nozzle, respectively.
The most widely used analytical tool for predicting the formation kinematics of vortex rings employing parallel wall vortex
generators is the slug model approximation for the vorticity flux
generated by the starting flow 1,4. By assuming the vorticity in
the flow to be confined into a very thin layer near the wall, one
can derive a simple expression relating the circulation flux to the
spatially averaged axial velocity Ue of the jet efflux d / dtt
= 1 / 2U2e t, where t is the time, and is the accumulated circulation. The slug model has had reasonable success in matching
theoretical predictions with empirical measurements of vortex formation from parallel-walled tubes, as long as the discharge duration is short, and the Reynolds number is large 4. However, its
prominent role in current modeling efforts is more closely related
to its ease of use rather than accuracy.
The case of the starting flow through conic nozzle geometries is
inherently beyond the purview of the slug model, due to the existence of nonzero transverse velocity components in the ejected
jet, resulting from the two-dimensional nature of the flow at the
exit. These components are also present during the early stages of
vortex formation from parallel-walled tubes 1,19. However, in
the latter case, their contribution represents a nearly constant shift
in the circulation production, and thus, can be resolved by the
addition of an empirical constant. By contrast, in conic nozzle
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 2 Comparison of the circulation evolution of the numerical solid line and experimental symbols results De / Dp
= 0.6
Methods
2.1 Numerical Model. A sketch of the axisymmetric computational domain is given in Fig. 1. The vortex generator has two
sections: a straight tube of length Ls = 40 cm and a constant diameter of D p = 2.5 cm, and a conic nozzle of length Ln = 5.1 cm the
dimensions are chosen to match the experimental setup and exit
diameter De of De / D p = 0.2, 0.4, 0.6, 0.8, or 1. The computational
domain downstream of the nozzle exit has a length of L / D p = 32,
and the outer boundary is at a radial distance of H / D p = 4. Numerical experiments have verified that the downstream and outer
boundaries are placed far enough from the region of interest.
The moving piston is modeled as a uniform velocity inlet. Undocumented simulations performed in the course of the present
study proved that this model produces results that are very similar
to a moving piston, as long as the maximal stroke length is
smaller than the length of the vortex generator. An impulsive velocity program of the piston is used with a constant piston velocity
of U p = 2 cm/ s in all the cases of the present study. Two velocity
programs are used: i a continuous inflow velocity velocity program no. 1 and ii the constant inflow is stopped at a given
dimensionless time velocity program no. 2. On the outer and
downstream boundaries, zero gauge pressure is specified. On all
the other boundaries, wall conditions are given, except on the axis
of symmetry.
The axisymmetric, incompressible, time dependent, and laminar NavierStokes equations have been solved using the finitevolume package of FLUENT 6.23, ANSYS, Inc. Second-order accurate temporal and spatial schemes have been used with pressureimplicit with splitting of operators PISO pressure-velocity
coupling. Mesh-size and time-step independence tests have been
performed, resulting in a mesh of 180,000 nodes with clustering in
the shear layer and near the walls, and a maximal dimensionless
time step of te = Uet / De0.02, where t is the uniform time
step. The difference in the circulation of the selected mesh relative
to a finer mesh with 685,000 nodes is less than 1%, and relative to
a time step, twice as fine, is less than 0.1%.
2.2 Dimensionless Parameters. The choice of nondimensional parameters is not obvious in the present case since two
velocity scales existone is the piston velocity that remains constant in all the cases and the other is the spatial average exit
velocity Ue that depends on the diameter of the nozzle exit. The
corresponding length scales are the piston diameter and the
nozzle-exit diameter, respectively. Gharib et al. 5 suggested for a
straight tube vortex generator a dimensionless formation time defined by
091204-2 / Vol. 131, SEPTEMBER 2009
t =
U pt
Dp
1a
U pD p
1b
U et
De
and
e =
U eD e
where the spatial average velocity is Ue = U pD p / De2. The relationship between the two formation time definitions is te
= t / De / D p3. As De / D p decreases, the inflow duration to reach a
given te is shorter. For example, to reach a formation time of te
= 4 requires a piston stroke of t = L p / D p = 4, 2.1, 0.85, 0.25, and
0.03 for a nozzle with De / D p = 1, 0.8, 0.6, 0.4, and 0.2,
respectively.
The Reynolds number, based on nozzle-exit parameters, is Re
= 2500, 1250, 830, 625, and 500 for the nozzles with De / D p
= 0.2, 0.4, 0.6, 0.8, and 1.0, respectively. Consequently, in all the
simulations, a laminar flow is assumed, although the case with
Re= 2500 might be transitional 2,20.
2.3 Validation With Experiments. A set of digital particle
image velocimetry DPIV experiments with a nozzle of De / D p
= 0.6 have been performed for validating the numerical simulations. The experimental setup is described in Ref. 7. A nearly
impulsive velocity program, as recorded in the experiments with
a mean velocity of U p 5.7 cm/ s, was imposed in the simulations. The circulation is a major quantity of interest in the present
study, and therefore, the measured and calculated evolution of the
circulation are compared in Fig. 2. Excellent agreement is obtained.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 3 The evolution of the axial velocity profile at the nozzle exit velocity program no. 1
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
t =
t
U
De
D
U
e
= U
2 is the running mean velocity
where the velocity scale U
cl
2 1 t 2
rameters. To test these normalizations of the circulation, the evolution of the total circulation is depicted in Fig. 5 for the
continuous velocity program no. 1. Normalization by the piston or
the exit parameters are plotted in Figs. 5a and 5b, respectively.
The total circulation is calculated by integrating the vorticity field
in the region outside of the straight tube vortex generator.
The circulation that is also proportional to the dimensional
circulation can be significantly increased by decreasing the exit
diameter of the nozzle, Fig. 5a. The dependence of e on De / D p
is smaller. In the latter normalization, the smallest De / D p exhibits
the lowest normalized circulation for a given te . The circulation,
as predicted by the slug model, is also plotted in Fig. 5b. In all
the nozzle cases, the circulation is larger than the slug model
prediction. The deviation decreases as De decreases due to the
more uniform exit velocity profile Fig. 3 that conforms with the
slug model approximation.
None of these normalizations could collapse the evolution of
the total circulation of all the nozzle cases into a single line,
although the same velocity program of the piston is employed. In
an attempt to alleviate this, an alternative normalization of the
circulation and formation time is suggested hereof. It will also
allow separating between two contributions to the ejected circulation, a slug-type one-dimensional flow contribution, and an inherently two-dimensional flow contribution.
4.3 The Two-Dimensional Contribution to Circulation.
The circulation flux ejected out of the nozzle neglecting diffusive
fluxes can be decomposed 1,19 into
d
=
dt
De/2
1 2
+
udr = Ucl
2
De/2
v
dr
x
1
t
+
2 U
D
e
De/2
v
drd
x
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
pcl / Ue versus te , Fig. 7b, yields a very similar plot to Fig. 7a.
The persisting two-dimensional contribution of the constricted
conic nozzles can be also noticed in the nonvanishing asymptotic
value of pcl for De / D p 1. One can notice for De / D p 0.6 that
for a brief period in the formation phase te 1, the 2D effects
contribute a negative circulation flux d / dt 1 / 2 in Fig.
7a, and pcl / U2e 0 in Fig. 7b. For larger exit diameter cases,
the negative two-dimensional contribution cancels out a portion of
the circulation added in the flow initiation phase. This occurs
when the vortex ring is developed enough to induce a significant
negative radial velocity at the exit, reversing the sign of v / x.
With the translation of the vortex ring farther downstream, the
The evolutions of total and vortex ring circulations both normalized by the exit parameters are shown in Fig. 9 for velocity
program no. 2 the inflow is stopped at te = 8. The maximal vortex
ring circulation is given in Table 1 for the various De / D p cases.
For De / D p 0.8, the enhanced circulation efflux increases the
normalized vortex ring circulation as well. For the straight tube
case De / D p = 1, it is e,V = 3.8, significantly larger than the vortex circulation of e,V 2.7 found by Rosenfeld et al. 10 for a
uniform exit velocity profile. In the case of an imposed parabolic
exit velocity profile, they have obtained a larger vortex circulation
of e,V = 3.5, close to the value calculated in the present study for
the straight tube. The relatively low Reynolds number flow in the
present straight tube case Re= 500 results in a parabolic-like exit
velocity profile Fig. 3, and hence, the increase in the vortex
circulation. The augmented vortex circulation is accomplished by
an axis-touching vortex ring that resembles a Hills spherical vortex, Fig. 4. For smaller exit diameter cases De / D p 0.6, the
vortex circulation reduces to e,V = 2.7, identical to the values obtained for a uniform exit velocity profile 10,12.
To further study the formation process, the evolution of the
normalized energy ejected out of the nozzle, as well as the normalized energy of the vortex ring, are plotted in Fig. 10. Following Gharib et al. 5, the normalized energy is defined as e
= E / I3, where E is the kinetic energy, and I is the impulse. The
normalized vortex ring energy is very useful in studying the dy-
Fig. 7 The evolution of a the normalized circulation flux versus the formation time t, and
b the centerline exit pressure versus te
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 9 The evolutions of the total lines and vortex ring lines
with symbols circulations
6
Table 1 The dependence of the vortex ring parameters on the
exit diameter
De / D p
e,V
0.2
0.4
0.6
0.8
1.0
2.73
2.69
2.79
3.10
3.80
4.1
3.9
3.8
3.9
4.1
0.27
0.24
0.22
0.22
0.21
So far, only one factor of the conic nozzle has been considered:
the exit diameter, while the length of the nozzle was kept constant.
The shortening of the nozzle is expanded to increase the twodimensional contribution. A limiting case is a zero-length
nozzlei.e., an orifice, see Fig. 11a. Such a case with an orifice
diameter of De = 1.5 cm, equal to the exit diameter of the conic
nozzle De / D p = 0.6, has also been simulated. It should be noticed
that in most previous vortex ring formation studies that refer to an
orifice case, e.g., Refs. 5,10, the geometry is different. It consists
of a flat plate, flushed with the exit of a straight tube with a small
diameter D p rather than 2H, as in the present case. The most
significant difference is that in the latter case, the flow at the exit
is mostly parallel to the wall of the tube, while in the present case,
there is a large radial component of the flow, upstream of the
orifice, that leads to a large axial gradient of the radial velocity
component, and consequently, to a significant two-dimensional
contribution to the circulation production.
Yet another limiting case that has been simulated is that of a
straight tube vortex generator with a constant diameter equal to
the diameter of the orifice De = D p = 1.5 cm, Fig. 11b. The inflow rate identical to previous conic nozzle cases is employed in
these two additional simulations as well, maintaining the same
exit Reynolds number as of the conic nozzle with De / D p = 0.6
Re= 830.
The exit centerline velocity determines the magnitude of the
one-dimensional contribution, Eq. 3. The evolution of the exit
centerline velocity for these three identical exit diameter cases is
shown in Fig. 12. In all the cases, the initiation of the flow te
1 is accompanied by an overshoot of the exit velocity profile
near the walls, leading to a centerline velocity lower than the
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Straight tube
Conic nozzle
Orifice
3.99
2.79
3.83
4.2
3.8
3.7
0.18
0.22
0.24
e,V
Conclusions
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Nomenclature
D
E
e
F
H
I
L
Lp
Ln
Ls
p
Re
t
te
diameter
kinetic energy g cm/ s2
normalized kinetic energy
formation number
distance of the outer boundary from the axis
cm
impulse g cm/ s
length of the computational domain from the
nozzle exit cm
piston stroke cm
length of the conic segment of the tube cm
length of the straight segment of the tube cm
pressure dyn/ cm2
Reynolds number
time s
formation time based on exit parameters Eq.
2
formation time based on piston parameters
Eq. 1
double-star formation time Eq. 4
the axial and radial velocity components cm/s
t
u, v
U
= U
2 , where U
2 = 1 / tt U2 d
velocity scale U
0 cl
cl
cl
U nozzle exit axial velocity cm/s
x, r axial and radial coordinates cm
Greek Symbols
density g / cm3
azimuthal vorticity 1/s
circulation ejected out of the tube cm2 / s
normalized circulation based on piston parameters Eq. 1
e normalized circulation based on exit parameters Eq. 2
References
1 Didden, N., 1979, Formation of Vortex RingsRolling-Up and Production of
Circulation, Z. Angew. Math. Phys., 301, pp. 101116.
2 Glezer, A., 1988, The Formation of Vortex Rings, Phys. Fluids, 3112, pp.
35323542.
3 Glezer, A., and Coles, D., 1990, An Experimental-Study of a Turbulent Vortex Ring, J. Fluid Mech., 211, pp. 243283.
4 Shariff, K., and Leonard, A., 1992, Vortex Rings, Annu. Rev. Fluid Mech.,
24, pp. 235279.
5 Gharib, M., Rambod, E., and Shariff, K., 1998, A Universal Time Scale for
Vortex Ring Formation, J. Fluid Mech., 360, pp. 121140.
6 Krueger, P. S., and Gharib, M., 2003, The Significance of Vortex Ring Formation to the Impulse and Thrust of a Starting Jet, Phys. Fluids, 155, pp.
12711281.
7 Dabiri, J. O., and Gharib, M., 2004, Fluid Entrainment by Isolated Vortex
Rings, J. Fluid Mech., 511, pp. 311331.
8 Krueger, P. S., Dabiri, J. O., and Gharib, M., 2006, The Formation Number of
Vortex Rings Formed in Uniform Background Co-Flow, J. Fluid Mech., 556,
pp. 147166.
9 Pullin, D. I., 1979, Vortex Ring Formation at Tube and Orifice Openings,
Phys. Fluids, 223, pp. 401403.
10 Rosenfeld, M., Rambod, E., and Gharib, M., 1998, Circulation and Formation
Number of Laminar Vortex Rings, J. Fluid Mech., 376, pp. 297318.
11 Mohseni, K., and Gharib, M., 1998, A Model for Universal Time Scale of
Vortex Ring Formation, Phys. Fluids, 1010, pp. 24362438.
12 Zhao, W., Frankel, S. H., and Mongeau, L. G., 2000, Effects of Trailing Jet
Instability on Vortex Ring Formation, Phys. Fluids, 123, pp. 589596.
13 Mohseni, K., Ran, H. Y., and Colonius, T., 2001, Numerical Experiments on
Vortex Ring Formation, J. Fluid Mech., 430, pp. 267282.
14 Linden, P. F., and Turner, J. S., 2001, The Formation of Optimal Vortex
Rings, and the Efficiency of Propulsion Devices, J. Fluid Mech., 427, pp.
6172.
15 Gharib, M., Rambod, E., Kheradvar, A., Sahn, D. J., and Dabiri, J. O., 2006,
Optimal Vortex Formation as an Index of Cardiac Health, Proc. Natl. Acad.
Sci. U.S.A., 10316, pp. 63056308.
16 Dabiri, J. O., Colin, S. P., and Costello, J. H., 2006, Fast-Swimming Hydromedusae Exploit Velar Kinematics to Form an Optimal Vortex Wake, J.
Exp. Biol., 20911, pp. 20252033.
17 Allen, J. J., and Naitoh, T., 2005, Experimental Study of the Production of
Vortex Rings Using a Variable Diameter Orifice, Phys. Fluids, 17, p. 061701.
18 Dabiri, J. O., and Gharib, M., 2005, Starting Flow Through Nozzles With
Temporally Variable Exit Diameter, J. Fluid Mech., 538, pp. 111136.
19 Krueger, P. S., 2005, An Over-Pressure Correction to the Slug Model for
Vortex Ring Circulation, J. Fluid Mech., 545, pp. 427443.
20 Naitoh, T., Fukuda, N., Gotoh, T., Yamada, H., and Nakajima, K., 2002, Experimental Study of Axial Flow in a Vortex Ring, Phys. Fluids, 141, pp.
143149.
21 Shusser, M., Gharib, M., Rosenfeld, M., and Mohseni, K., 2002, On the
Effect of Pipe Boundary Layer Growth on the Formation of a Laminar Vortex
Ring Generated by a Piston/Cylinder Arrangement, Theor. Comput. Fluid
Dyn., 155, pp. 303316.
22 Heeg, R. S., and Riley, N., 1997, Simulations of the Formation of an Axisymmetric Vortex Ring, J. Fluid Mech., 339, pp. 199211.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Bu-Geun Paik
e-mail: ppaik@moeri.re.kr
Kyung-Youl Kim
Jong-Woo Ahn
Maritime and Ocean Engineering Research
Institute,
KORDI,
Jang-dong 171,
Yuseong-gu,
Daejeon 305-343, Korea
The objective of the present study is to investigate propeller wake using particle image
velocimetry (PIV) technique with bubble type of tracers, naturally generated by the
decrease in the static pressure in a cavitation tunnel. The bubble can be grown from the
nuclei melted in the water tunnel and the size of bubbles is changed by varying the tunnel
pressure. A series of experiments are conducted in the conditions of the uniform and high
velocity gradient flows to find out the characteristics of bubble tracers and compared the
measurement results using bubbles with those using solid particles. Bubbles showed good
trace ability in the region of 15 ReS 75; however, some discrepancies showed at high
velocity gradient region of ReS 1000. The fitted vorticity reduction rate would give
reference for the prediction in a real flow when bubble tracers are utilized in PIV measurements of a vortical flow. In addition, the characteristics of bubble slip velocity can
provide information on the vortex core center and the reduction in the Reynolds shear
stress caused by bubbles deformability. DOI: 10.1115/1.3192136
Keywords: PIV (particle image velocimetry), bubble, slip velocity, cavitation tunnel,
vorticity, turbulence
Introduction
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
3
120
1.0
0.8
0.5
0.3
0.2
5
100
1.0
0.8
0.5
0.3
0.2
8
80
1.0
0.8
0.5
0.3
0.2
In the present study, the propeller wake as a high velocity gradient flow was investigated using bubbles generated naturally by
the decrease in tunnel pressure. In particular, velocity fields were
extracted to confirm the PIV technique with bubble tracers from
both the uniform flow and the propeller wake within the condition
of high Reynolds number over 106. In addition, the reason why
some discrepancies between two tracer cases occurred was investigated in the viewpoint of the slip velocity. Reference data were
also provided to correct vorticity values in the results of bubbly
flow.
ship model has five blades, a mean pitch ratio of 0.990, and a
diameter of 250 mm. The wake screen was arranged in front of the
propeller model to consider the bare hull wake. The significant
decrease in tunnel pressure for the sake of the generation of
bubbles may result in a reduction in the thrust, mismatching the
design condition of the propeller. Therefore, bubble generation
has to be considered after finding appropriate revolution condition. The present study satisfied well the thrust identity.
The analysis of measurement uncertainties or precision errors
was carried out for actual particle images whose flow field was
known in advance. The quiescent flow was tested for the uncertainty evaluation of the present PIV system by following the procedure recommended by Raffel et al. 9 to estimate the measurement errors of PIV systems under the same condition as explained
in the experimental apparatus and method. The standard errors
encountered in measuring the displacement vector were calculated, and the results were summarized in Table 2. To minimize
the measurement uncertainties, 150 instantaneous velocity fields
were ensemble-averaged in the postprocessing. As a result, the
uncertainty levels measured by PIV system were 0.038% and
0.063% for the axial and radial velocity vectors, respectively. Instantaneous velocity fields were obtained for each measurement
phase. They were ensemble-averaged to obtain the time-averaged
in-plane velocity, vorticity, turbulence intensity, and Reynolds
shear stress in terms of tracers for PIV.
In addition, the bubble size was measured by using the shadowgraph technique with the images of bubbles as tracers. Although the shadowgraph technique may give rather less accurate
sizing values than phase Doppler or global rainbow methods, we
employed it because of the absence of other bubble sizing equipments. Two light bulbs were utilized to illuminate the bubbles at
the bottom and side wall. The visualization system consisted of a
high-speed camera Photron, FASTCAM APX-RS, a Nikon 50
mm lens f = 1.4, metal lamps Photron, HVC-SL, an image processor, and a PC. Figure 3 shows one image of the shadowgraph
technique at 0.38 atm. The camera frame rate was varied from
2000 frames/s fps to 6000 fps, since the flow speed was changed
in the range 38 m/s. The spatial resolution and the size of the
measurement plane were 512 512 pixels and 15 15 cm2, respectively. The bubble size was determined from the occupied
pixels by a bubble, the measurement area, and the spatial resolution of the camera.
x
y
root-mean-square
rms
Min.
Max.
0.050
0.082
6.317
8.052
6.367
8.4562
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
PIV analysis. The increment in the flow speed raised the number
of the error vectors in the vector map as shown in Fig. 5c. The
empty spaces in the velocity fields mean errors caused by bubbles.
The velocity field of the flow speed of 8 m/s with the tunnel
pressure of 0.3 atm was similar to that of 0.2 atm and 5 m/s. As
the condition of 8 m/s and 0.2 atm largely increased the error
SEPTEMBER 2009, Vol. 131 / 091301-3
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Table 3 Results obtained from the velocity fields in the uniform flows
U0
Pt Bubble size Error ratio U mean V mean Standard error
m/s atm
mm
%
pixel pixel 103 pixel
3
5
8
1.0
0.2
1.0
0.2
1.0
0.2
0.74 0.24
1.03 0.13
1.60 0.11
0.5
1.0
0.8
1.5
2.8
3.5
3.080
3.077
4.210
4.215
5.487
5.490
0.03
0.02
0.03
0.02
0.05
0.03
9.51
9.98
12.72
10.43
14.16
18.25
vectors caused by the light scattered from the bubble, the appropriate bubble images should be considered in advance according
to the experimental conditions.
Table 3 shows the results of the 400 velocity fields in each
uniform flow case. Here, U0 and Pt mean the free stream velocity
and tunnel pressure, respectively. Only the error ratio equal to the
number of error vectors/number of total vectors was obtained
from the original velocity fields, and other parameters such as
mean value and standard error were extracted after the postprocessing of vector maps. The mean displacement in the vertical
direction shows a nearly constant value of about 0.03 pixels,
which were independent of the tunnel flow speed and the size of
bubbles. In other words, the buoyancy effect did not occur in the
uniform flow conditions. At the atmospheric pressure condition,
the error ratio was below 2%, indicating a good uniformity in the
tracer distribution. In the case of 3 m/s, the error ratio and standard error had small values regardless of the variation in the tunnel pressure. With the flow speed of 5 m/s, the error ratio and the
standard error had values larger than those of 3 m/s because of the
increase in the bubble tracer. However, the PIV measurement
could be employed in the 5 m/s condition because the error ratio
of less than 2% means the acquisition of good particle images.
Although the error ratio and the standard error at 8 m/s were also
larger than those of at 3 m/s, a fine control of the time interval
between the two laser sheets, the laser intensity and the aperture
of the CCD camera would be helpful in acquiring better images.
In the case of the 8 m/s and 0.2 atm tunnel pressure, the error ratio
was larger than the 1.0 atm tunnel pressure. The increase in the
error ratio is attributed to the secondary effects by the light
sources scattered from the rather larger bubbles, which illuminate
the surrounding flow region needlessly. Most of all, it is necessary
to treat the bubble tracers carefully in the case of high flow speed
and low tunnel pressure so that large bubbles would not fill up the
whole interrogation window area.
To measure the high velocity gradient flow trace ability of
bubbles in a complicated flow by using the 2D PIV technique, the
propeller wake was selected because it contains the vortical structure with a high velocity gradient. The wake screen in front of the
rotating propeller provided many tiny bubbles through the small
meshes on it, and this made it possible to reduce the interrogation
window for PIV measurements to 50% overlapped 32 32
pixels. In the condition of 5 m/s and 0.25 atm, the mean diameter
of the bubble was 1.03 mm and the size of the measurement plane
was 20 10 cm2. Since the propeller is working in fixed conditions of above pressure and flow speed, the variation in bubble
size was not available in the propeller wake case. The propeller
wake has a quite complicated flow structure that consists of hub
vortices, tip vortices, and trailing vortices. The hub vortices were
excluded in this study because the diffused reflection from them
was too strong to capture the tracer images. The tip vortices are
generated from the pressure difference between the suction and
pressure sides at the blade tip region. The trailing edge of a blade
produced the trailing vortices, which looked like a curved layer
connected to a tip vortex. Figures 6 and 7 showed the particle
images and the instantaneous velocity fields behind the propellers
obtained at two tunnel pressures and 5 m/s free stream. The tip
091301-4 / Vol. 131, SEPTEMBER 2009
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
vortices and the wake sheets generated from the blade trailing
resemble each other in both of the instantaneous velocity fields. At
the atmospheric pressure, the aperture of the CCD camera was
fully widened to capture the images of the tiny solid particles. On
the other hand, the aperture was narrowed tightly so that it would
not damage the CCD cells at the condition of 0.25 atm and cavitation number of 1.4 because the light intensity scattered from the
bubbles was very strong.
Figure 8a shows the contour plots of phase-averaged axial
velocity at = 0 deg when the solid particles are put into the
tunnel water. The axial velocity component has large values in the
slipstream, which means that the region has high axial momentum
inside the trajectories of the tip vortices and relatively smaller
values near the blade tip and the propeller axis. Along the tip
region 0.5 Y / D 0.4 of each blade, there are pairs of negative and positive isocontours with the shape of a circle because of
the tip vortices. The velocity cores of Rankin type may give the
information on vortex size from the characteristics of vortex. This
tendency is very similar to the contour plots of bubble tracers
Fig. 8b. The result of the bubble tracers shows the contours
location and shapes similar to that of the solid particles. Although
the magnitudes of axial velocity in the slipstream region are a
little different between solid and bubble cases, the bubbles seem
to show a trace ability as good as the solid tracers.
Figure 9 shows the axial velocity profile at several radial locations at = 0 deg. There was a velocity defect, which looks like a
hump in the region 0.3 r / R 0.7 around X / D = 0.3 and 0.6, and
Journal of Fluids Engineering
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 8 The comparison of the phase-averaged axial velocity contours: a solid particles and b
bubble tracers
UmagS =
Ububble Uwater
Ufree
Ububble UwaterDbubble
Here, Dbubble and are the mean diameter of the bubbles and the
kinematical viscosity of water, respectively. The slipstream region
in Fig. 10a gives the range of 15 ReS 75 when Dbubble is
supposed to be about 1 mm. This Reynolds number range is
higher than that 0 ReS 9 of the rising bubble case see Ref.
1 in spite of the similar diameter of bubbles. In the range 4
ReS 100, the drag coefficient reported by Clift et al. 11 is as
follows:
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 10 Contour plots of slip velocity between bubble and solid tracers: a for velocity magnitude
and b for vertical velocity
CD = 13.725 ReS0.74
and the rising bubble case gives 1 CD 20. However, it decreased to a range of 0.56 CD 1.85 when a bubble appeared in
a high Reynolds number flow above the transition region, experiencing a small velocity gradient.
In the case of the small velocity gradient flow, the drag force on
the bubbles was too small to require serious correction of a bubbly
flow into a real flow. On the other hand, the high velocity gradient, especially the tip vortex region 0.50 Y / D 0.35 shown
in Fig. 10a, gives different characteristics of bubbles. The ReS
was about 1000 because the slip velocity magnitude increased to
0.2U0 at a tip vortex. The rotational velocity components were
very strong around the core of a tip vortex. The solid tracers
followed the rotational trace around the core better than the
bubble tracers because the skin friction of the solid tracers is
higher than that of bubbles, and it kept the inertial force continuously even at the region of very high velocity gradient. The lift
force on individual bubble see Ref. 12 could be one potential
source of the high slip velocity and it provides bubbles slip in a
shear and rotational region. The pressure gradient force owing to
density effect could be another source because it accelerates
bubbles in the radial direction of vortices. Those made bubbles not
to follow the rotational trajectory exactly and have a rather low
axial velocity component instead around the tip vortex core,
which was different from the real flow behaviors.
In addition, each tip vortex has a similar vertical slip velocity
VS pattern though the wake moves downstream. The area in which
the bubble has a more vertical velocity than water is sandwiched
between those of VS 0 as shown in Fig. 10b. The solid and
Journal of Fluids Engineering
bubble tracers circulated around the core of the tip vortex tube
because of significant rotational motion and centrifugal force. In
particular, the bubbles have a greater vertical velocity magnitude
because they may move more rapidly and more easily than the
solids in the state of cavitating vortical flow. The region of the
highest value of VS means the center of the tip vortex core, which
can be confirmed in the Z-directional vorticity contours, as shown
in Fig. 11.
The Z-directional vorticity values of the trailing vortices in the
bubble tracer case are similar to those in solid tracers, and the
contours in the bubble case show a slightly broader distribution
than those of solids as shown in Fig. 11. Near the tip vortices, the
vorticity values in the bubble case became smaller because it lost
some of its trace ability owing to an increase in slip velocity.
Actually, bubble concentration around tip vortices is different
from that of the other regions. In the tip vortex region the bubble
concentration was about 184/ 1202 pixels and it reduced to
75/ 1202 pixels at the sandwiched region between tip vortices.
The bubbles accumulated near vortex cores and resulted with the
reduction in velocity magnitudes of bubbles, which is closely related to the decrease in vorticity value at the tip vortex region of
the bubble case. However, the locations of the tip and trailing
vortices are very similar in both tracer cases. These findings led us
to learn that a flow visualization using bubble tracers can be carried out in terms of the spatial configuration of a complicated
flow. For example, the variation in the blade pitch of a propeller
changes the spatial spacing between the tip and trailing vortices in
a wake flow. Although the PIV technique using bubble tracers
cannot measure the vorticity value accurately in a high velocity
SEPTEMBER 2009, Vol. 131 / 091301-7
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 11 The comparison of the phase-averaged vorticity contours: a solid particles and b bubble
tracers
gradient region such as a vortex core, it would be helpful in diagnosing the parameters affecting the spatial evolution of the flows.
Figure 12 shows the phase-averaged vorticity profiles at several
Fig. 12 The comparison of the phase-averaged vorticity profiles at several radial locations
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 14 The comparison of the Reynolds shear stress profiles at several radial locations
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Conclusions
Acknowledgment
This research was sponsored by the Ministry of Knowledge
Economy, Korea under the project code PNS132B and also supported by the basic research project at the MOERI/KORDI
PES128B. In addition, the authors would like to thank Professor
Yuichi Murai Hokkaido University for his valuable comments
and kind advices.
References
1 Brenn, G., Braeske, H., and Durst, F., 2002, Investigation of the Unsteady
Two Phase Flow With Small Bubbles in a Model Bubble Column Using PhaseDoppler Anemometry, Chem. Eng. Sci., 57, pp. 51435159.
2 Lane, G. L., Schwarz, M. P., and Evans, G. M., 2005, Numerical Modeling of
Gas-Liquid Flow in Stirred Tanks, Chem. Eng. Sci., 60, pp. 22032214.
3 Guet, S., Ooms, G., Oliemans, R. V. A., and Mudde, R. F., 2004, Bubble Size
Effect on Low Liquid Input Drift-Flux Parameters, Chem. Eng. Sci., 59, pp.
33153329.
4 Lindken, R., and Merzkirch, W., 2002, A Novel PIV Technique for Measurements in Multiphase Flows and Its Application to Two-Phase Bubbly Flows,
Exp. Fluids, 33, pp. 814825.
5 Jansen, P. C. M., 1986, Laboratory Observation of the Kinematics in the
Aerated Region of Breaking Waves, Coastal Eng., 9, pp. 453477.
6 Govender, K., Mocke, G. P., and Alport, M. J., 2002, Video-Imaged Surf
Zone Wave and Roller Structures and Flow Fields, J. Geophys. Res., 107, p.
3072.
7 Ryu, Y., Chang, K. A., and Lim, H. J., 2005, Use of Bubble Image Velocimetry for Measurement of Plunging Wave Impinging on Structure and Associated Greenwater, Meas. Sci. Technol., 16, pp. 19451953.
8 Paik, B. G., Kim, K. Y., Ahn, J. W., and Kim, K. S., 2008, Velocity Field
Measurements Using Bubble Tracers in a Cavitation Tunnel, Proceedings of
the 18th International Offshore and Polar Engineering Conference, Vancouver, BC, Canada, Jul. 611, pp. 198203.
9 Raffel, M., Willert, C., and Kompenhans, J., 1998, Particle Image Velocimetry,
Springer, New York.
10 Murai, Y., Oishi, Y., Takeda, Y., and Yamamoto, F., 2006, Turbulent Shear
Stress Profiles in a Bubbly Channel Flow Assessed by Particle Tracking Velocimetry, Exp. Fluids, 41, pp. 343352.
11 Clift, R., Grace, J. R., and Weber, M. E., 1978, Bubbles, Drops and Particles,
Academic, New York.
12 Auton, T. R., 1984, The Dynamics of Bubbles, Drops and Particles in Motion
in Liquids, Ph.D. thesis, University of Cambridge, Cambridge, UK.
13 Kitagawa, A., Hishida, K., and Kodama, Y., 2005, Flow Structure of
Microbubble-Laden Turbulent Channel Flow Measured by PIV Combined
With the Shadow Image Technique, Exp. Fluids, 38, pp. 466475.
14 Moriguchi, Y., and Kato, H., 2002, Influence of Microbubble Diameter and
Distribution on Frictional Resistance Reduction, J. Mar. Sci. Technol., 7, pp.
7985.
15 Risso, F., and Fabre, J., 1998, Oscillations and Breakup of a Bubble Immersed in a Turbulent Field, J. Fluid Mech., 372, pp. 323355.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Yoshiki Yoshida1
Kakuda Space Center,
Japan Aerospace Exploration Agency,
Koganezawa 1, Kimigaya,
Kakuda, Miyagi 981-1525, Japan
e-mail: yoshida.yoshiki@jaxa.jp
Yoshifumi Sasao
Institute of Fluid Science,
Tohoku University,
Katahira 2-1-1, Aoba,
Sendai, Miyagi 980-8577, Japan
Mitsuo Watanabe
Tomoyuki Hashimoto
Kakuda Space Center,
Japan Aerospace Exploration Agency,
Koganezawa 1, Kimigaya,
Kakuda, Miyagi 981-1525, Japan
Yuka Iga
Toshiaki Ikohagi
Institute of Fluid Science,
Tohoku University,
Katahira 2-1-1, Aoba,
Sendai, Miyagi 980-8577, Japan
Thermodynamic Effect on
Rotating Cavitation in an Inducer
Cavitation in cryogenic fluids has a thermodynamic effect because of the thermal imbalance around the cavity. It improves cavitation performances in turbomachines due to the
delay of cavity growth. The relationship between the thermodynamic effect and cavitation
instabilities, however, is still unknown. To investigate the influence of the thermodynamic
effect on rotating cavitation appeared in the turbopump inducer, we conducted experiments in which liquid nitrogen was set at different temperatures (74 K, 78 K, and 83 K)
with a focus on the cavity length. At higher cavitation numbers, supersynchronous rotating cavitation occurred at the critical cavity length of Lc / h 0.5 with a weak thermodynamic effect in terms of the fluctuation of cavity length. In contrast, synchronous rotating
cavitation occurred at the critical cavity length of Lc / h 0.9 1.0 at lower cavitation
numbers. The critical cavitation number shifted to a lower level due to the suppression of
cavity growth by the thermodynamic effect, which appeared significantly with rising
liquid temperature. The unevenness of cavity length under synchronous rotating cavitation was decreased by the thermodynamic effect. Furthermore, we confirmed that the fluid
force acting on the inducer notably increased under conditions of rotating cavitation, but
that the amplitude of the shaft vibration depended on the degree of the unevenness of the
cavity length through the thermodynamic effect. DOI: 10.1115/1.3192135
Keywords: thermodynamic effect, cavitating inducer, rotating cavitation, cavity length,
fluid force
Introduction
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
was basically the same as that mentioned in the previous experiment 5. When cavitation develops on the suction side of the
blade at the tip, the pressure sensor shows a particular wave form
in which a flat region appears at the bottom. This region extends
from the suction side to the pressure side within the interblade the
cavitation number decreases. The pressure of this region can be
considered to be the vapor pressure due to the cavitation. Thus, we
judged this domain as the cavity region.
Figure 4 is drawn by the rearrangement of the wave forms of all
the sensors. The dark area in these figures shows estimated cavity
region. Using these figures at each cavitation number, the cavity
length as an indication of cavitation can be observed by this indirect visualization. In a previous work 6, results of this indirect
visualization using the pressure sensors agreed with the direct
optical observation results in water. Uncertainty in Lc / h nondimensional cavitation length measured by these pressure sensors
is 0.03.
The shaft displacement sensor is installed at the backend of the
rotor, as shown in Fig. 2. The sensor for shaft vibration should be
installed near the inducer. However, there is no space around the
inducer due to the limitation of the pump structure
Rotating cavitation can be divided into two patterns based on
the aspect of cavity fluctuations. Figure 4a indicates the cavity
behavior for typical supersynchronous rotating cavitation, and
Fig. 4b indicates typical synchronous rotating cavitation. Unevenness of the cavity length is observed in both figures. The
longer/shorter cavity propagates from blade to blade on the sequence of rotation under supersynchronous rotating cavitation in
Fig. 4a. In contrast, the cavity length on each blade is uneven
but not unsteady under synchronous rotating cavitation in Fig.
4b. It is steady with no propagation during the sequence of rotation.
3
Fig. 2 Schematic diagram of the test inducer showing the installed pressure sensors and shaft displacement sensor
3.1 Rotating Cavitation. As we have already mentioned, rotating cavitation is cavitation instability and has two patterns. Figure 5 shows examples of the cavitation state under conditions of
supersynchronous Fig. 5a and synchronous rotating cavitation
Fig. 5b by direct optical observation in water experiments 7.
In Fig. 5a, the unevenness of cavity length medium, short and
long on each blade blades 13 can be clearly observed. However, supersynchronous rotating cavitation is an unsteady phenomenon in which the asymmetric pattern of the cavity length
propagates from blade to blade at a speed 1.11.2 times higher
than that of inducer rotation. In contrast, the unevenness of cavity
length long, long and short on each blade blades 13 can be
also observed in Fig. 5b. Though the cavitation state in the inducer is surely asymmetric and uneven, it does not change
from blade to blade in rotation. From that viewpoint, synchronous
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
RR + 3 R 2 + R t = C p +
2
2
where
=
Fig. 5 a Photographs of supersynchronous rotating cavitation in water from Ref. 7. b Photographs of synchronous
rotating cavitation in water from Ref. 7.
rotating cavitation can be defined as a steady anomalous phenomenon. These direct optical observations in water agree with
the indirect observations in liquid nitrogen shown in Figs. 4a
and 4b.
3.2 Thermodynamic Effect. The introduction of the temperature depression into the simple RayleighPlesset equation for
a spherical cavity yields 3
3
pvT p
RR + R2 + Rt =
2
l
The third term on the left is the thermal term caused by the
thermodynamic effect. Thermodynamic function , which depends only on temperature T, was originally introduced by Brennen 8. The dimension of thermodynamic function is in m / s3/2
T =
vL2
l2C plTl
C
U3
Experimental Results
Nitrogen
Temperature
K
Rotation N
rpm
m / s3/2
= C / U31/2
74
78
83
18,300
18,300
18,300
1.4 10+04
3.3 10+04
8.6 10+04
4.3
10.2
26.6
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 8 a Fluctuations of cavity lengths for each channel in the case of supersynchronous
rotating cavitation left: 74 K, right: 83 K uncertainty in Lc / h = 0.03 and b fluctuations of
cavity lengths for each channel in the case of synchronous rotating cavitation left: 74 K, right:
83 K uncertainty in Lc / h = 0.03
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 9 Influence of temperature on cavity length and the region where the rotating cavitations occurs at 74 K, 78 K, and 83
K uncertainty in / 0 = 0.02 and Lc / h = 0.03
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 11 a Instant vector orbits of fluid force under typical supersynchronous rotating cavitation at 74 K and 83 K uncertainty in F / Fref = 0.03 and b time-averaged vector orbits of fluid
force during synchronous rotating cavitation at 74 K and 83 K uncertainty in F / Fref = 0.03
Discussions
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Conclusions
By experiments with liquid nitrogen set at different temperatures and considerations based on the cavity length as an indication of cavitation, the following points were clarified.
1 Rotating cavitation occurs at the same cavity length independent of temperature. The critical cavity length at the
onset of supersynchronous rotating cavitation was Lc / h
0.5, and that at the onset of synchronous rotating cavitation was Lc / h 0.9 1.0 in the present inducer.
2 The thermodynamic effect suppresses the growth of cavity
length. Growth of the cavity is distinctly suppressed with
increasing liquid temperature.
3 Due to the above, the thermodynamic effect shifts the critical cavitation number of the rotating cavitations lower.
4 The differences in the behavior of the cavity fluctuation and
oscillating frequency under supersynchronous rotating
cavitation were not prominent at any of temperatures.
5 However, unevenness of the cavity length under synchronous rotating cavitation was greatly affected by the thermodynamic effect in the present experiments of flow rate
Q / Qd = 1.06.
6 The favorable thermodynamic effect on cavitation instability becomes greater when the cavity length extends over the
throat at lower cavitation numbers.
7 The fluid force acting on the inducer significantly increases
under supersynchronous/synchronous rotating cavitation.
The thermodynamic effect influences considerably the unevenness of the cavity length under synchronous rotating
cavitation. Therefore, synchronous shaft vibration under
synchronous rotating cavitation is much affected by the
temperature of working fluid.
Acknowledgment
Authors would like to thank Mr. Yusuke Kazami of Tohoku
University for his great help with data processing.
Nomenclature
C
Cp
C pl
F
Fref
h
L
Lc
p
pv
R
Tc
T
t
U
l
v
l
chord of inducer
pressure coefficient
specific heat of liquid
fluid force
reference value of fluid force
blade spacing
latent heat
cavity length
pressure
vapor pressure
radius of bubble
critical temperature
ambient temperature
time
tip velocity of inducer
thermal diffusivity of liquid
vapor density
liquid density
thermodynamic function defined by Eq. 2
nondimensional thermodynamic parameter by
Eq. 4
cavitation number
reference value of cavitation number upper
limit of cavitation number in the present
experiment
transit time =Lc / U
head coefficient
reference value of head coefficient
angular velocity of rotor
angular velocity of rotating cavitation
References
1 Rosenmann, W., 1965, Experimental Investigations of Hydrodynamically Induced Shaft Forces With a Three-Bladed Inducer, Proceedings of the Symposium on Cavitation in Fluid Machinery, Winter Annual Meeting, ASME, New
York, pp. 172195.
2 Kamijo, K., Yoshida, M., and Tsujimoto, Y., 1993, Hydraulic and Mechanical
Performance of LE-7 LOX Pump Inducer, AIAA J., 96, pp. 819826.
3 Franc, J. P., Rebattet, C., and Coulon, A., 2004, An Experimental Investigation of Thermal Effects in a Cavitating Inducer, ASME J. Fluids Eng., 126,
pp. 716723.
4 Cervone, A., Testa, R., and dAgostino, L., 2005, Thermal Effects on Cavitation Instabilities in Helical Inducers, AIAA J., 215, pp. 893899.
5 Yoshida, Y., Sasao, Y., Okita, K., Hasegawa, S., Shimagaki, M., and Ikohagi,
T., 2007, Influence of Thermodynamic Effect on Synchronous Rotating Cavitation, ASME J. Fluids Eng., 1297, pp. 871876.
6 Yoshida, Y., Kikuta, K., Hasegawa, S., Shimagaki, M., and Tokumasu, T.,
2007, Thermodynamic Effect on a Cavitating Inducer in Liquid Nitrogen,
ASME J. Fluids Eng., 1293, pp. 273278.
7 Maekawa, Y., 1996, Experimental Study of Unsteady Cavitation on an Inducer, MS thesis, Osaka University, Osaka in Japanese.
8 Brennen, C. E., 1973, The Dynamic Behavior and Compliance of a Stream of
Cavitating Bubbles, ASME J. Fluids Eng., 95, pp. 533541.
9 Watanabe, S., Hidaka, T., Horiguchi, H., Furukawa, A., and Tsujimoto, Y.,
2007, Analysis of Thermodynamic Effects on Cavitation Instabilities,
ASME J. Fluids Eng., 1299, pp. 11231130.
10 Kato, H., 1984, Thermodynamic Effect on Incipient and Development of
Sheet Cavitation, Proceedings of International Symposium on Cavitation Inception New Orleans, LA, Dec. 914, ASME, New York, Vol. 16, pp. 127
136.
11 Zoladz, T., 2000, Observations on Rotating Cavitation and Cavitation Surge
From the Development of the Fastrac Engine Turbopump, Proceedings of the
36th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, Huntsville, AL, Jul.
1719, Paper No. AIAA200-3403.
12 Kimura, T., Yoshisa, Y., Hashimoto, T., and Shimagaki, M., 2008, Numerical
Simulation for Vortex Structure in a Turbopump Inducer: Close Relationship
With Appearance of Cavitation Instabilities, ASME J. Fluids Eng., 130, p.
051104.
13 Watanabe, S., Furukawa, A., and Yoshida, Y., 2008, Theoretical Analysis of
Thermodynamic Effects on Cavitation in Cryogenic Inducer Using Singularity
Method, Proceedings of the 12th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Honolulu, HI, Feb. 1722,
Paper No. ISROMAC12-2008-20251.
14 Kobayashi, S., 2006, Effects of Shaft Vibration on Occurrence of Asymmetric
Cavitation in Inducer, JSME Int. J., Ser. B, 494, pp. 12201225.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
A. K. Das
e-mail: arup@mech.iitkgp.ernet.in
P. K. Das1
e-mail: pkd@mech.iitkgp.ernet.in
Department of Mechanical Engineering,
IIT Kharagpur,
721302, India
J. R. Thome
LCTM,
EPFL,
Lausanne CH-1015, Switzerland
e-mail: john.thome@epfl.ch
The two fluid model is used to simulate upward gas-liquid bubbly flow through a vertical
conduit. Coalescence and breakup of bubbles have been accounted for by embedding the
population balance technique in the two fluid model. The simulation enables one to track
the axial development of the voidage pattern and the distribution of the bubbles. Thereby
it has been possible to propose a new criterion for the transition from bubbly to slug flow
regime. The transition criteria depend on (i) the breakage and coalescence frequency, (ii)
the bubble volume count below and above the bubble size introduced at the inlet, and (iii)
the bubble count histogram. The prediction based on the present criteria exhibits excellent agreement with the experimental data. It has also been possible to simulate the
transition from bubbly to dispersed bubbly flow at a high liquid flow rate using the same
model. DOI: 10.1115/1.3203205
Keywords: bubbly flow, population balance, transition, coalescence, breakage
Introduction
fluid model 5. There have been many variations of this technique. An important idealization has been made by the researchers
particularly for the cases when one phase is thoroughly dispersed
in the other. This is known as interpenetrating continua or mixed
continua 6. According to the hypothesis of the interpenetrating
continua, both the phases are simultaneously present at each and
every point of the computational domain, while their influence on
the hydrodynamics and the transport processes is given by the
local phase fraction.
In multiphase flow, the phase distribution changes both with
time and space even if the input and the operating conditions are
kept invariant. Though the two fluid model is a considerable success toward the modeling of multiphase flow, it needs augmentation to take care of the variable phase distribution. Several methodologies have been suggested for this. The interfacial area
transport equation was suggested by Ishii 7. In 1998 Wu et al.
8 developed one grouped interfacial area transport equation for
spherical bubbly flow pattern. Further, it is realized that the shape
of the bubbles distorts during motion causing a change in the drag
force and the fluid particle interaction mechanisms. Two-group
spherical/distorted bubble and cap/churn turbulent bubble area
transport equation was developed 9 to circumvent this.
An alternate approach of simulation through population balance
model PBM 10 has also been adopted. PBM essentially considers the secondary phase as the assemblage of discrete class or
subgroup having different characteristics lengths. Additionally,
PBM can meticulously track the birth and death of the members
of these subgroups. Various methodologies, such as method of
discrete classes 10, quadrature method of moments 11, least
square methods 12, etc., have been developed to implement
PBM numerically.
As this model has a wide applicability with significant accuracy, a number of researchers 5,8,13 have adopted the technique
for the simulation of bubbly flow with different modifications.
Chen et al. 14 showed the applicability of the model in case of
bubble columns of different diameters. Similarly, Yeoh and Tu
15 and Cheung et al. 16 solved the population balance equation along with the two fluid model using the method of discrete
classes 10 for pipe flow. An extensive review of the works related to population balance modeling for bubbly flows in the
bubble column can be found in Ref. 17, while Krepper et al. 18
reviewed the work for pipe flow.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Furthermore, for a better prediction of the assemblage of different sized entity multiple-sized group models MUSIGs has also
been tried by different researchers 1820. In these models, different groups based on size of bubbles are assumed to have
different velocities. A separate set of conservation equations are
solved for each group of the bubbles, considering each one of
them as one distinct continuum. However, the increase in computation time and the need for a large number of closure relationships cannot be overlooked in this type of simulations.
Efforts have also been made to directly solve the two phase
interacting momentum balance equation by using sufficiently fine
grids and smaller time steps in the problem domain. Esmaelli and
Tryggvason 21 applied direct numerical simulation DNS for
laminar bubbly flow in a vertical channel. After that, continuous
efforts 2224 have been made to improve the numerical scheme
for a better prediction. Nevertheless, due to the volume of computations involved, only small systems with fewer complexities
can be modeled at present 25. As a consequence, models based
on averaging both in time and space, such as the two fluid
model, are very important for simulating problems of practical
importance.
It is clear from the available literature that the application of a
two fluid model along with PBM or interfacial area transport
equation has produced reasonable success in the prediction of
bubbly flow. However, to date not much effort has been made to
predict the transition of cocurrent bubbly flow through a circular
conduit using computational fluid dynamics CFD. On the other
hand, CFD is rigorously used by researchers 2629 to determine
the transition of homogeneous to heterogeneous regime for bubble
column reactors. This motivated us for the numerical investigations of transition regimes of bubbly flow situations. Bubbly flow
through a vertical tube is bound by two neighboring flow regimes.
At the low and moderate liquid flow rate, bubbly flow transforms
into slug flow with the increase in air flow rate, whereas transition
from bubbly to dispersed bubbly flow is observed at high liquid
flow rates 30. Classically, these transition criteria have been derived from the control volume approach using a number of simplifying assumptions 31,32. From these criteria of transition are
predicted based on the phase superficial velocities without taking
other process parameters such as tube diameter, inlet bubble size
distribution, etc. into cognition. Furthermore, a considerable
variation has been reported in the transition boundaries observed
in different experiments. This has motivated the present work to
apply CFD technique to identify the suitable transition criteria,
which are rigorously based on hydrodynamics. The bubbly flow
has been simulated using the two fluid model along with the population balance technique. Though the basic structure of the model
has been retained, a few important modifications have been made
to simplify the computation and to improve the simulation. The
model simulation has been used to investigate the dynamic development of coalescence and to break up processes along the axial
direction. This has been studied noting the breakage and coalescence frequency, bubble volume account below and above the
bubble size at the inlet, and bubble count histogram. It has been
possible to propose unique transition criteria from bubbly to slug
flow. The transition boundary predicted by the present model
matches very well with the experimental data. Finally, the present
model has been extended to the zone of high liquid flow rate to
predict the transition to dispersed bubbly flow. Experimentally
observed bubbly to dispersed bubbly transition is predicted satisfactorily.
1.1 Model Development. Figure 1 depicts the typical up flow
of gas-liquid mixture with low superficial velocities, where the
lighter phase is dispersed in the form of bubbles of different sizes
in a round vertical conduit. The constitutive equations of each
phase are developed based on the volume average properties of
the phase. The total population of bubbles is discretized into subgroups based on volume to avoid intermediate bubble generation
during coalescence. In the present model each group of bubbles
091303-2 / Vol. 131, SEPTEMBER 2009
1
riiui + iiwi = 0
i i +
t
r r
z
For r momentum,
+ i ig r
iiui + iiui2 + iiuiwi = i
r
t
r
z
+ i
2
2
i iu i i
iu i + i 2 iu i
+
2 iu i +
2
r
z
r
r r
For z momentum,
+ i ig z
iiwi + iiuiwi + iiwi2 = i
z
t
r
z
+ i
2
2
i
iw i + i 2 iw i
2 iw i +
r
z
r r
Here i = l stands for liquid phase, and i = g stands for gas phase.
P is assumed as the average pressure of the mixture, and it can be
evaluated using the equation of state in the following manner:
P = ll + ggRT
Constitutive relationships for the force terms used in the momentum equations are selected from the published literature. As the
generic nature of the radial and axial forces is the same, directions
of the forces are omitted in the respective equations Eqs.
510. For the vector quantities associated with the forces, we
Transactions of the ASME
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
2CFL
gug ulug ul + lug ug ul
r
D
2
The first term of the force represents the drag force, and the second term is the apparent mass force due to the change in kinetic
energy induced by the bubble motion in the liquid.
In the above expression CFL is the interfacial friction factor and
is given by
CFL = CD1 1.7
l D
g DB
Here, PL0 is the single phase pressure drop due to friction and
n = 0.25 for air water system. In Eq. 7 the empirical constant B is
dependent on the fluid pair associated. It can be calculated as
follows:
B=
where
C=
ul
ug
CY c 22n + 2
Y 2c 1
u2gg
l
1 2
g
ul l
Xc =
1+
1 l ul
g ug
and
Yc =
f G0l
f L0g
0.5
The wall friction factor for the gaseous phase is also calculated
using Eq. 7 after employing the gaseous phase frictional pressure drop PG0 as the multiplier in place of PL0.
The turbulent dispersion force of gaseous phase is considered in
the form of Favre averaged variables 36
Fdisp = CTDCD
tg l g
Sctg l
g
10
nr,z,t;d
nr,z,t;d
nr,z,t;d
+ wg
= BBr,z,t;d
+ ug
r
z
t
DBr,z,t;d + BCr,z,t;d DCr,z,t;d
11
In the present model, the entire population of bubble is discretized into different subgroups of equal volume interval. Each
subgroup has its own span of bubble characteristic length. The
mean of this span is considered to be the pivotal length scale
associated with the subgroup. As the bubbles are equispaced in
volume, a newly born bubble during coalescence explicitly falls in
one of the subgroup nodes. This exactly satisfies the conservation
of bubble volume. On the other hand, during breakage the daughter bubbles may not have a size corresponding to the pivotal
length scale. Such bubbles are distributed among the neighboring
subgroups, keeping the bubble mass and number fixed. Nucleation
of new bubbles, evaporation, condensation, shrinkage, and elongation are not considered in the present model. Present model
only recognizes the appearance of newly born bubbles and the
disappearance of the extinct bubbles due to dynamic interactions
in a fixed control volume. However, their redistribution is considered to be random and solely guided by the advection of two
phase mixture. Furthermore, only binary coalescence and binary
breakage of the bubbles are considered.
Collision with turbulent eddies are considered as the mechanism of breakup and is schematically shown in Fig. 2a. Breakage of a bubble is possible only when the turbulence kinetic energy of the striking eddy supersedes the effect of surface energy of
the interacting bubble 38. This leads to the surface rupture of the
bubbles and the formation of daughter bubbles. Moreover,
breakup occurs when an eddy of size comparable to the bubble
characteristic length collides with it. If the eddy size is larger than
the bubble characteristic length, collision between them is unable
to make any craters on the bubble. Thus, in Fig. 2a only eddies
denoted by 1, 2, and 3 are effective for the breakup process. The
size of the daughter bubbles due to breakage depends on the
strength of the eddies 39, which, in turn, depend on the local
hydrodynamics.
Birth and death of bubbles due to the breakup process can be
calculated as follows:
BBr,z,t;d =
d d,ddgdnr,z,t;ddd 12
DBr,z,t;d = nr,z,t;dgd
13
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
gd = k1
2
d
1/3
min
3
f vU m
2D
15
Um is calculated as follows:
Um = u2g + w2g2 + ul2 + wl22
16
3/5
dmin = 1.15 3/5 2/5
l
18
19
d dcm + dm
c
d1,d2 =
1
1
2z 1 6I
+
+
d1
b + 0.5 d32
d1
+a 1 +b
d2
d2
20
0.5
z1
ln1 + a lnb +
b + 0.5
and
z=
a
4b1 + b1 a
a and b are parameters that define the shape of the daughter drop
size distribution function. In the present model, the values a
= 0.1 and b = 1 45 are used that signifies U shaped bubble size
distribution.
For coalescence of bubbles to occur in the turbulent flow field
of a gas-liquid mixture, the bubbles must first collide with each
other and then remain in contact for sufficient time so that the
processes of liquid film drainage, film rupture, and finally coalescence may occur. After the collision between two bubbles, a successful coalescence through film rupture is depicted in Fig. 2b.
The following expressions can be taken for the birth and death
due to coalescence process, respectively
BCr,z,t;d =
1
2
vd/2
dvv,dvnr,z,t;dvvnr,z,t;dvdv
21
dv,dvnr,z,t;dvdv 22
17
gd = b
DCr,z,t;d = nr,z,t;dv
hd1,d2 = c2
2/3 1/2
d1 + d22d2/3
1 + d2
1+
23
pd1,d2 = exp
d 1d 2
c 3 l l
21 + 3 d1 + d2
24
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
nr,0,t;d =
for d din
25
and
nr,z,0;d = 0
v.
26
where din is the bubble size entering at the inlet plane and
can be selected as desired. This mimics a situation where
the gas phase is introduced into the conduit by a large
number of nozzles uniformly distributed at the inlet section.
Overall possible bubble sizes min d / 2 D / 2 are divided into 40 equal volumes to simulate a realistic flow
phenomenon.
Fig. 5 Comparison of void fraction of present model and experimental data of Serizawa et al. 50
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
wall, it is clear from Fig. 5 that the model predicts the radial void
fraction distribution satisfactorily. Near the wall, an improved
modeling of the motion of the continuous phase, including the
effect of boundary layer and buffer layer, is expected to give a
better result without the stringent requirement of the finer grids
everywhere in the computational domain.
Ohnuki and Akimoto 56 used an optical probe for the measurement of void fraction and a dual sensor resistivity probe for
the measurement of bubble velocity and its fluctuation. In Fig. 6,
radial distribution of void fraction obtained by numerical simulation is compared with the results of Ohnuki and Akimoto 56 for
identical conditions liquid velocity of 0.18 m/s and gas velocity
of 0.11 m/s. A satisfactory agreement has been observed. Estimated void fraction profile for a higher liquid flow rate of 0.35
m/s keeping the air flow rate the same 0.11 m/s is also depicted
in the same figure. The same phase velocities were used by
Ohnuki and Akimoto 56. A good agreement has also been obtained in this case.
Comparisons of the predicted radial profile of volume averaged
bubble diameter with experiment at same initial situations 56 are
depicted in Fig. 7. The developed model predicts the correct trend
though it underpredicts the average bubble diameter toward the
tube wall. Present model is also used to investigate the distribution of bubble Sauter mean diameter in a radial plane. When the
surface area and the volume of the averaged diameter bubbles are
equal to those of the polydisperse bubbles, this averaged diameter
is called the Sauter diameter volume-surface area mean diam-
Fig. 7 Radial distribution of volume averaged bubble diameter; a comparison between present model and experimental
data of Ohnuki and Akimoto 51 at high flow rates
1
fi
i=1 di
N
27
It has been mentioned earlier that the gas phase at the inlet
boundary is considered as a uniform distribution of equal sized
bubbles. Therefore, it is obvious that any transformation from
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
group has been entered. Moreover, a wide gap between the maximum possible bubble sizes and the bubbles to be entered will not
bias the transition criteria. There is bubble growth due to coalescence; rather smaller size bubbles are formed as one moves toward the downstream. This is clear from the bubble count histograms estimated at 10D Fig. 10, 30D Fig. 11, and 60D Fig.
12. The number of tenth subgroup bubbles introduced at the
inlet reduces, while the number of bubbles belonging to smaller
subgroups increases gradually. The coalescence of bubbles is
clearly discernable from Figs. 1315 at a high phase velocities
liquid velocity of 3 m/s and gas velocity of 1 m/s even at an
axial length of 10D. At this position one can see formation of
bubble up to a subgroup of 30 that is considerably larger compared with the bubbles introduced at the inlet plane tenth subgroup. The number of bubbles just double and triple of inlet
bubble size group is high compared with the other higher sized
bubble group as there is a rich tendency of coalescence of two
inlet sized bubbles. At further downstream the number of large
size bubbles Vi V p increases and new bubbles with even larger
diameter come into existence.
To confirm the generation of larger bubbles at higher gas and
liquid flow rates, a total volume count of bubbles above and below
those introduced at the inlet are depicted in Table 1. The generation of larger bubbles at higher phase flow rates can also be confirmed by taking an account of the volume of the bubbles of
different subgroups. For these a nondimensional parameter,
namely, Vhigh has been defined as follows:
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
v in i
Vhigh =
i=p+1
28
v in i
i=1
Table 1 Liquid velocity 0.3 m/s and gas velocity 0.1 m/s
Axial
location
10D
30D
60D
Vhigh
0.00004588
0.00004689
0.00004701
0.002962
0.00296
0.00298
0.01525
0.01559
0.015526
Vhigh
0.000667
0.001399
0.001895
0.002331
0.001609
0.001104
0.2225
0.4663
0.6317
Fig. 17 Evolution of breakage and coalescence with axial location at lower flow rates
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 18 Evolution of breakage and coalescence with axial location at higher flow rates
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Fig. 21 Transition of bubbly flow; present simulation and experimental results of Olmos et al. 27 and Mercadier 62
Fig. 22 Radial bubble count histogram for the high liquid flow
rate showing dispersed bubbly flow
cU2c D P
Dsurf
= 0.55
DP
0.6
m f
c1
0.4
29
Conclusion
Fig. 23 Indication of different flow regimes by the bubble volume histogram obtained with a variation in liquid velocity for a
fixed gas flow rate
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
Nomenclature
b
B
BB
BC
c1
c2
c3
CD
CFL
CTD
dc
DB
DC
D
d
f Go
f Lo
fv
Fdispr
Fdispz
FLGz
FLGr
FGIz
FGIr
FLIz
FLIr
References
1 Zuber, N., and Findlay, J. A., 1965, Average Volumetric Concentration in
Two Phase Flow Systems, ASME J. Heat Transfer, 87, pp. 453468.
2 Ishii, M., and Zuber, N., 1970, Thermally Induced Flow Instabilities in TwoPhase Mixtures, Proceedings of the Fourth International Heat Transfer Conference, Paris, France.
3 Saha, P., and Zuber, N., 1978, An Analytical Study of the Thermally Induced
Two Phase Flow Instabilities Including the Effects of Thermal NonEquilibrium, Int. J. Heat Mass Transfer, 21, pp. 415426.
4 Kim, C., and Roy, R. P., 1981, Two Phase Flow Dynamics by a Five Equation
Drift Flux Model, Lett. Heat Mass Transfer, 8, pp. 5768.
5 Hibiki, T., and Ishii, M., 2000, Two-Group Interfacial Area Transport Equations at Bubbly-to-Slug Flow Transition, Nucl. Eng. Des., 202, pp. 3976.
6 Rakhmatulin, K. A., 1956, Fundamentals of Gas Dynamics of Interpenetrating Motions of Compressible Media, Prikl. Mat. Mekh., 202, pp. 184195.
7 Ishii, M., 1977, One Dimensional Drift Flux Model and Constitutive Equations for Relative Motion Between Phases in Various Two-Phase Flow Regimes, Argonne National Laboratory, Technical Report No. ANL-77-47.
8 Wu, Q., Ishii, M., and Uhle, J., 1998, Frame Work of Two-Group Model for
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
46294637.
39 Luo, H., and Svendsen, H. F., 1996, Theoretical Model for Drop and Bubble
Breakup in Turbulent Dispersions, AIChE J., 42, pp. 12251233.
40 Kostoglou, M., and Karabelas, A. J., 2005, Towards a Unified Framework for
the Derivation of Breakage Functions Based on the Statistical Theory of Turbulence, Chem. Eng. Sci., 60, pp. 65846595.
41 Ioannou, K., Hu, B., Matar, O. K., Hewitt, G. F., and Angeli, P., 2004, Phase
Inversion in Dispersed LiquidLiquid Pipe Flows, Proceedings of the Fifth
International Conference on Multiphase Flow, Yokohama, Japan.
42 Tsouris, C., and Tavlarides, L. L., 1994, Breakage and Coalescence Models
for Drops in Turbulent Dispersions, AIChE J., 40, pp. 395406.
43 Troshko, A. A., and Zdravistch, F., 2009, CFD Modeling of Slurry Bubble
Column Reactors for FisherTropsch Synthesis, Chem. Eng. Sci., 64, pp.
892903.
44 Carrica, P. M., and Clausse, A. A., 1993, Mathematical Description of the
Critical Heat Flux as Nonlinear Dynamic Instability, Instabilities in Multiphase Flow, G. Gouesbet and A. Berlemont, eds., Plenum, New York.
45 Kostoglou, M., and Karabelas, A. J., 1998, Theoretical Analysis of Steady
State Particle Size Distribution in Limited Breakage Process, J. Phys. A, 31,
pp. 89058921.
46 Tomiyama, A., Nakahara, Y., Adachi, Y., and Hosokawa, S., 2003, Shapes
and Rising Velocities of Single Bubbles Rising Through an Inner Subchannel,
J. Nucl. Sci. Technol., 40, pp. 136142.
47 Coulaloglou, C. A., and Tavlarides, L. L., 1977, Description of Interaction
Processes in Agitated Liquid-Liquid Dispersions, Chem. Eng. Sci., 32, pp.
12891297.
48 Chesters, A. K., 1991, The Modelling of Coalescence Processes in FluidLiquid Dispersions: A Review of Current Understanding, Trans. Inst. Chem.
Eng., 69, pp. 259270.
49 Lovick, J., 2004, Horizontal Oil-Water Flows in the Dual Continuous Flow
Regime, Ph.D. thesis, University College London, England.
50 Hu, H. G., and Zhang, C., 2007, A Modified k Turbulence Model for the
Simulation of Two-Phase Flow and Heat Transfer in Condensers, Int. J. Heat
Mass Transfer, 50, pp. 16411648.
51 Ekambara, K., Sanders, R. S., Nandakumar, K., and Masliyah, J. H., 2008,
CFD Simulation of Bubbly Two-Phase Flow in Horizontal Pipes, Chem.
Eng. J., 144, pp. 277288.
52 Prosperetti, A., and Tryggvason, G., 2007, Introduction: A Computational
Approach to Multiphase Flow, A. Prosperetti and G. Tryggvason, eds., Computational Methods for Multiphase Flow, Cambridge University Press, Cambridge, England.
53 Tomiyama, A., and Shimada, N., 2001, A Numerical Method for Bubbly
Flow Simulation Based on a Multi-Fluid Model, ASME J. Pressure Vessel
Technol., 123, pp. 510516.
54 Schumann, U., 1975, Subgrid Scale Model for Finite Difference Simulations
of Turbulent Flows in Plane Channels and Annuli, J. Comput. Phys., 18, pp.
376404.
55 Serizawa, A., Kataoka, I., and Michiyoshi, I., 1975, Turbulence Structure of
Air-Water Bubbly FlowII. Local Properties, Int. J. Multiphase Flow, 2, pp.
235246.
56 Ohnuki, A., and Akimoto, H., 2000, Experimental Study on Transition of
Flow Pattern and Phase Distribution in Upward Air-Water Two-Phase Flow
Along a Large Vertical Pipe, Int. J. Multiphase Flow, 263, pp. 367386.
57 Shen, X., Mishima, K., and Nakamura, H., 2005, Two-Phase Phase Distribution in a Vertical Large Diameter Pipe, Int. J. Heat Mass Transfer, 481, pp.
211225.
58 Nakoryakov, V. E., Kashinsky, O. N., Randin, V. V., and Timkin, L. S., 1996,
Gas Liquid Bubbly Flow in Vertical Pipes. Data Bank Contribution, ASME
J. Fluids Eng., 118, pp. 377382.
59 Lucas, D., Krepper, E., and Prasser, H. M., 2005, Development of Co-Current
Air-Water Flow in a Vertical Pipe, Int. J. Multiphase Flow, 31, pp. 1304
1328.
60 Radovicich, N. A., and Moissis, R., 1962, The Transition From Twophase
Bubble Flow to Slug Flow, MIT Technical Report No. 7-7633-22.
61 Bilicki, Z., and Kestin, J., 1987, Transition Criteria for Two Phase Flow
Patterns in Vertical Upward Flow, Int. J. Multiphase Flow, 13, pp. 283294.
62 Mercadier, Y., 1981, Contribution aletude des propagations de perturbations
de taux de vide dans les ecoulements diphasiques eau-air a bulles, Ph.D.
thesis, Institut National Polytechnique de Grenoble, Universite Scientifique et
Medicale, France.
63 Matuszkiewicz, A., Flamand, J. C., and Boure, J. A., 1987, The Bubble Slug
Flow Pattern Transition and Instabilities of Void Fraction Waves, Int. J. Multiphase Flow, 13, pp. 199217.
64 Sun, B., Wang, R., Zhao, X., and Yan, D., 2002, The Mechanism for the
Formation of Slug Flow in Vertical Gas-Liquid Two Phase Flow, Solid-State
Electron., 46, pp. 23232329.
65 Krussenberg, A. K., Prasser, H. M., and Schaffrath, A., 1999, A New Criterion for the Bubble Slug Transition in Vertical Tubes, Proceedings of the
Ninth International Topical Meeting on Nuclear Reactor Thermal Hydraulics
NURETH-9, San Francisco, CA.
66 Hinze, J., 1955, Turbulence, McGraw-Hill, New York.
67 Brauner, N., 2001, The Prediction of Dispersed Flows Boundaries in LiquidLiquid and Gas-Liquid Systems, Int. J. Multiphase Flow, 27, pp. 885910.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
A. K. Das
e-mail: arup@mech.iitkgp.ernet.in
P. K. Das
e-mail: pkd@mech.iitkgp.ernet.in
Department of Mechanical Engineering,
IIT Kharagpur,
721302, India
J. R. Thome
LCTM,
EPFL,
Lausanne CH-1015, Switzerland
e-mail: john.thome@epfl.ch
Introduction
bution with groups of constant mass. Their model can predict the
presence of wall and core peakings but remains silent about the
possibility of any other voidage profile. Carrica and Clausse 10
solved the behavior of polydispersed bubble formation using a
multigroup model along with a two fluid model. Recently,
Kreeper et al. 11 used the N-M MUSIG model for solving bubbly flow having N velocity group and M bubble group size. Yeoh
and Tu 12 solved the population balance equation in axial direction along with a two fluid model using the method of discrete
classes in order to incorporate the effect of bubble size.
In a companion paper 13, using the two-fluid model and the
population balance technique, the hydrodynamics of bubbly flow
has been simulated. Emphasis has been given on the analysis of
axial flow development. The population balance technique enables
one to detect the bubble population at any axial level and to assess
the process of coalescence and breakage. An objective and comprehensive criterion for the transition from bubbly to slug flow has
been proposed. The simulation shows excellent agreement with
published results.
Certainly the importance of axial development cannot be overlooked as it has direct bearing to the transition. Nevertheless, hydrodynamics of bubbly flow strongly depends also on the phase
distribution and bubble population along the radial direction. This
gives rise to different peaked structure 1,3,14 of the void fraction and also influences the transition. To date the prediction of
bubbly-slug transition based on superficial phase velocities
15,16 is well accepted by many researchers. However, such a
criterion does not take into cognition the effect of inlet bubble size
17 and tube diameter 18 on the transition from bubbly to slug
flow. In the present paper, the basic two fluid and population balance model 13 for bubbly flow has been employed to investigate
the above issues. In the first part of the paper, we have made a
thorough study on the radial distribution of dispersed phase. In the
second part, the effect of operating variables other than the phase
velocities on the flow regime transition has been investigated.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
In a number of experiments on gas liquid upflow through vertical conduits, nonuniform void distribution has been observed at
the downstream while gas bubbles were introduced uniformly at
the inlet plane. Theoretical simulations of the bubbly flow also
corroborate the peaked characteristic of the void profile 15,19. It
is important to note that not only the absolute value of the void
fraction varies in the different flow regimes, each of the regimes
has a typical pattern for radial distribution of the void fraction.
Therefore, in the present paper, we want to investigate further into
the transformation of the void profile due to the operating and
inlet conditions.
The present model is compared with the experimental void
fraction obtained by Zun 3, Song et al. 4, and Lucas et al. 20
in the case of a two-phase flow mixture of air and water in vertical
ducts. To simulate the experiments, the ducts were discretized into
40 radial nodes with finer discretization near the wall of the tube.
In the axial direction 500 uniform spaced grids were generated to
track the evolution of the bubbly flow. Bubble mass is divided into
40 equal groups by volume to eradicate the complexity of the
volume discretization during the coalescence of bubbles. In case
of bubble breakage, if the daughter bubbles do not fall in an exact
size group they are distributed among its immediate neighbors,
maintaining the conservation of mass and number.
An extensive study of the bubbly flow 20 has been made in
the MTLoop facility. The test section in this loop has an inner
diameter of 52.3 mm and a length of 3.5 m. Wire mesh sensors
have been used for the measurement of void fractions at a distance
of 60D from the inlet. Figure 1a shows the simulation results at
an air flow rate of 0.0096 m/s and a water flow rate of 0.0405 m/s
as has been used by Lucas et al. 20. Most of the bubbles are
gathered near the wall of the tube. Results from the present model
also show a similar pattern of the void fraction distribution. Similar experimental studies were made by Prasser et al. 21. They
used a larger tube with 195.3 mm diameter and 8.78 m length. At
an air flow rate of 0.53 m/s and water flow rate of 1.02 m/s, they
observed that the void distribution pattern changed Fig. 2 into a
parabolic shape exhibiting core peak. Simulation results for the
same conditions satisfactorily agree Fig. 2 with the experimental
observations. This shows that the present model can handle both
inward and outward lateral migrations of the bubbles, depending
on the flow velocities and tube diameter.
Zun 3 proposed that uniform input of big bubbles at the inlet
produces a core peak while smaller bubbles at the inlet produce a
wall peak. Interestingly, when a mixture of big and small bubbles
091304-2 / Vol. 131, SEPTEMBER 2009
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
rate is kept at 0.079 m/s and air flow rate is maintained at 0.0082
m/s, which are identical to the experimental conditions of Song et
al. 4. Generated results at an axial length of 60D were presented
to minimize the entrance effect. Agreement between the simulation and experiment is reasonably well.
Occurrence of different peaked structures of the voidage profile
in the bubbly flow has been observed in different independent
experiments 1. The reason for this has loosely been attributed to
the forces acting on individual bubbles and the interactions between neighboring bubbles. As the evolution of the voidage profile is a result of the interplay of various complex phenomena, its
prediction is not possible through a control volume based model
prevalent in the analysis of the two-phase flow. On the contrary,
a computational fluid dynamics CFD simulation with an inbuilt
capability to model the bubble interactions can simulate this phenomenon reasonably well. The capability of the present model for
predicting bubbly-slug transition has been presented in an earlier
communication 13. The model prediction depicts a good match
with the regime boundary proposed by Taitel et al. 22 experimentally. It is interesting to note that a unique regime boundary
was not observed by different investigators. In Fig. 5 bubbly-slug
transitions as observed in different experiments point out a substantial difference between themselves. While the trend of the
transition boundary of Griffith and Wallis 23, Dukler and Taitel
15, and Mishima and Ishii 24 are similar which also matches
with the prediction of the present simulation, the nature of the
other transition boundaries 22,2528 are different. It may be
quite possible that the difference in the inlet condition of the spe-
cific experiments is mainly responsible for this mismatch. Unfortunately, in most of the experimental results the inlet conditions,
particularly the bubble size distribution, have not been specified
explicitly. Therefore, it is difficult to substantiate the above proposition directly. As an alternative, we have tried to predict the transition boundary for different inlet conditions. The results are discussed below.
To see the effect of bubble size at the inlet on the flow pattern
map, results were generated using the present model for homogeneous inputs of various bubble sizes. Figure 6 shows the change in
the flow regime map due to the change in the bubble size at the
inlet from 0.5 mm to 1 mm and then to 3 mm. Bubbly flow turns
into slug flow at a lower flow rate of liquid for a fixed gas flow
rate with the increase in inlet bubble size. It is obvious that
bubbles of higher size at the input takes less effort to merge into a
big gaseous slug. Results generated from the developed model
efficiently support the intuitive concept about bubble coalescence.
Flow pattern map of Bilicki and Kestin 29 for different sized
bubble inputs are also plotted in Fig. 6 to support the results
generated from the present model. A plexiglass tube with 20 mm
diameter and 1.5 m length has been used in Ref. 29 to view the
two-phase flow patterns for air and water. Tests were carried out
for bubble diameters of 1 mm and 3 mm. The observations are in
close concordance with the simulation results.
In most of the experiments monosized bubbles are rarely injected at the inlet. In general there could be a distribution of
bubble size. Simulations are made to study the effect of a mixed
bubble population at the inlet on the flow regime transition. Different inlet populations consisting mainly of a 0.5 mm bubble
diameter, b 1 mm bubble diameter and c 2 mm bubble diameter as well as d 50% bubbles of 0.5 mm diameter+ 50% of 1
mm diameter, and e 50% bubbles of 1 mm diameter+ 50% of 2
mm diameter were considered at the inlet, as shown in Fig. 7. In
all the cases, the inlet volume fraction is kept constant. It can be
seen that with the mixing of higher sized bubbles the flow pattern
map for 1 mm diameter bubbles at the inlet shifts downward,
causing the transition at a lower air flow rate. Similarly, when 0.5
mm diameter bubbles are mixed with 1 mm diameter bubbles,
slug bubbles appear at a higher gas flow rate for a fixed liquid
flow rate compared with the homogeneous input of 1 mm diameter bubbles. This analysis can be extended further to investigate
the effect of polydispersed bubble input at the inlet plane on the
void distribution profiles and flow pattern map. It is also clear
from the results generated from the simulation that along with the
superficial velocities, the bubble diameter is an important parameter for the transition of bubbly flow to slug flow. It is interesting
to note that curves a and d, and curves c and e intersect each other,
indicating an implicit effect of the inlet bubble size on the transition boundary. Effect of mixing two different bubble sizes cannot
SEPTEMBER 2009, Vol. 131 / 091304-3
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
be predicted even qualitatively from the logic of linear averaging. The input bubble diameter distribution influences the hydrodynamics in a number of ways. It has got direct effects on the rate
of breakage and coalescence as well as on the drag force. The
laws governing these phenomena are strongly nonlinear. The regime boundaries depicted in Fig. 7 is a reflection of that fact.
It is a well recognized fact that flow regime transition is not a
discrete event. It is more likely that a fully developed bubbly flow
transforms into a fully developed slug flow through a gradual
change. In this regard the different peaked structures of the voidage distributions in bubbly flow may have a role to play. The
present model has been used to investigate this. For a case study,
uniform distribution of bubbles of 0.05D has been considered at
the inlet plane and the void fraction profile is investigated at 60D
from the inlet. Figure 8 shows the boundary obtained between the
wall and core peaking void fraction profiles based on the liquid
and gaseous phase flow velocities. To make a comparison with the
available void distribution, map results of Serizawa and Kataoka
30, Ohnuki and Akimoto 14, and Lucas et al. 31 are depicted
in the same figure. Serizawa and Kataoka 30 proposed a band of
flow velocities beyond which a void fraction profile with core
peaking changes into a profile with wall peaking. Present results
of the simulation falls well in the band proposed by Serizawa and
Kataoka 30. The boundaries proposed by Ohnuki and Akimoto
14 and Lucas et al. 31, though fall within the band, do not
match well with the present simulation. In their model, Lucas et
al. 31 assumed that the whole bubble mass is divided into three
separate bubble groups, whereas in the present simulation 40 subgroups have been considered. This may be one of the reasons in
the discrepancy of the results.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
3
60 grids have been taken for the 76.2 mm diameter tube. Axial
length is kept constant for all the tube diameters at 5 m. All the
results depicted in the present paper are taken at a distance of 40D
from the inlet plane. Results are depicted in Figs. 1012. Void
fraction distribution at an air flow rate of 0.1 m/s and a liquid flow
rate of 1 m/s is depicted in Fig. 10 for all the three tube sizes. In
all the cases, sharp peak can be visible at the wall of the conduit.
At a 25.4 mm diameter the rise of the peak is highest at the
vicinity of the wall while the deepest valley is observed at the
central region. With the increase in the diameter the peak sharpness decreases, making the voidage profile relatively uniform. The
position of the peak near the wall also shifts toward the center of
the tube geometry. Results of the void fraction distribution for
various tube diameters at higher liquid 3 m/s and air 2 m/s flow
rates are depicted in Fig. 11. In all the cases, core peaking is
observed. However, the profile becomes flatter with the increase
in the tube diameter. The figures bring out the diminishing wall
effect on the gaseous phase with the increase in the tube diameter.
As the change in tube diameter causes a change in the voidage
profile, it is not unwise to expect that the tube diameters may have
an effect on the flow regime transition. Flow regime boundary has
been constructed using the present simulation for the three different tube diameters, as shown in Fig. 12. There is a distinct upward
shift of the transition line with the increase in the conduit size.
The data available for the 25.4 mm 18 and 50.8 mm 22 tube
diameters exhibit reasonable agreement with the simulation results.
Conclusion
References
1 Serizawa, A., Kataoka, I., and Michiyoshi, I., 1975, Turbulence Structure of
Air-Water Bubbly FlowII Local Properties, Int. J. Multiphase Flow, 2, pp.
235246.
2 Sekoguchi, K., Fukui, H., and Sato, Y., 1981, Flow Characteristics and Heat
Transfer in Vertical Bubble Flow, Proceedings of Japan-U.S. Seminar on
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
19 Cheung, S. C. P., Yeoh, G. H., and Tu, J. Y., 2007, On the Modelling of
Population Balance in Isothermal Vertical Bubbly FlowsAverage Bubble
Number Density Approach, Chem. Eng. Process., 46, pp. 742756.
20 Lucas, D., Krepper, E., and Prasser, H. M., 2007, Use of Models for Lift,
Wall and Turbulent Dispersion Forces Acting on Bubbles for Poly-Disperse
Flows, Chem. Eng. Sci., 6215, pp. 41464157.
21 Prasser, H. M., Beyer, M., Carl, H., Gregor, S., Lucas, D., Pietruske, H.,
Schtz, P., and Weiss, F. P., 2007, Evolution of the Structure of a Gas-Liquid
Two Phase Flow in a Large Vertical Pipe, Nucl. Eng. Des., 237, pp. 1848
1861.
22 Taitel, Y., Bornea, D., and Dukler, A. E., 1980, Modelling Flow Pattern
Transitions for Steady Upward Gas-Liquid Flow in Vertical Tubes, AIChE J.,
263, pp. 345354.
23 Griffith, P., and Wallis, G. B., 1961, Two-Phase Slug Flow, ASME J. Heat
Transfer, 83, pp. 307320.
24 Mishima, K., and Ishii, M., 1984, Two-Fluid Model and Hydrodynamic Constitutive Relations, Nucl. Eng. Des., 8223, pp. 107126.
25 Sternling, V. C., 1965, Two-Phase Theory and Engineering Decision, Award
Lecture, AIChE Annual Meeting.
26 Govier, G. W., and Aziz, K., 1972, The Flow of Complex Mixtures in Pipes,
Van Nostrand Reinhold, New York, pp. 388389.
27 Gould, T. L., 1974, Vertical Two-Phase Steam Water Flow in Geothermal
Wells, J. Pet. Technol., 26, pp. 833842.
28 Oshinowo, T., and Charles, M. E., 1974, Vertical Two-Phase Flow. Part 1:
Flow Pattern Correlations, Can. J. Chem. Eng., 52, pp. 2535.
29 Bilicki, Z., and Kestin, J., 1987, Transition Criteria for Two Phase Flow
Patterns in Vertical Upward Flow, Int. J. Multiphase Flow, 13, pp. 283294.
30 Serizawa, A., and Kataoka, I., 1987, Phase Distribution in Two-Phase Flow,
Proceedings of ICHMT International Seminar on Transient Two-Phase Flow,
Dubrovnik, Yugoslavia, pp. 2430.
31 Lucas, D., Krepper, E., and Prasser, H. M., 2005, Development of Co-Current
Air-Water Flow in a Vertical Pipe, Int. J. Multiphase Flow, 31, pp. 1304
1328.
32 Rao, N. M., 2002, Investigations on Buoyancy Induced Circulation Loops,
Ph.D. thesis, IIT, Kharagpur, India.
33 Serizawa, A., and Kataoka, I., 1988, Phase Distribution in Two Phase Flow,
Transient Phenomena in Multiphase Flow, N. Afgan ed., Hemisphere, Washington, DC, pp. 179224.
34 Hibiki, T., and Ishii, M., 2000, Two-Group Interfacial Area Transport Equations at Bubbly-to-Slug Flow Transition, Nucl. Eng. Des., 202, pp. 3976.
35 Dukler, A. E., and Taitel, Y., 1986, Flow Pattern Transitions in Gas-Liquid
Systems: Measurement and Modeling, Multiphase Science and Technology,
Vol. 2, G. F. Hewitt, J. M. Delhaye, and N. Zuber, eds., Hemisphere, Washington, DC, pp. 194.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
M. Firouzi
e-mail: hashemabadi@iust.ac.ir
Computational Fluid Dynamics Research Laboratory,
School of Chemical Engineering,
Iran University of Science and Technology,
Narmak, Tehran 16846, Iran
In this paper, the motion equation for steady state, laminar, fully
developed flow of Newtonian fluid through the concave and convex ducts has been solved both numerically and analytically.
These cross sections can be formed due to the sedimentation of
heavy components such as sand, wax, debris, and corrosion products in pipe flows. The influence of duct cross section on dimensionless velocity profile, dimensionless pressure drop, and friction
factor has been reported. Finally based on the analytical solutions
three new correlations have been proposed for the product of
Reynolds number and Fanning friction factor C f Re for these
geometries. DOI: 10.1115/1.3184026
Keywords: analytical, numerical, Newtonian, concave, convex,
noncircular passage, Fanning friction factor, sedimentation
Introduction
A range of noncircular duct configurations is employed for internal flow and heat transfer applications in chemical, petroleum,
pharmaceutical, food, and plastics industries. It is of great practical interest to be able to predict the pressure drop and heat transfer
characteristics of fluid flows in different duct geometries. One of
the applications of fluid flow in irregular ducts is in the piping
systems that encountered with sedimentation, and most fluid flows
in environmental and industrial applications are directly influenced by sediment. Sediment transport is important to river,
shoreline, and harbor projects in piping systems. It can cause to
form geometric domains that are far from ideal in shape. It is not
difficult to imagine subsea pipelines blocked by accumulated wax
or by hydrate plugs, bundled pipes with debris settlement, or
heavily clogged cross sections with grease or hard water scaling
Fig. 1 1. Few exact solutions are available only for simple fluid
models and circular pipe cross sections. Many attempts were
made in modeling Newtonian and non-Newtonian fluids in various noncircular cross sectional ducts, which most of them are
numerical or experimental. Several investigations have been reported for flow characteristics in trapezoidal cross section 24,
in rectangular cross section 57, non-Newtonian fluid flow
through triangular channels 8,9, friction factor, pressure drop,
and Nusselt number of Newtonian fluid flowing in ducts of arbitrary cross sections 10, and an approximate solution for determining the pressure drop in microchannels of arbitrary cross section 11.
Here we focused on analytical solution for laminar fully developed Newtonian fluid flow through a family of noncircular passage with concave and convex cross sections, which is selected
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received May 20, 2008; final manuscript
received June 21, 2009; published online August 14, 2009. Assoc. Editor: Ye Zhou.
Mathematical Modeling
2.1
S. H. Hashemabadi
Analytical Solution
R sin 1 sinh
,
cosh cos
y=
R sin 1 sin
cosh cos
cosh cos
sin 1
2v 2v
+
=4
2 2
2 sin 1sin1
cosh cos
2vh 2vh
+
=0
2
2
vh=2 +
vh=1 = 0
vh= = 0
2 sin2 1sin1
=0
cosh cos 2
sinh 1cosd
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
2cosh cos
cosd = 1
sinh
14
15
2K1 + cosh
sinhsinh 1
16
sinh2 1cosd +
2 sin2 1sin1
=0
cosh cos 2
v , =
cosd = 1
11
where
= sinh2 1 ,
K =
sinh 1cosd
+ 2 sin1
10
sin 1
cosh cos
17
2 sin2 1sin1
cosh cos 2
Q=
12
v,
R sin 1
cosh cos
dd
R4 sin2 1 dp
F1, 2
dz
4
18
19
where
+
F1, 2 =
v ,
dd
cosh cos 2
The flow pressure drop through the noncircular ducts can be calculated from Eq. 18:
4Q
dp
F1, 2
=
dz R4 sin2 1
20
dp
dp
dz
=
=
dz
dp
2 sin21F1, 2
dz pipe
21
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
z =
1 cosh cos v,
2
sin 1
22
C f Re = 4
sin sin 1
cosh cos
cosh cos
2 sin 1
cosh 1d 23
v =
D2h
R
z , 1
Dh =
z, 2
2R 1 + 1/2 sin21
1 + sin1
28
1
2 + sin22
2
sin 1
1 + 2
sin 2
v ,
z, 2
24
R sin 1
d
cosh cos 1
R sin 1
d
cosh cos 2
R sin 1
cosh cos
R sin 1
d
cosh cos
dd
v =
25
dd
26
v
,
R2 dp
4 dz
x =
x
,
D
y =
27
y
D
32
sin 1
1
2R 1 + sin21
2
sin 2
1
A
R sin 1
v ,
cosh cos
Dh
R
where A and P are the area and the perimeter of the duct cross
section, respectively. By substituting Eqs. 25 and 26 into Eq.
24, it yields the following:
R sin 1
d +
cosh cos
z, 1
C f Re =
1
P
+
Then the friction factor, which is one of the most important characteristics in fluid flows, can be calculated. The definition of Fanning friction factor is 13 as follows:
w
v
z = cos 1
29
2.2 Numerical Solution. In this section, to recheck the results
of the analytical solution because of insufficiency of related established results for more validation of obtained analytical results,
the momentum equation is solved numerically with finite element
method FEM. The momentum equation is preferred to be solved
in the Cartesian coordinate system, which, by considering the
mentioned assumptions in analytical solution, can be simplified as
follows:
zx zy dp
+
+
=0
x
y
dz
30
v
v
= 16
+
x x
y x
where
Journal of Fluids Engineering
31
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
= 0.5
D
1 cos 2
1
tan
2
33
Fig. 7 Variation in the dimensionless pressure drop with aspect ratio h / D for various duct cross sections different 1
and 2
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
C f Re = 16 1 +
n
n 1
n=1
C f Re = 16 1 +
n
n 1
n=1
34
35
Conclusion
Nomenclature
A
Cf
D
Dh
F
h
L
p
P
Q
flow area
Fanning friction factor
pipe diameter
hydraulic diameter
defined by Eq. 19
duct height
pipe length
pressure
perimeter of the duct
flow rate
SEPTEMBER 2009, Vol. 131 / 094501-5
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm
2 =
2
2
R
Re
v
x
y
Y
z
0.180
2.170
0.191
1.576
0.216
0.356
4.970
0.096
0.003
7.467
1.449
0.200
6.021
0.252
0.043
2.673
0.136
0.112
0.616
-
0.057
-
pipe radius
Reynolds number
velocity
related to the Cartesian coordinate
related to the Cartesian coordinate
distance from the bottom wall
axial direction
Greek Symbols
dynamic viscosity
ratio of the radius vectors in the bipolar
coordinate
density
stress
view angle of the interface in the bipolar
coordinate
Subscripts
1 upper wall
2 bottom wall
max maximum
Superscripts
dimensionless property
- average value
h homogeneous
p particular
References
1 Barat, R., Bouzas, A., Mart, N., Ferrer, J., and Seco, A., 2009, Precipitation
Assessment in Wastewater Treatment Plants Operated for Biological Nutrient
Removal: A Case Study in Murcia, Spain, J. Environ. Manage., 90, pp.
850857.
2 Wu, P. Y., and Little, W. A., 1983, Measurement of Friction Factor of the
Flow of Gases in Very Fine Channels Used for Microminiaturize JouleThomson Refrigerators, Cryogenics, 248, pp. 273277.
3 Wu, H. Y., and Cheng, P., 2003, Friction Factors in Smooth Trapezoidal
Silicon Microchannels With Different Aspect Ratios, Int. J. Heat Mass Transfer, 46, pp. 25192525.
4 Weilin, Q., Mala, G. M., and Dongqing, L., 2000, Pressure-Driven Water
Flows in Trapezoidal Silicon Microchannels, Int. J. Heat Mass Transfer,
433, pp. 353364.
5 Peng, X. F., Peterson, G. P., and Wang, B. X., 1994, Frictional Flow Characteristics of Water Flowing Through Microchannels, Exp. Heat Transfer, 74,
pp. 249264.
6 Pfund, D., Rector, D., Shekarriz, A., Popescu, A., and Welty, J., 2000, Pressure Drop Measurements in a Micro Channel, AIChE J., 468, pp. 1496
1507.
7 Hrnjak, P., and Tu, X., 2007, Single Phase Pressure Drop in Microchannels,
Int. J. Heat Fluid Flow, 281, pp. 214.
8 Etemad, S. Gh., Mujumdar, A. S., and Nassef, R., 1996, Simultaneously
Developing Flow and Heat Transfer of Non-Newtonian Fluids in Equilateral
Triangular Duct, Appl. Math. Model., 2012, pp. 898908.
9 Hashemabadi, S. H., Etemad, S. Gh., Golkar Naranji, M. R., and Thibault, J.,
2003, Laminar Flow of Non-Newtonian in Right Triangular Ducts, Int.
Commun. Heat Mass Transfer, 301, pp. 5360.
10 Etemad, S. Gh., and Bakhtiari, F., 1999, General Equations for Fully Developed Fluid Flow and Heat Transfer Characteristics in Complex Geometries,
Int. Commun. Heat Mass Transfer, 262, pp. 229238.
11 Bahrami, M., Yovanovich, M., and Richard Culham, J., 2007, A Novel Solution for Pressure Drop in Singly Connected Microchannels of Arbitrary CrossSection, Int. J. Heat Mass Transfer, 501314, pp. 24922502.
12 Firouzi, M., and Hashemabadi, S. H., 2008, Analytical Solution for
Newtonian-Bingham Plastic Two-Phase Pressure Driven Stratified Flow
Through the Circular Ducts, Int. Commun. Heat Mass Transfer, 35, pp. 666
673.
13 Bird, R. B., Armstrong, R. C., and Hassager, O., 1978, Dynamics of Polymeric
Liquids, Vol. 1, Wiley, New York, pp. 592595.
14 Brauner, N., Rovinsky, J., and Maron, D. M., 1996, Analytical Solution for
Laminar-Laminar Two Phase Stratified Flow in Circular Conduit, Chem. Eng.
Sci., 141, pp. 103143.
15 Biberg, D., and Halvorsen, G., 2000, Wall and Interfacial Shear Stress in
Pressure Driven Two-Phase Laminar Stratified Pipe Flow, Int. J. Multiphase
Flow, 2610, pp. 16451673.
16 Syrjl, S., 2002, Accurate Prediction of Friction Factor and Nusselt Number
for Some Duct Flows of Power-Law Non-Newtonian Fluid, Numer. Heat
Transfer, Part A, 411, pp. 89100.
Downloaded 03 Jun 2010 to 171.66.16.158. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm