J Fluids Engineering 2009 Vol 131 N4
J Fluids Engineering 2009 Vol 131 N4
J Fluids Engineering 2009 Vol 131 N4
Fluids Engineering
Published Monthly by ASME
RESEARCH PAPERS
Flows in Complex Systems
041101
041102
041103
041104
041202
PUBLISHING STAFF
041203
041204
041205
PUBLICATIONS COMMITTEE
Chair, B. RAVANI
OFFICERS OF THE ASME
President, THOMAS M. BARLOW
Executive Director, THOMAS G. LOUGHLIN
Treasurer, T. D. PESTORIUS
Multiphase Flows
041301
041302
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
APRIL 2009
TECHNICAL BRIEFS
044501
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
M. Brian Thomas
Department of Industrial and Manufacturing
Engineering,
Cleveland State University,
2121 Euclid Avenue, SH 221
Cleveland, OH 44115-2214
e-mail: m.thomas.84@csuohio.edu
Gary P. Maul
Department of Industrial, Welding, and Systems
Engineering,
The Ohio State University,
210 Baker Systems,
1971 Neil Avenue,
Columbus, OH 43210
e-mail: maul.1@osu.edu
Considerations on a Mass-Based
System Representation of a
Pneumatic Cylinder
Pneumatic actuators can be advantageous over electromagnetic and hydraulic actuators
in many servo motion applications. The difficulty in their practical use comes from the
highly nonlinear dynamics of the actuator and control valve. Previous works have used
the cylinders position, velocity, and internal pressure as state variables in system models.
This paper replaces pressure in the state model with the mass of gas in each chamber of
the cylinder, giving a better representation of the system dynamics. Under certain circumstances, the total mass of gas in the cylinder may be assumed to be constant. This
allows development of a reduced-order system model. DOI: 10.1115/1.3089533
Keywords: pneumatics, servo motion, control
Pneumatic actuators present an alternative to the use of electromagnetic and hydraulic actuators in servo motion applications.
Like electromagnetic actuators i.e., motors, pneumatic actuators
are generally clean and reliable in operation. Like hydraulic actuators, pneumatic actuators may be directly coupled to a payload
without the need for a transmission. Servopneumatics offer high
power-to-weight ratios and can provide a cost benefit as high as
10:1 over competing technologies 1,2. In certain niche applications, pneumatic actuators provide unique advantages over competitive technologies. Unlike electromagnetic actuators, pneumatic actuators are naturally compliant due to the compressibility
of air but are well adapted to wash-down environments. The natural compliance of a pneumatic actuator also makes it advantageous when personnel must be inside the work envelope of the
machine 3.
The principal barrier to the more widespread application of servopneumatic technology in industry has been the problem of control. The dynamics of a pneumatic actuator are highly nonlinear
due to the effects of friction, air compressibility, and time delays
and are compounded by the nonlinear flow characteristics of flow
control valves 3. Linear control methods such as proportionalintegral-derivative PID, while well understood by both engineers
and technicians, perform poorly when coupled to such highly nonlinear systems. With the recent availability of low-cost powerful
microprocessors, computationally expensive nonlinear control
strategies may now be performed in real time. Advanced control
strategies applied to servopneumatics include fuzzy control 4,
sliding mode control 5, model-based adaptive control 1,6, neural networks 7, and pulse-width modulation schemes for solenoid valves 8,9.
The principal contributions of this paper are as follows:
1 Presentation of a mass-based state representation for a
pneumatic cylinder. Previous works use position, velocity,
and the pressure in two chambers of a cylinder as state
variables, as these quantities are readily measured. A massbased representation more closely links the state variables
with the process dynamics.
2 Formal development and verification of the constant mass
assumption. This assumption states that under certain conContributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received July 25, 2007; final manuscript
received December 22, 2008; published online March 6, 2009. Review conducted by
Joseph Katz.
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
Festo DNG-50200-P-A
Servovalve
Position feedback
Flow measurement
Cylinder stroke, L
Cylinder diameter, d
Rod diameter, drod
Payload mass, M
Static friction, Fstat
Dynamic friction, Fdyn
Viscous friction, b
Orientation
Control
Festo MPYE-5a
Gemco Blue Ox 952QD
None
175 mmb
50.0 mm
25 mm
11.0 kgc
21.8 N 4.90 lb
6.9 N 1.55 lb
2.19 N s / m0.125 lbs/ in.
Vertical; rod end down
PI with velocity feed-forward; implemented
through Allen-Bradley M02AE
motor control module
2.1 Conventional Representation. The conventional mathematical representation of a pneumatic actuator uses four state
variables x, x, P1, and P2 and three nonlinear differential equations to describe the dynamics of the system. Though previous
works deriving these equations differ on their underlying assumptions, they all use the general approach presented by Shearer in
Ref. 10. This approach begins with the internal energy of volume of gas,
Table 1 Experimental system properties
Cylinder
Servovalve
Position feedback
Flow measurement
Cylinder stroke, L
Cylinder diameter, d
Rod diameter, drod
Payload mass, M
Supply pressure, PS
Orientation
Control
RT
P 1A 1
1
m
x
A1x + VX,1
A1x + VX,1
RT
P 2A 2
2+
m
x
A2L x + VX,2
A2L x + VX,2
P1 =
System Models
A typical servopneumatic system consists of a pneumatic actuator, one or more control valves to regulate flow, a position sensor,
and a controller. The components may be selected from an extensive list of models and vendors, depending on the particular requirements of the application. This work considers a system consisting of a single-rod, double-acting pneumatic cylinder with a
proportional spool valve controlling the flow of air to and from
the cylinder Fig. 1. Mass flow meters measure flow in and out of
the cylinder. Table 1 details the configuration of the system. Previous work by the authors contains a simulation of a similar system 21; its configuration is detailed in Table 2.
E = cVVT
P2 =
where VX,1 and VX,2 are the excess volumes of gas in either end of
the cylinder and in the plumbing between the cylinder and the
valve. Derivation of Eqs. 2 and 3 is provided in Appendix A.
The third governing equation for the pneumatic actuator addresses
the mechanical dynamics of the system,
x = Mg + P1A1 P2A2 PatmArod bx Ffric/M
RT
P
m
V
P =
V
V
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
ences between an isothermal and an adiabatic model are not significant in their simulation and recommended using the mathematically less-complicated isothermal model.
The state equationsEqs. 24do not address the compu1
tation of the mass flow rates into each chamber of the cylinder, m
2. These flow rates are dependent on the design, construcand m
tion, and operation of the valve or valves, the upstream and
=
m
C DA
Pup
RT1
2
1
C DA
2/
Pdn
Pup
Pup
RT1
0.6847
Pdn
Pup
2
+1
crit
/1
QSCFM =
22.48 CV
Pdn
Pdn
Pup PdnPdn
:
Pup
Pup
T
11.22 CV
Pup
Pdn
Pdn
:
Pup
Pup
crit
crit
QSCFM =
22.67 CV Pup 1
15.11 CV Pup
where
X=
X
3XT
X
:X XT
T
XT
T
Pup Pdn
Pdn
=1
Pup
Pup
:X XT
Pdn
Pup
+1/
Pdn
Pdn
Pup
Pup
crit
Pdn
Pdn
:
Pup
Pup
crit
12
m2RT
A2L x + VX,2
13
P1 =
P2 =
10
11
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
Orientation
Incoming flow
Payload direction of travel
Supply pressure
Representative test velocity
Velocities tested
Acceleration
Coulomb friction force
Flow rate source
Potential energy rate, E P = Mgx
Kinetic energy rate, EK = Mx x
Energy rate from Coulomb friction, EF = Ffricx
System studied
Rod end up
Blind end
Up
340 kPa 50 psi gauge
200 mm/s
201200 mm/s
20 m / s2
N/A
75 SLPM data
3.1 J/s
6.4 J/s
N/A
N/A
0.14 J/s
338 J/s
588 J/s
x =
1
m1RT
m2RT
Mg +
PatmArod bx Ffric
M
x + xT,1 L x + xT,2
14
where xT,1 and xT,2 represent the equivalent cylinder lengths of the
excess volumes in the cylinder and hose,
xT,1 =
VX,1
A1
15
xT,2 =
VX,2
A2
16
17
2RT
P2A2x = m
18
19
20
where the error term from atmospheric pressure acting over the
area of the rod, e, is defined as
041101-4 / Vol. 131, APRIL 2009
e=
PatmArodx
RT
21
22
Because the mass flow rates are equal in magnitude and opposing
in direction, the sum m1 + m2 must be constant,
m1 + m2 = m
23
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
Rather, they are associated with the transient spike in the air flow
at the beginning of each cycle, occurring concurrently with the
acceleration of the cylinder.
Table 5 lists the slope of the best-fit lines for the experiments
conducted, and Table 6 lists the R2 values. Figures 5 and 6 chart
these data. While agreement with the constant mass assumption is
good, the general trend is that the error increases as a function of
velocity. Equation 21 predicts this trend, as it contains an error
term proportional to velocity. Also contributing to the error is the
acceleration of the cylinder. As the command velocity increases, a
larger fraction of the duty cycle occurs with the cylinder under
acceleration, violating the preconditions set for the constant mass
assumption.
The constant mass assumption can remain valid even when the
cylinder undergoes stick-slip motion, characterized by short periods of rapid motion with longer periods in which friction prevents
motion. Stick-slip is a common phenomenon in compliant systems
with friction at low command velocities. In Fig. 7a the experimental system is observed to undergo stick-slip behavior on retraction; Fig. 7b shows the constant mass assumption remaining
Fig. 3 Simulated mass flow rates, from Ref. 21, 254 mm/s
command velocity, and 655 kPa supply. a Flow versus time
and b flow versus 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
Table 5 Slope of best-fit line, flow from rod end versus flow from blind end of cylinder
Commandvelocity
mm/s
170 kPa,
rod end up
170 kPa,
rod end down
340 kPa,
rod end up
340 kPa,
rod end down
520 kPa,
rod end up
520 kPa,
rod end down
20
50
100
200
300
400
800
1200
0.8857
0.9312
0.9397
0.9013
N/Aa
0.8733
0.8564
0.8561
0.8992
0.9506
0.9431
0.9108
N/Aa
0.8824
0.8770
0.8764
0.9686
0.9976
0.9791
0.9359
0.9076
0.8928
N/Ab
N/Ab
0.9996
1.0117
0.9880
0.9501
0.9233
0.8976
N/Ab
N/Ab
0.9452
1.0082
0.9942
0.9572
0.9338
0.9164
N/Ab
N/Ab
1.0178
1.0172
1.0063
0.9795
0.9443
0.9081
N/Ab
N/Ab
No data collected.
Flow sensors saturated at this combination of pressure and velocity.
Table 6 R2 values, flow from rod end versus flow from blind
end of cylinder
Command 170 kPa, 170 kPa, 340 kPa, 340 kPa, 520 kPa, 520 kPa,
velocity rod end rod end rod end rod end rod end rod end
down
up
down
up
down
up
mm/s
20
50
100
200
300
400
800
1200
a
0.9465
0.9289
0.9374
0.9215
N/Aa
0.9056
0.9088
0.8927
0.9468
0.9224
0.9201
0.9139
N/Aa
0.9013
0.8754
0.8539
0.9710
0.9674
0.9696
0.9546
0.9409
0.9315
N/Ab
N/Ab
0.9644
0.9678
0.9635
0.9486
0.9392
0.9300
N/Ab
N/Ab
0.9838
0.9808
0.9784
0.9689
0.9612
0.9412
N/Ab
N/Ab
0.9776
0.9784
0.9746
0.9679
0.9473
0.9334
N/Ab
N/Ab
xeq =
24
1m 2 m
2m 1
m
m1 + m22
25
1
m
m
26
1
Mg kcylx xeq bx PatmArod Ffric
M
27
P 1A 1
P 2A 2
+
xeq + xT,1 L xeq + xT,2
28
No data collected.
Flow sensors saturated at this combination of pressure and 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
without external loads or friction and does not reach either limit of
travel. The constant mass assumption was verified by both simulation and experimentation.
Table 7 State equations for pneumatic cylinder models
Conventional model x , x , P1 , P2
RT
P1A1
1
m
x
A1x + VX,1
A1x + VX,1
RT
P2A2
2+
m
x
P2 =
A2L x + VX,2
A2L x + VX,2
1
x = Mg + P1A1 P2A2 PatmArod bx Ffric
M
2 = fvalve
1,m
m
P1 =
Conclusions
3
4
711
Mass-based model x , x , m1 , m2
x =
1
m1RT
m2RT
Mg +
PatmArod bx
M
x + xT,1 L x + xT,2
14
Ffric
2 = fvalve
1,m
m
711
Reduced-order model x , x , xeq
k x x 1
cyl eq
+ Mg bx PatmArod Ffric
M
M
1
m
xeq = L + xT,1 + xT,2
m
1 = fvalve
m
x =
29
26
711
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
The authors wish to acknowledge the Nestl R&D Center, Inc.,
in Marysville, OH, for their support of this project.
Nomenclature
A
b
CD
c p, c V
CV
d
E
F
Fdyn
Fstat
k
L
M
m
area
coefficient of viscous damping
discharge coefficient
specific heats
flow coefficient
diameter
energy
force
Coulomb friction
static friction stiction
stiffness or gain
stroke length
combined mass of payload, piston, and rod
gas mass
P
QSCFM
R
T
V
X
x
Subscripts
1
2
air
atm
crit
cyl
dn
eq
fric
rod
S
spring
T
up
valve
X
pressure
volumetric flow rate, SCFM
gas constant
temperature
volume
pressure drop ratio
cylinder position
ratio of specific heats
density
Appendix A
The state equation for a gas, also known as the ideal gas law,
P = RT
A1
A2
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 PT PV
cVVT = m
dt
A3
d cV
c PT
PV + PV = m
dt R
A4
cP
cV
A6
P1A1
P1A1
P1A1 VX,1
xeq +
x +
x
x
x
A1
1
k1x
2
B8
3
2
B9
A7
k =
1
m1RT
3
x
2
B10
RT
P 1A 1
1
m
x
A1x + VX,1
A1x + VX,1
A8
A similar analysis may be performed on the rod end of the cylinder. In this case, the deflection is x as it acts in the opposite
direction on chamber 2 than on chamber 1,
RT
P 2A 2
2+
m
x
A2L x + VX,2
A2L x + VX,2
A9
k = F = P2A2
2
x
x
P1 =
P2 =
0 = P 1A 1 +
A5
RT
P
m
V
P =
V
V
where VX,1 and VX,2 are the excess volumes of gas in either end of
the cylinder and in the plumbing between the cylinder and the
valve.
B11
Appendix B
Assuming adiabatic compression and expansion of the gas in
the cylinder, the initial equilibrium energy is equal to the displaced internal energy state of the gas plus the potential energy
stored in the compressed gas,
Eeq = Eair + Espring
B1
0 = P 2A 2 +
+
Eeq,1 =
P1A1xeq + VX,1
1
Eair,1 =
P1 + P1A1xeq + x + VX,1
1
B4
The potential energy stored by the gas in the blind end, acting as
a spring, is given by
1
Espring,1 = 2k1x2
B5
where k1 is the effective nonlinear spring rate for the blind end,
defined as the change in force per unit displacement.
k = F = P1A1
1
x
x
Combining Eqs. B1 and B3B5,
Journal of Fluids Engineering
B6
B13
k =
2
m2RT
B3
3
P2A2 = k2 L xeq + xT,2 + x
2
B2
P2A2
P2A2
P2A2 VX,2
L xeq
x +
x
x
x
A2
1
k2x
2
B12
3
x
2
B14
B15
m1
3
x
2
m2
3
x
2
B16
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
k RT
cyl
m1
m2
+
xeq + xT,12 L xeq + xT,22
B17
P 1A 1
P 2A 2
+
xeq + xT,1 L xeq + xT,2
B18
References
1 Bobrow, J. E., and Jabbari, F., 1991, Adaptive Pneumatic Force Actuation and
Position Control, ASME J. Dyn. Syst., Meas., Control, 113, pp. 267272.
2 Richardson, R., Plummer, A. R., and Brown, M. D., 2001, Self-Tuning Control of a Low-Friction Pneumatic Actuator Under the Influence of Gravity,
IEEE Trans. Control Syst. Technol., 92, pp. 330334.
3 Moore, P. R., and Pu, J. S., 1996, Pneumatic Servo Actuator Technology,
IEE Colloq. on Actuator Technology: Current Practice and New Developments
110, pp. 3/13/6.
4 Shih, M., and Ma, M.-A., 1998, Position Control of a Pneumatic Cylinder
Using Fuzzy PWM Control Method, Mechatronics, 8, pp. 241253.
5 Richer, E., and Hurmuzlu, Y., 2000, A High Performance Pneumatic Force
Actuator System: Part II-Nonlinear Controller Design, ASME J. Dyn. Syst.,
Meas., Control, 122, pp. 426434.
6 McDonell, B. W., and Bobrow, J. E., 1993, Adaptive Tracking Control of an
Air Powered Robot Actuator, ASME J. Dyn. Syst., Meas., Control, 115, pp.
427433.
7 Gross, D. C., and Rattan, K. S., 1998, An Adaptive Multilayer Neural Network for Trajectory Control of a Pneumatic Cylinder, IEEE International
Conference on Systems, Man, and Cybernetics, San Diego, CA, Vol. 2, pp.
16621667.
8 Barth, E. J., Zhang, J., and Goldfarb, M., 2003, Control Design for Relative
Staility in a PWM-Controlled Pneumatic System, ASME J. Dyn. Syst.,
Meas., Control, 1253, pp. 504508.
9 Shih, M.-C., and Ma, M.-A., 1998, Position Control of a Pneumatic Cylinder
Using PWM Control Method, Mechatronics, 8, pp. 241253.
10 Shearer, J. L., 1956, Study of Pneumatic Processes in the Continuous Control
of Motion With Compressed AirI, Trans. ASME, 78, pp. 233242.
11 Shearer, J. L., 1956, Study of Pneumatic Processes in the Continuous Control
of Motion With Compressed AirII, Trans. ASME, 78, pp. 243249.
12 Liu, S., and Bobrow, J. E., 1988, An Analysis of a Pneumatic Servo System
and Its Application to a Computer-Controlled Robot, ASME J. Dyn. Syst.,
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
Zongchang Qu
Hua Yang
Bingfeng Yu
Institute of Compressor,
Xian Jiaotong University,
Xian 710049, China
The synchronal rotary compressor (SRC) has been developed to resolve high friction and
severe wear that usually occur in conventional rotary compressors due to the high relative velocity between the key tribo-pairs. In this study, the working principle and structural characteristics of the SRC are presented first. Then, the kinematic and force models
are established for the key componentscylinder, sliding vane, and rotor. The velocity,
acceleration, and force equations with shaft rotation angle are derived for each component. Based on the established models, numerical simulations are performed for a SRC
prototype. Moreover, experiments are conducted to verify the established models. The
simulated results show that the average relative velocity between the rotor and the cylinder of the present compressor decreases by 8082% compared with that of the conventional rotary compressors with the same size and operating parameters. Moreover, the
average relative velocity between the sliding contact tribo-pairs of the SRC decreases by
9394.3% compared with that of the conventional rotary compressors. In addition, the
simulated results show that the stresses on the sliding vane are greater than those on the
other components. The experimental results indicate that the wear of the side surface of
the sliding vane is more severe than that of the other components. Therefore, special
treatments are needed for the sliding vane in order to improve its reliability. These
findings confirm that the new SRC has lower frictional losses and higher mechanical
efficiency for its advanced structure and working principle. DOI: 10.1115/1.3089534
Keywords: synchronal rotary compressor, dynamic model, motion, forces
Introduction
Rotary type compressors are used more widely than reciprocating compressors in refrigeration and air-conditioning systems because they have advantages such as less components, simpler
structure, better dynamic equilibrium characteristics, and higher
reliability. However, owing to the high relative velocity between
the rotor, the cylinder, and the sliding vane of the conventional
rotary compressors, the frictional loss is high and then limits their
performances. Much work has been devoted to reducing the friction and wear of the conventional rotary compressors.
Jeon and Lee 1 discussed the tribological characteristics of
sliding surfaces between the vane and the flange in the rotary
compressor and conducted experiments to evaluate the effects of
different hard coatings on the compressor performance. Suh et al.
2 and Demas and Polycarpou 3 conducted compressor tribological experiments by using a high tribometer to explore the contact conditions. Oh et al. 4 and Oh and Kim 5 investigated the
friction and wear characteristics of sliding surfaces with various
coatings for a rotary compressor. They conducted experiments
with different working media and pressures in order to find out the
best surface treatment method to improve the tribological properties. Huang and Shiau 6 described a method to improve the
rotary vane compressor performance by employing extended rods
on both edges of each vane and guide slots on both cover plates.
In addition, Huang and Li 7 also established an optimum model
and found the optimal tolerance allocation for a vane rotary compressor. Cai et al. 8 proposed a perfect profile of the vane tip for
a rotary vane compressor to reduce the friction and wear of the
vane. Ooi 9 predicted that a 50% reduction in mechanical loss
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received October 24, 2007; final manuscript received January 9, 2009; published online March 6, 2009. Assoc. Editor:
Chunill Hah.
can be achieved by employing a multivariable, direct search, constrained optimization technique. Lee and Oh 10 conducted tribological experiments under various operation conditions and proposed an optimum initial surface roughness to prolong the wear
life of sliding surfaces.
The research mentioned above on reducing the friction and
wear is mostly focused on surface treatment and structural optimization for rotary compressors. However, the research cannot
reduce the friction and wear caused by high relative velocity between the key moving parts in the conventional rotary compressors radically. This paper develops an innovative synchronal rotary mechanism, in which the high friction and severe wear
caused by high relative velocity in the conventional rotary compressors are effectively reduced by the method of the cylinder and
the rotor rotating around their own axis synchronously.
Recently, we have performed the research on kinematics characteristics of the vane 11 and leakage model 12 of the SRC.
However, no studies have been published on the whole system
dynamic model of the SRC. The purpose of this study is to investigate the performance of the SRC by establishing the dynamic
model and testing the dynamic performance.
Figure 1 shows the working principle and structural characteristics of the proposed SRC. It can be seen from Fig. 1 that the
machines main components consist of a rotor, a cylinder, a sliding vane, a shaft, and two end covers.
The outside wall of the rotor and the inside wall of the cylinder
are tangent to each other, and they form the working chamber
together. The rotor is driven by the shaft, and it rotates around its
center OR at the angular velocity . The cylinder is driven by the
connecting sliding vane, and it rotates around its own center OC at
the angular velocity C. The sliding vane serves as a connector of
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 rotor and the cylinder, and it separates the working chamber
into a suction chamber and a compression chamber.
We define the starting point as the tangent point, as the shaft
rotation angle, and as the discharge angle. The typical stages of
the working process are shown in Fig. 1. It can be seen in Fig. 1
that when = 0, the discharge stage is finished. When 0 ,
the compression chamber volume decreases with the increase in
, and the gas is compressed; meanwhile, the suction chamber
volume increases and the gas is sucked in. During this stage, the
suction chamber volume equals the compression chamber volume
when = . When = , the discharge process starts and the gas
is expelled from the compression chamber. When = 2, the discharge is completed; the compression chamber closes and the suction chamber is full of fresh gas again. Thus, one operating cycle
is completed.
Because the relative velocity between the sliding contact tribopairs of the SRC is much smaller than that of the conventional
rotary compressors, the SRC has many advantages such as lower
friction and wear, easier to achieve seals, and lower vibration and
noise level.
e2 sin cos
dx
= e sin
dt
R + RS2 e2 sin2
4
According to the theorem of acceleration composition, the absolute acceleration aSG of point G equals the vector sum of the
convected acceleration aSGe, the relative acceleration aSGr, and the
Coriolis acceleration aSGk,
aSG = aSGe + aSGr + aSGk
where
aSGe = 2G = 2 e cos + R + RS2 e2 sin2 + RS L/2
aSGr = x = d2x/dt2
= 2e cos 2e2cos 2R + RS2 e2 sin2 0.5
+ 1/4e2 sin22R + RS2 e2 sin2 1.5
aSGk = 2 vSGr
3.2 Velocity and Acceleration of the Cylinder. Unlike the
conventional rotary compressors, the cylinder of the SRC is not
fixed but driven by the sliding vane and rotates synchronously
with the rotor around its own center. In order to analyze the kinematics characteristics of the cylinder, a rotational coordinate system OCXY, which rotates synchronously with the cylinder, is
established.
Figure 3 shows the velocity triangle of the cylinder. Because
the hinge joint hole center of the cylinder coincides with the head
center OS of the sliding vane, the absolute velocity vC for the
hinge joint hole center of the cylinder equals the absolute velocity
vS for the point OS on the sliding vane. On the basis of the above
motion model of the sliding vane, the absolute velocity vC at shaft
rotation angle can be expressed as
vC = vS = vSe + vSr
where
Fig. 2 Velocity triangle of the sliding vane
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
vSr = x =
e2 sin cos
dx
= e sin
dt
R + RS2 e2 sin2
Accordingly, the cylinder angular velocity C and angular acceleration C can be expressed as
3.3 Kinematic Simulation and Analysis. Based on the kinematic models established above, simulations of the relative velocity between the key moving components can be performed. Figure
4 shows the comparison of average relative velocity between the
rotor and the cylinder of the SRC and the conventional rotary
compressors with the same design and operating parameters. In
Fig. 4, it can be seen that the average relative velocity between the
rotor and the cylinder of the SRC decreases by approximately
8082% compared with that of the conventional rotary compressors. As a result, the disadvantages caused by high relative velocity between the rotor and the cylinder such as severe wear, difficulty to control the meshing clearance between the rotor and the
cylinder, and difficulty to achieve seals can be overcome.
For a SRC and a rolling piston compressor RPC, the vane and
the rotor are in sliding contact and the large frictional loss exists
between them. For a rotary vane compressor RVC, the vane tips
and the cylinder are in sliding contact and the large frictional loss
exists between the vane tips and the cylinder. The frictional loss
would decrease if the relative velocity between these sliding contact tribo-pairs could decrease effectively. Figure 5 shows the
comparison of average relative velocity between the sliding contact tribo-pairs of the SRC and the conventional rotary compressors with the same design and operating parameters. It can be seen
in Fig. 5 that the average relative velocity between the sliding
contact tribo-pairs is 0.8 m/s for the SRC, 11.3 m/s for the RPC,
and 13.6 m/s for the RVC. Thus, the average relative velocity
between the sliding contact tribo-pairs of the SRC decreases by
92.894.1% compared with that of the conventional rotary compressors. Therefore, the disadvantages caused by high relative velocity between the sliding contact tribo-pairs such as large frictional loss and severe wear could be overcome completely.
4.1.1.1 Gas force FCg and moment M Cg. The working chamber is separated into a suction chamber and a compression chamber by the sliding vane, and the gas force FCg acting on the cylAPRIL 2009, Vol. 131 / 041102-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
LCSn = OC P = e sin
LCS = OS P = e cos + R + RS2 e2 sin2 + e cos
14
15
16
= arccose2 + 2 R + RS2/
Pd Ps
2
10
11
4.1.1.3 Forces from lubricating oil at the meshing point. Lubricating oil exists in the meshing clearance of the rotor and the
cylinder. Because of the relative motion between the rotor and the
cylinder, there are normal fluid hydrodynamic forces FCMn and
tangential frictional forces FCM from the lubricating oil acting on
the cylinder. FCMn is transmitted to the shaft and causes moment
of flexion, which can be ignored because it is much smaller than
that of the other moments on the shaft.
According to the lubrication theory, the tangential frictional
force FCM and frictional moment M CM at the meshing point can
be obtained by
FCM =
or2H C
o
or2RH C
M CM =
o
18
4.1.1.4 Frictional force from the head of the sliding vane. The
sliding vane rotates in the hinge joint hole of the cylinder and
causes the sliding frictional force FCSf , which is given by
FCSf = f CS FCS = f CS FCSn2 + FCS2
19
20
1/2
The gas force crosses the center OC of the cylinder, and thus the
moment caused by gas at the center OC is zero.
4.1.1.2 Driving force FCS and moment M CS from the sliding
vane. The driving force on the cylinder by the sliding vane can be
divided into two components, namely, the normal component FCSn
and the tangential component FCS, which is normal and tangential to the cylinder surface, respectively. The arms of FCSn and
FCS at the cylinder center OC are given as
041102-4 / Vol. 131, APRIL 2009
17
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
end part by the vane slot, and FSRf2 is the corresponding frictional
force.
Both FSR1 and FSR2 are normal to the vane side surface. The
moments at the head center OS of the sliding vane caused by FSR1
and FSR2 are given by
M SR1 = FSR1 LSR1 = FSR1 e cos + R + RS2 e2 sin2 r
27
M SR2 = FSR2 LSR2 = FSR2 L RS
FSRf1 and FSRf2 act along the contact surface of the sliding vane
and the vane slot, and their directions are opposite to the velocity
vSr. FSRf1 and FSRf2 are given by
FSRf1 = f SR FSR1
FSRf2 = f SR FSR2
22
Accordingly, the frictional force and moment acting on the cylinder are given by
FCBf = f CBFCB
M CBf = FCBf RB
23
4.1.2 Differential Equation of Rotation of the Cylinder. According to the differential equations of rotation of a rigid body
with a fixed axis, the differential equation of rotation of the cylinder can be expressed as
J C C =
M F = M
i
CS
+ M CM + M CBf
24
4.2 Force Model for the Sliding Vane. The forces acting on
the sliding vane are the gas force, the contact forces and frictional
forces from the vane slot, the constrained forces from the cylinder,
and the inertia forces due to the vane motion, etc. Figure 7 shows
all the forces acting on the sliding vane. The forces are assumed to
be positive when they deviate from the rotor, and the moments are
assumed to be positive when they inhibit the sliding vane rotation.
4.2.1.1 Gas force FSg . The pressure difference between the
suction chamber and the compression chamber exerts gas force to
the sliding vane. The gas force is given by
FSg = Pd Ps l H
= Pd Ps e cos + R2 e2 sin2 r H
25
Thus, the moment M Sg at the head center OS of the sliding vane
by the gas force can be given as
26
4.2.1.2 Constrained force FSC by the cylinder. The constrained force FSC on the sliding vane by the cylinder can be
divided into the normal component FSCn and the tangential component FSC. FSCn and FSC are equal and opposite to the driving
forces FCS and FCSn on the cylinder by the sliding vane, respectively.
4.2.1.3 Contact forces FSR1, FSR2 and frictional forces FSRf1,
FSRf2 from the rotor. As shown in Fig. 7, FSR1 is the contact force
caused by the vane slot up margin, and FSRf1 is the corresponding
frictional force. FSR2 is the contact force acting on the sliding vane
Journal of Fluids Engineering
29
28
FSk = mC aSGk
30
FSr = mC aSGr
The corresponding moments caused by the inertia forces may be
expressed as
M Se = 0
M Sk = FSk LSk = FSk L/2 RS
31
M Sr = 0
A set of simultaneous equations including the normal and tangential force balance equations of the sliding vane, the moment
balance equation for the head of the sliding vane, and the differential equation of the cylinder rotation are obtained as
FCSn + FSe + FSr FSR1 f SR FSR2 f SR = 0
FCS + FSa + FSk + FSR1 FSR2 = 0
FSgLSg + FSkLSk + FSR1 LSR1
FSR2 LSR2 FSR1 f SR
32
b
b
FSR2 f SR = 0
2
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
2C r4 rin4
Or rin
M RCf =
f SR
f SR
LCS LCSn
0
FSg FSk
FSe FSr
J C C
FSgLSg FSkLSk
FCSn
FSR2
33
or2HC
O
or3HC
M RM =
O
FSR1
4.3.1.1 Gas force FRg. The gas force acting on the rotor is
caused by the pressure difference between the suction chamber
and the compression chamber. The gas force FRg crosses the rotor
center OR and is given by
34
The moment M Rg created by the gas force is zero because the gas
force crosses the center of the rotor.
4.3.1.2 Reaction forces FRS1, FRS2 and frictional forces FRSf1,
FRSf2 from the sliding vane. These four forces acting on the rotor
from the sliding vane are force couples with the corresponding
FSR1, FSR2, FSRf1, and FSRf2 described previously in the force
model for the sliding vane.
The moments at the rotor center OR caused by these four forces
are given by
37
4.3 Force Model for the Rotor. In the SRC, the rotor rotates
at a uniform angular velocity . The forces and moments acting
on the rotor are shown in Fig. 8, including the gas force, the
reaction forces and frictional forces on the vane slot from the
sliding vane, the frictional forces from the two end covers, the
normal dynamic flow force and the tangential frictional force created by the lubricating oil at the meshing point M of the rotor and
the cylinder, and the driving moment from the shaft.
FRM =
FCS
All the four unknown forces can be obtained by solving the above
matrix. Based on these four forces, all the other unknown forces
acting on the cylinder and the sliding vane could be obtained.
36
Likely, the frictional force and moment acting on the rotor by the
oil at the meshing point can be calculated by
C r4 rin4
38
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 connector that drives the cylinder to rotate and at the same time
is driven by the rotor. In Fig. 10, it can be seen that the principal
forces acting on the sliding vane are the contact forces by the vane
slot, the tangential constrained forces by the cylinder, and the gas
force. The other forces are relatively much smaller, varying over a
range from 38 N to 79 N. Moreover, in Fig. 10, the tangential
constrained force by the cylinder fluctuates from 247.9 N to
237.9 N and changes from negative to positive at the rotation
angle , which not only increases the machine vibration but also
aggravates the fatigue damage of the sliding vane. It is suggested
that the inertia of the cylinder should be reduced in order to mitigate the negative impacts on the sliding vane.
5.3 Simulation Results for Forces on the Rotor. Figure 11
shows the simulation results of the forces acting on the rotor.
From Fig. 11, it can be noted that the gas force is the greatest in
all forces, and its maximum value occurs at the shaft rotation
angle of 220 deg, where the discharge process just starts. All the
other forces are less than 346 N. In addition, it can be seen that the
resistant moment consists mainly of the acting forces by the sliding vane, and the frictional moments are relatively smaller.
Model Verification
Model Verification.
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
Operating parameters
Evaporating temperature k
Suction pressure MPa
Discharge pressure MPa
Subcooling K
Superheat K
Refrigerant
0.120
0.108
0.070
0.008
0.029
0.006
ment within the rotation angle range of 190270 deg, and the
former is greater than the latter within the other rotation angle
range. The error between the measured torque input and the simulated resistant moment reaches the maximum value 11 N.m at
the shaft rotation angle of 220 deg, which may be caused by the
buffering of the rotor and the shaft coupling.
Conclusions
The dynamic model has been established, and numerical simulations have been conducted for a novel SRC in order to study its
dynamic characteristics. Motion and stress for its components
have been discussed in detail. From the research, the following
conclusions could be drawn:
a
c
Fig. 13 Wear of the sliding vane
273
0.3
1.5
7
10
R134a
Acknowledgment
The authors acknowledge the National Natural Science Foundation of China NSFC for funding this research project as part
of the Design for Fluid Mechanism Program Grant No.
50476053.
Nomenclature
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
r rotor radius
rin rotor inside radius
R
RB
RS
t
vC
vS, vSe, vSr,
V
V d
x
C
C
cylinder radius
cylinder bearing radius
sliding vane radius
time
absolute velocity of the cylinder
absolute, convected, and relative velocity for
point OS of the vane
absolute, convected, and relative velocity
for point G of the vane
total volume of the working chamber
compression chamber volume at rotation angle
displacement
angle, = ORBOC
bow angle of the oil film
shaft rotation angle
included angle between of the FCS with axis
X
angle, = MOCB
length of OSOR
length of ORB
length of ORG
rotor angular velocity
cylinder angular velocity
cylinder angular acceleration
kinematic viscosity of the lubricating oil
oil film thickness
References
1 Jeon, H-.G., and Lee, K-.S., 2007, Friction and Wear of Lubricated Sliding
Surfaces of Coated Vane and Flange in Rotary Compressor With Carbon Dioxide as a Refrigerant, Asian Pacific Conference for Fracture and Strength,
pp. 17851788.
2 Suh, A. Y., Patel, J. J., and Polycarpou, A. A., 2006, Scuffing of Cast Iron and
Al390-T6 Materials Used in Compressor Applications, Wear, 260, pp. 735
744.
3 Demas, N. G., and Polycarpou, A. A., 2006, Tribological Investigation of
Cast Iron Air-Conditioning Compressor Surfaces in CO2 Refrigerant, Tribol.
Lett., 22, pp. 271278.
4 Oh, S.-D., Lee, K.-S., and Kim, J.-S., 2006, Experimental Investigation on
Friction and Wear of Coated Vane Under the Environments of Lubricants and
CO2 as a Refrigerant, Key Eng. Mater., 326328II, pp. 11851188.
5 Oh, S. D., and Kim, J. W., 2004, Friction and Wear Characteristics of TiN
Coated Vane for the Rotary Compressor in a R410a Refrigerant, Tribol.
Trans., 481, pp. 2933.
6 Huang, Y. M., and Shiau, C.-S., 2006, Optimal Tolerance Allocation for a
Sliding Vane Compressor, ASME J. Mech. Des., 128, pp. 98107.
7 Huang, Y. M., and Li, C. L., 2005, The Stress of Sliding Vanes in a Rotary
Compressor, 2005 World Tribology Congress III, pp. 99100.
8 Cai, H., Li, L., and Guo, B., 2005, Research on Tip Profile of Vane for Rotary
Vane Compressor, Fluid Machinery Group-International Conference on Compressors and Their Systems, pp. 215222.
9 Ooi, K. T., 2005, Design Optimization of a Rolling Piston Compressor for
Refrigerators, Appl. Therm. Eng., 25, pp. 813829.
10 Lee, Y.-Ze., and Oh, S.-D., 2003, Friction and Wear of the Rotary Compressor VaneRoller Surfaces for Several Sliding Conditions, Wear, 255, pp.
11681173.
11 Zhou, H., Qu, Z., and Feng, J., 2005, Sliding Vane Kinetic Characteristics
Analysis for a SRC, J. Shanghai Jiaotong Univ., 399, pp. 982984.
12 Zhou, H., Qu, Z., and Yu, B., 2007, Leakage Research in Synchronal Rotary
Gas Compressor, China J. Mech. Eng., 182, pp. 205208.
13 Solzak, T. A., and Polycarpou, A. A., 2006, Tribology of WC/C Coatings for
Use in Oil-Less Piston-Type Compressors, Surf. Coat. Technol., 201, pp.
42604265.
14 Xu, H., Cai, C., Yan, J., Wang, K., and Zhou, S., 2000, Machine Design
Handbook 2, China Machine Press, Beijing, China, Chap. 16, pp. 1623.
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
Angelo Cervone
Postdoctoral Fellow
e-mail: angcervone@mbox.me.es.osaka-u.ac.jp
Yoshinobu Tsujimoto
Professor
e-mail: tujimoto@me.es.osaka-u.ac.jp
Engineering Science,
Osaka University,
1-3 Machikaneyama, Toyonaka,
Osaka 560-8531, Japan
The paper will present an analytical model for the evaluation of the pressure and flow
rate oscillations in a given axial inducer test facility. The proposed reduced order model
is based on several simplifying assumptions and takes into account the facility design and
the dynamic properties of the tested inducer. The model has been used for evaluating the
dynamic performance of a prototype of the LE-7 engine liquid oxygen (LOX) inducer, in
tests carried out under given external flow rate excitations. The main results of these
calculations will be shown, including the expected oscillations under a wide range of
operational conditions and the influence of facility design. Calculations showed that the
only way to obtain the two linearly independent test conditions, necessary for evaluating
the inducer transfer matrix, is by changing the facility suction line: Any other changes in
the facility design would result ineffective. Some other important design indications provided by the analytical model will be presented in the paper. DOI: 10.1115/1.3089535
Yutaka Kawata
Professor
Osaka Institute of Technology,
5-6-1 Omiya, Asahi-Ku,
Osaka 535-8505, Japan
e-mail: kawata@med.oit.ac.jp
Introduction
Space rocket turbopumps are one of the most crucial components of all primary propulsion concepts powered by liquid propellant engines. Generally these pumps include an axial inducer,
upstream of the centrifugal stages, in order to improve the suction
performance and reduce the propellant tank pressure and weight.
Severe limitations are associated with the design of high power
density, dynamically stable machines capable of meeting the extremely demanding suction, pumping, and reliability requirements
of space transportation systems 1.
It is well known that many of the flow instabilities acting on
axial inducers are significantly influenced by the so-called dynamic matrix, which relates the pressure and flow oscillations at
the pump inlet with those at the pump outlet 24. The first steps
in the analytical and experimental characterizations of the dynamic matrix date back to the work of Brennen and co-worker in
the 1970s 5,6. However, more recent works have given important contributions by evaluating the previously obtained results
through a careful analysis of the successive experimental and numerical data see, for example, Refs. 7,8.
Conventionally, the dynamics of hydraulic systems is treated in
terms of lumped-parameter models, which assume that the distributed physical effects between two measuring stations can be
represented by lumped constants. This assumption is usually considered valid when the geometrical dimensions of the system are
significantly shorter than the acoustic wavelength at the considered frequency. As a direct consequence of this assumption, the
dynamic matrix of a generic system like an axial inducer can be
written as
pd
Q
d
H11 H12
H21 H22
pu
Q
u
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
Experimental Apparatus
assumption is not valid in real pumps under cavitating conditions 7,8,11, it could represent a good starting point for
the simplified reduced order analysis presented in this paper. As it will be shown in the following, the results and the
main conclusions drawn by the analyses are not significantly affected by this assumption, especially for what concerns the design of the experiment.
Under the above assumptions, pressure and flow rate can be
written in a complex form, as functions of time, as follows:
pt = p + p eit
e it
+Q
Qt = Q
usually real are the pressure and flow rate
where p and Q
usually complex are the pressure and
steady values, p and Q
flow rate oscillating components, and is the frequency of the
oscillations.
The relevant components for the evaluation of the dynamic matrix are the oscillating ones. Under the above assumptions, the
matrices of all the components of the facility can be easily evaluated as follows.
Model Description
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.1 Pipes. With a suitable manipulation of the generic continuity and momentum equations, the dynamic matrix can be written as follows:
H=
1 R iL
1
Q
u
A u2
L=
A d2
A u2
Q
u
iC 1
V
p
1
1
rT d
+ Q
u
2
Au d
Ad
A u2
1 S + iX R iL
iC
1 i M
+ CLind2 i
rT
R + iL =
Ad
i M
A u2
Q
rT
u
2
Ad
Ad
C=
2r T2
Q
rT
u
C
Ad
A d2
Q
u
VC
+ V C c
in which p and V are the mean pressure and volume of the air at
the top of the tank, and is the specific heat ratio of air.
R=
1
VC
AurT
V C , = a 1 c 4 + a 2 c 3 + a 3 c 2 + a 4 c
M=
A d2
where the constants a1, a2, a3, and a4 can be evaluated using
the four experimental values of C provided by Brennen 13.
In the above equation, c is the choked cavitation number of the
BrennenAcosta partially cavitating cascade model i.e., the cavitation number for which the cavity becomes infinitely long 13,
which depends on the inducer geometry and the incidence angle
and, as a consequence, the flow coefficient.
The cavity volume at the choked cavitation number, VCc, is
evaluated assuming the following.
1 At the choked cavitation number c, complete breakdown
occurs i.e., the head rise of the inducer falls to zero.
2 At c, the blade channels are completely filled by bubbly
cavitation, having a void fraction 0, which can be evaluated using the relationship shown by Brennen 5,
=
p2
+ MLind
02 02
21 02 sin2 b
Q
2
HM 11 HM 12
HM 21 HM 22
p1
Q
1
+ iLind
HM =
1 Rout iLout
0
1
1 Rin iLin
0
1 S + iX Rind iLind
iCind
1 i M
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
Considering the components of the facility between the downstream measurement section and the flow rate fluctuator see Fig.
1, it is possible to write the oscillations immediately after the
fluctuator as
+Q
=Q
Q
3
2
F
0.255
for
= 1 Hz
0.057
for
= 2 Hz
p3 = p2 RA + iLAQ
2
0.023
for
= 3 Hz
0.007
for
= 5 Hz
0.002
for
= 10 Hz
p1
Q
1
HB11 HB12
HB21 HB22
p3
Q
3
1 RC iLC
1
iCT 1
1 RB iLB
0
RB + iLBQ
F
=Q
=
Q
1
2
Rm + RA + RB + RC + iLm + LA + LB + LC
p1 =
RC + iLCRB + iLBQ
Q
F
F
+
iCT Rm + RA + RB + RC + iLm + LA + LB + LC
p2 =
+ R + iL
+ L R + iL Q
R
Q
C
C
B
B
F
F
+
iCT Rm + RA + RB + RC + iLm + LA + LB + LC
p2a
Q
2a
p2b
Q
2b
HM 11 HM 12
HM 21 HM 22
HM 11 HM 12
HM 21 HM 22
p1a
Q
1a
p1b
Q
1b
where
Rm = Rin + Rind + Rout
Lm = Lin + Lind + Lout
The first term in the expressions of p1 and p2 is usually very
small due to the large value of the tank compliance CT and can
therefore be neglected, at least for sufficiently high values of the
exciting frequency . As an example, for the Osaka University
facility with LE-7 engine LOX inducer under nominal operating
conditions = 0.078, = 3000 rpm, mean tank pressure p
det T =
2
RBa + iLBaRBb + iLBbRCb RCa + iLCb LCaQ
F
Rm + RAa + RBa + RCa + iLm + LAa + LBa + LCaRm + RAb + RBb + RCb + iLm + LAb + LBb + LCb
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
4.2 Evaluation of Alternative Facility Designs. The effectiveness of different facility designs has been evaluated by means
of the backward calculation. The results provided by the model
have therefore been used as a starting point for the evaluation of
the dynamic matrix elements, in order to compare this evaluation
to the actual nominal values provided by the model. The dynamic
matrix elements have been evaluated at different excitation frequencies and flow coefficients .
In order to better simulate the real experimental conditions expected in the facility, three different artificial error sources have
been intentionally added to the nominal data. The first two error
sources are a random white noise at all the frequencies included
in the spectrum of the measured signal which simulates the typical white noise usually present in laboratory environments and a
defined sinusoidal noise at 60 Hz frequency of electrical disturbance in Japan. The third error source consists in taking only a
finite number of decimals from the numerical solutions provided
by the model, in order to simulate the experimental sensors accuracy.
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 Amplitude and phase of the upstream pressure oscillations, as a function of the excitation frequency cavitating
conditions, for = 3000 rpm, = 0.1, = 0.078, and amplitude
of the flow rate oscillations generated by the fluctuator equal to
1 l/s
Fig. 8 Amplitude and phase of the downstream pressure oscillations, as a function of the excitation frequency cavitating conditions, for = 3000 rpm, = 0.1, = 0.078, and amplitude of the flow rate oscillations generated by the fluctuator
equal to 1 l/s
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 Amplitude and phase of the upstream flow rate oscillations, as a function of the excitation frequency cavitating
conditions, for = 3000 rpm, = 0.1, = 0.078, and amplitude
of the flow rate oscillations generated by the fluctuator equal to
1 l/s
Fig. 10 Amplitude and phase of the downstream flow rate oscillations, as a function of the excitation frequency cavitating conditions, for = 3000 rpm, = 0.1, = 0.078, and amplitude of the flow rate oscillations generated by the fluctuator
equal to 1 l/s
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
A.C. would like to acknowledge the support of Japan Society
for the Promotion of Science JSPS, which has funded his research program at Osaka University. He also expresses his gratitude to all the colleagues at Pisa University in Italy, and, in particular, Professor Luca dAgostino and Professor Mariano
Andrenucci, for their kind and constant encouragement.
Nomenclature
ai constants in cavity volume equation,
i = 1 4 m3
A pipe section m2
C compliance m4 s2 / kg
k loss coefficient dimensionless
L inertance kg/ m4
M mass flow gain factor s1
p pressure Pa
H generic dynamic matrix
HB dynamic matrix outside the measurement
section
HM dynamic matrix of the measurement section
Q flow rate m3 / s
QF externally imposed flow oscillation m3 / s
R resistance kg/ m4 s
rT inducer tip radius m
S real part of the pressure gain factor
dimensionless
t time s
T characteristic matrix for backward calculation
V volume m3
VC cavity volume m3
X imaginary part of the pressure gain factor s
0 void fraction dimensionless
b inducer tip blade angle rad
head loss at breakdown dimensionless
flow coefficient dimensionless
specific heat ratio dimensionless
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
inertial length m
head coefficient dimensionless
density kg/ m3
cavitation number dimensionless
choked cavitation number dimensionless
frequency of oscillations Hz
inducer rotating speed rad/s
Superscripts
oscillating component of a given quantity
steady component of a given quantity
Subscripts
a , b linearly independent experimental conditions
A pipework from measurement section to
fluctuator
B discharge line
C suction line
d downstream section generic
in measurement section upstream of inducer
ind inducer
m total measurement section including inducer
out measurement section downstream of inducer
T main tank
u upstream section generic
1 upstream measurement point
2 downstream measurement point
3 point immediately after the fluctuator
References
1 Stripling, L. B., and Acosta, A. J., 1962, Cavitation in TurbopumpsPart 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
Computational fluid dynamics (CFD) was used to evaluate the flow field and thrust
performance of a promising concept for reducing the noise at take-off of dual-stream
turbofan nozzles. The concept, offset stream technology, reduces the jet noise observed on
the ground by diverting (offsetting) a portion of the fan flow below the core flow, thickening and lengthening this layer between the high-velocity core flow and the ground
observers. In this study a wedge placed in the internal fan stream is used as the diverter.
Wind, a Reynolds averaged NavierStokes (RANS) code, was used to analyze the flow
field of the exhaust plume and to calculate nozzle performance. Results showed that the
wedge diverts all of the fan flow to the lower side of the nozzle, and the turbulent kinetic
energy on the observer side of the nozzle is reduced. This reduction in turbulent kinetic
energy should correspond to a reduction in noise. However, because all of the fan flow is
diverted, the upper portion of the core flow is exposed to the freestream, and the turbulent
kinetic energy on the upper side of the nozzle is increased, creating an unintended noise
source. The blockage due to the wedge reduces the fan mass flow proportional to its
blockage, and the overall thrust is consequently reduced. The CFD predictions are in
very good agreement with experimental flow field data, demonstrating that RANS CFD
can accurately predict the velocity and turbulent kinetic energy fields. While this initial
design of a large scale wedge nozzle did not meet noise reduction or thrust goals, this
study identified areas for improvement and demonstrated that RANS CFD can be used to
improve the concept. DOI: 10.1115/1.3089536
Introduction
test points required, additional small scale testing, and large scale
testing 9.63 in. 24.5 cm fan exit diameter in the NASA Glenn
Nozzle Acoustic Test Rig NATR 4.
There are four other papers from the OST study that form a
complete description of the program. A companion CFD study
examining the vane deflector and s-duct nozzles was performed
by Dippold et al. 5. The experimental results from NATR, acoustics and particle image velocimetry PIVbased velocity measurements in the exhaust plume, were reported by Brown et al.
6, and a comparison of experimental results between UCI and
NASA Glenn was detailed by Zaman et al. 7. The results of the
DOE study were reported by Henderson et al. 8.
The experimental results at NATR indicate that the wedge deflectors do reduce the noise on the lower side of the nozzle at
downstream angles. However, at other angles on the lower side
and at all angles on the sideline and upper side, the noise was
increased. In general, the overall estimated perceived noise level
EPNL increased relative to the baseline due to the higher noise
at the other locations. One configuration, BPR= 5.0, M = 0.0,
showed a small noise reduction of 0.8 EPN dB. Based on measurements of turbulent kinetic energy, Brown et al. 6 concluded
that the wedge designs tested were too aggressive.
Thrust performance is a critical factor in assessing low noise
nozzle concepts. Papamoschou 13 estimated the thrust loss in
nozzles with vane deflectors based on airfoil drag and nozzle flow
deflection 2, and Papamoschou et al. 9 experimentally examined the aerodynamics of a wedge deflector in a simplified configuration. But up to this point, there has been no detailed examination of the performance of these nozzles. The NATR tests did
not include force measurements. A key element of the CFD calculations presented here was to characterize the thrust of these
nozzles and compare them to a baseline nozzle without deflectors.
This work represents the highest fidelity estimates of Offset
Stream nozzle performance to date.
This paper describes the CFD analyses performed to support
the design and testing of turbofan nozzles with the wedge shaped
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
core flow nozzle section was replaced to affect the bypass ratio
change. The 4BB nozzle has an internal plug in the core flow and
has a bypass ratio of 8. This nozzle was studied computationally
but not tested in the experiment. All three nozzles have the same
fan diameter of 9.63 in. 24.5 cm.
The shape of the wedge deflector studied here is the same as the
wedge used by Papamoschou 13 in his initial studies. The
wedge has a half angle of 11 deg and a base width of 3.44 in.
8.74 cm. The upper and lower surfaces of the wedge conform to
the internal contours of the fan stream. The wedge reduces the exit
area of the fan stream by 14%.
It is anticipated that the wedge deflector would be incorporated
as a pair of flaps into the existing structure of the engine. For
underwing mounted engines found on most large transport aircraft, the flaps could be placed in the pylon. For other installations, side mounted or more highly integrated installations, the
flaps could be incorporated into a strut. In either case the flaps
would deploy to form a wedge on take-off and retract for cruise.
The pylons for underwing mounted engines already deflect some
fan flow. However, additional flow deflection is probably necessary to maximize noise reduction.
The calculations performed for the pretest predictions and performance analyses were done using cycle data for a representative
take-off condition. Details of the operating conditions are in Table
1. The experiment was run subsequent to the CFD analyses, and
the data were taken at slightly different conditions. The analyses
were rerun at these conditions for code validation purposes Table
2. For the 3BB validation case, the nozzle was tested with static
freestream conditions. This analysis was run with a freestream
Mach number of 0.05 for numerical stability. Two 5BB cases were
run, one with heated streams and one with cold streams. Both
cases had a freestream Mach number of 0.20. Because the wedge
deflectors will most likely be removed from the flow after takeoff, only that flight condition was analyzed.
Analysis Procedure
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
Core
Freestream
BPR
NPR
T0
K
U
m/s
NPR
T0
K
U
m/s
p
kPa
T
K
5.0
8.0
8.0
1.83
1.57
1.62
364
347
356
340
290
303
1.68
1.52
1.42
833
844
832
480
437
399
98.6
98.6
98.6
294
294
294
0.29
0.29
0.29
Core
Freestream
BPR
NPR
T0
K
U
m/s
NPR
T0
K
U
m/s
p
kPa
T
K
5.0
8.0
8.0
1.83
1.62
1.42
356
356
298
336
303
239
1.68
1.42
1.86
833
832
298
480
399
311
98.6
98.6
98.6
298
298
298
0.00a
0.20
0.20
ing the following: 1 the internal primary flow path, 2 the external primary flow path including the plug, 3 the internal secondary flow path, 4 the external secondary flow path, 5 the
nacelle, and 6 the jet plume region Fig. 4. Dimensions of the
grids are given in Table 3. At block interfaces the grids are not
contiguous, and a Roe-based zone coupling algorithm is used to
pass the fluxes between blocks.
The grids for the baseline nozzles with no wedge deflector were
planar, taking advantage of the axisymmetric geometry. A total of
66,734 points were used. The grids for the nozzles with the wedge
deflector were three dimensional: a 180 degree segment, taking
advantage of the symmetry about the xy-plane vertical streamwise plane. A total of 6,072,794 grid points were used.
3.2 Flow Solution. The flow solver used for all the computations was Wind version 5.0 12. Wind is a general purpose Reynolds averaged NavierStokes RANS equation solver. All the
calculations presented here used the same code options. The default numerical scheme was used: implicit time stepping and a
second-order physical Roe scheme for stretched grids. Turbulence
Baseline
Axial
Radial
Azimuthal
Total
73
53
65
45
69
253
45
165
45
109
45
173
91
91
91
91
91
91
298,935
795,795
266,175
446,355
282,555
3,982,979
6,072,794
66,734
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
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
contours in Figs. 510. Data are presented for all three nozzles
3BB, 4BB, and 5BB, and the results show that all three flow
fields behave similarly.
Velocity contours in the streamwise planes compare the wedge
nozzle to the baseline nozzle Figs. 57. These contours show
that the presence of the wedge directs the core flow upward and
significantly shortens the core streams potential core. The layer of
fan flow on the underside of the jet plume is substantially thicker.
Fig. 8 Streamwise xy-plane contours of turbulent kinetic energy, k / Uc2, for 3BB nozzles
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 Streamwise xy-plane contours of turbulent kinetic energy, k / Uc2, for 4BB nozzles
stream, closer to the nozzle exit. On the lower side of the exhaust
plume, the levels of k / U2c have been reduced. This is the desired
effect. A small local peak in turbulence levels is seen directly
behind the wedge, again indicating a separated region and potential noise source.
The movement of the fan flow induced by the wedge deflector
is more easily seen in cross-planes downstream of the nozzle exit
Figs. 1116. The fan flow moves to the lower half of the exhaust
Fig. 10 Streamwise xy-plane contours of turbulent kinetic energy, k / Uc2, for 5BB nozzles
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
plume, and at x / D f = 5.0 the fan stream is almost completely below the core flow. The upper portion of the core flow comes in
direct contact with the freestream flow and is not buffered by the
fan flow. This region of large velocity gradients results in high
turbulence levels and will certainly increase the noise on this side
of the nozzle. At x / D f = 3.0 high turbulence levels are present
from the top of the nozzle down to the sideline position 090
deg, measured from the vertical. The turbulent kinetic energy
levels resulting from the core stream interacting directly with the
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 Cross-plane yz-plane contours of turbulent kinetic energy, k / Uc2, for 3BB nozzles
Nozzle Performance
4.2.1 Calculation Procedure. Nozzle performance was computed using the Wind utility, CFPOST. Stream thrust was com-
puted at the fan and core exits. Pressure and viscous forces were
computed on the remaining surfaces: core nozzle lip, fan nozzle
lip, plug, splitter between fan and core, nacelle, and wedge base.
Nozzle performance will be reported in two ways. First, the
total thrust of the nozzle with wedge deflector will be expressed as
a fraction of the baseline nozzle thrust
F/Fbaseline
Fig. 15 Cross-plane yz-plane contours of turbulent kinetic energy, k / Uc2, for 4BB nozzles
Fig. 16 Cross-plane yz-plane contours of turbulent kinetic energy, k / Uc2, for 5BB nozzles
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
Second, it will be expressed in terms of a thrust coefficient. Typically nozzle performance results are reported using the gross
thrust coefficient,
C fg =
Factual
actualUideal
m
Using the actual mass flow in the denominator, C fg, removes the
deviation from the ideal mass flow from the performance measure.
Discharge coefficient is a separate measure for the reduction in
mass flow,
Cd =
actual
m
ideal
m
For subsonic nozzles the mass flow through the nozzle can vary
due to downstream effects, and the thrust coefficient does not give
a true measure of nozzle performance. In these circumstances, the
change in the measured thrust is caused by the change in the
nozzle flow rate and not a change in the exit velocity. This is
especially true for two-stream nozzles where the behavior of one
stream can affect the flow rate of the other. Nozzles with noise
suppressing features can also affect flow rates. Therefore, a thrust
coefficient that accounts for the change in mass flow rate is used
in this paper. It is simply a thrust coefficient based on the ideal
mass flow rather than the actual mass flow. It can be conveniently
written as the product of the gross thrust coefficient and discharge
coefficient,
Journal of Fluids Engineering
actual Factual
Factual
m
=
= CdC fg
idealUideal m
ideal m
actualUideal
m
For the cases where the wedge is installed in the nozzle, the ideal
mass flow is computed using the reduced nozzle exit area.
A baseline case, no wedge, was calculated for each nozzle geometry. Performance of the wedge nozzles is referenced to the
baseline and is reported as CdC fg,
CdC fg = CdC fgwedge CdC fgbaseline
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
exit area and therefore does not include the loss due to blockage.
The external plug nozzles have uniformly higher performance
than the internal plug nozzle. Figure 18b shows the decrement in
performance due to the wedge. The losses due to the wedge vary
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
from about 0.7% to 1.1% and appear to be relatively low, considering the large base area created by the wedge. A detailed breakout of the forces that comprise the thrust calculation is shown for
the 5BB baseline and wedge nozzles in Fig. 19. The figure confirms that the primary loss mechanism is the loss in fan stream
thrust due to the reduced mass flow. The presence of the wedge
reduces the positive force on the nozzle plug and increases the
drag on the nacelle. The base drag on the wedge surface is a
relatively minor component. This is somewhat surprising but can
be explained by the streamlines in Fig. 20. The figure shows that
the flow that separates and recirculates behind the base of the
wedge is the external flow over the nacelle. This also shows the
linkage between the presence of the wedge and the increased nacelle drag Fig. 19. The fan flow, the primary contributor to the
thrust, is turned by the wedge to create the offset stream but does
not separate behind the wedge, which would result in large thrust
losses.
By examining nozzle performance using the CdC fg parameter,
we can see that the losses created by the wedge deflector would be
acceptable if the reduction in mass flow could be removed or
accounted for in the nozzle design. The wedge deflector could be
Fig. 22 Comparison to experiment; u-velocity, u / Uc, for the 3BB nozzle with wedge, hot conditions
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. 23 Comparison to experiment; turbulent kinetic energy, k / Uc2, for the 3BB nozzle with wedge, hot
conditions
a viable noise reduction concept if the wedge size could be significantly reduced or if an external wedge that did not reduce the
fan flow rate was employed.
The wedge induces a downward movement of the fan flow and
a corresponding upward movement of the core flow. The resulting
vertical component of the nozzle thrust is presented in Fig. 21.
The data indicate that the wedge induces a very small vertical
component to the thrust. The external plug nozzles have a smaller
component, about 0.5% of the axial thrust, and the internal plug
has a slightly larger component, approximately 1.5% of the axial
thrust.
The wedge deflector was one of three devices used to offset the
fan stream examined in NASAs OST program. The performance
of the other two devices, vane deflectors and s-duct nozzles, were
also analyzed using a similar CFD method 5. The performance
of all three devices are compared in the Appendix of Ref. 17.
4.3 Comparison to Experiment. Once the experimental test
program was carried out, the CFD analyses were rerun to match
041104-12 / Vol. 131, APRIL 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
Fig. 24 Comparison to experiment; u-velocity, u / Uc, for the 5BB nozzle with wedge, hot conditions
compared with the predictions for three wedge nozzle cases listed
in Table 2. Comparisons for the baseline nozzle cases can be
found in Ref. 17.
The 3BB nozzle was tested at conditions similar to the pretest
predictions, but with a static freestream. The CFD analysis used a
small freestream velocity, M = 0.05, to prevent numerical instabilities. Plots showing the comparison between prediction and experiment for both axial velocity and turbulent kinetic energy are
shown in Figs. 22 and 23. Near the fan exit plane at x / D f = 1.0 and
3.0, the proximity of the nozzle plug to the measurement planes
caused a low signal to noise ratio in the experimental data. These
data were removed, resulting in the white regions in the center of
the plots. The white areas at the corners of the experimental data
are due to the orientation of the nozzle relative to the PIV measurement plane. Overall, the agreement between CFD and experiment is very good. The contour plots indicate that the CFD does a
good job of predicting the shape of the nozzle plume and the
movement of diverted fan flow relative to the core. One does see
that the CFD predicts a more diffuse plume shape compared with
experiment. Line plots along the y-axis vertical axis show that
Journal of Fluids Engineering
the CFD does a good job of predicting the levels of both velocity
and turbulent kinetic energy. For this case the Wind code predicts
higher peak velocities. This is due to the inclusion of the forward
velocity in the freestream for numerical reasons. This small coflow serves to reduce the nozzle mixing rate 19.
The 5BB nozzle was tested at two different flow conditions.
The quality of the PIV data near the nozzle exit was much better
with only small regions at x / D f = 1.0 where the data were removed. The first case had the same nozzle conditions as the pretest prediction, but the freestream Mach number was 0.20 instead
of 0.29. The comparison is presented in Figs. 24 and 25. For this
case agreement is very good for both velocity and turbulent kinetic energy. The contour plots show similar agreement in plume
shape. The line plots show that there is excellent agreement in the
prediction of both velocity and turbulence levels. The improved
agreement in the levels over the 3BB case is most likely due to the
fact that the CFD analysis matched the experimental freestream
velocity. The second 5BB case was tested at cold nozzle temperatures and a Mach 0.20 freestream. The comparison for this case is
presented in Figs. 26 and 27. Again, for all aspects agreement
APRIL 2009, Vol. 131 / 041104-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
Fig. 25 Comparison to experiment; turbulent kinetic energy, k / Uc2, for the 5BB nozzle with wedge, hot
conditions
Detailed Reynolds averaged NavierStokes calculations of turbofan nozzles with wedge deflectors for noise reduction were carried out to support NASAs Offset Stream Technology program.
The nozzles were designed to reduce jet noise at take-off condi041104-14 / Vol. 131, APRIL 2009
tions. The computational fluid dynamics calculations were intended to characterize the exhaust plume aiding in wedge design,
predict nozzle thrust performance, and validate the computational
technique against experimental data.
The predicted flow fields showed that the wedge deflector was
effective at redirecting the fan flow to the underside of the exhaust
plume, effectively reducing noise on this side of the nozzle. The
predictions indicate that the wedge design is too aggressive, deflecting too much flow and exposing the core flow to the
freestream on the top and sides of the plume. Predictions of turbulent kinetic energy suggest that noise levels at the sideline and
upper side of the nozzle may approach those of an isolated jet at
core flow conditions.
Nozzle performance calculations show that the thrust of the
nozzle is significantly reduced relative to an unmodified nozzle.
This thrust loss is caused by a reduction in mass flow of the fan
stream, which is directly proportional to the area reduction of the
fan exit caused by the wedge. Thrust losses imparted by the
wedge not due to the fan flow reduction are relatively small. If 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
Fig. 26 Comparison to experiment; u-velocity, u / Uc, for the 5BB nozzle with wedge, cold conditions
Nomenclature
k
m
p
u
x,y,z
BPR
Cd
C fg
D
F
M
NPR
T
U
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. 27 Comparison to experiment; turbulent kinetic energy, k / Uc2, for the 5BB nozzle with wedge, cold
conditions
Subscripts
actual
baseline
ideal
c
f
0
measured quantity
nozzle without flow deflectors
isentropically expanded to freestream condition
core primary
fan secondary
stagnation/total condition
freestream condition
References
1 Papamoschou, D., 2002, Noise Suppression in Moderate-Speed Multistream
Jets, AIAA Paper No. 2002-2557.
2 Papamoschou, D., 2004, New Method for Jet Noise Reduction in Turbofan
Engines, AIAA J., 4211, pp. 22452253.
3 Papamoschou, D., and Nishi, K., 2005, Jet Noise Suppression With Fan Flow
Deflectors in Realistic-Shaped Nozzle, AIAA Paper No. 2005-993.
4 Castner, R., 1994, The Nozzle Acoustic Test Rig; An Acoustic and Aerodynamic Free-Jet Facility, NASA Report No. TM 106495.
5 Dippold, V., Foster, L., and Wiese, M., 2007, Computational Analyses of
Offset Stream Nozzles for Noise Reduction, AIAA Paper No. 2007-3589.
6 Brown, C., Bridges, J., and Henderson, B., 2007, Offset Stream Technology
TestSummary of Results, AIAA Paper No. 2007-3664.
7 Zaman, K., Bridges, J., and Papamoschou, D., 2007, Offset Stream
TechnologyComparison of Results From UC and GRC Experiments, AIAA
Paper No. 2007-438.
8 Henderson, B., Norum, T., and Bridges, J., 2006, An MDOE Assessment of
Nozzle Vanes for High Bypass Ratio Jet Noise Reduction, AIAA Paper No.
2006-2543.
9 Papamoschou, D., Vu, A., and Johnson, A., 2006, Aerodynamics of WedgeShaped Deflectors for Jet Noise Reduction, AIAA Paper No. 2006-3655.
10 Janardan, B. A., Hoff, G. E., Barter, J. W., Martens, S., Gliebe, P. R., Mengle,
V., and Dalton, W. N., 2000, AST Critical Propulsion and Noise Reduction
Technologies for Future Commercial Subsonic Engines Separate-Flow Exhaust System Noise Reduction Concept Evaluation, NASA Report No. CR
2000-210039.
11 Pointwise, Inc., http://www.pointwise.com/gridgen/
12 The NPARC Alliance, http://www.grc.nasa.gov/www/winddocs/
13 Menter, F. R., 1994, Two-Equation Eddy-Viscosity Turbulence Models for
Engineering Applications, AIAA J., 328, pp. 15981605.
14 Papamoschou, D., 2003, A New Method for Jet Noise Reduction in Turbofan
Engines, AIAA Paper No. 2003-1059.
15 Khavaran, A., and Kenzakowski, D. C., 2007, Noise Generation in Hot Jets,
NASA Report No. CR 2007-214924.
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
16 Papamoschou, D., 2006, Fan Flow Deflection in Simulated Turbofan Exhaust, AIAA J., 4412, pp. 30883097.
17 DeBonis, J., 2008, RANS Analyses of Turbofan Nozzles With Wedge Deflectors for Noise Reduction, AIAA Paper No. 2008-0041.
18 Wernet, M., and Bridges, J., 2002, Application of DPIV to Enhanced Mixing
Heated Nozzle Flows, AIAA Paper No. 2002-0691.
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
e-mail: t.vanholten@tudelft.nl
Monique Heiligers
Annemie Jaeken
Department of Design, Integration and Operation
of Aircraft and Rotorcraft,
Faculty of Aerospace Engineering,
Delft University of Technology,
Kluyverweg 1,
2629 HS Delft, The Netherlands
The behavior of a vortex flow through a Laval nozzle was studied in connection with the
purification of natural gas. By creating a vortex and passing it through a Laval nozzle,
the gas will be cooled, and water droplets will form and will be centrifuged out of the
gas. This system is named the Condi-Cyclone. An analytical theory is developed to reveal
the most important phenomena of the flow, to first order accuracy. Experiments have been
performed with a prototype of the Condi-Cyclone. A Euler numerical simulation was
performed, using the geometry of the test channel. This paper presents an analytical
theory for a vortex flow through a Laval nozzle. It will demonstrate that when a vortex is
present the total velocity reaches sonic conditions upstream of the nozzle throat, that the
axial component of the velocity in the nozzle throat is equal to the local speed of sound
and that the mass flow through the Laval nozzle decreases with increasing vortex
strength. The predictions of the analytical theory have been compared with the results of
the experiments and the Euler numerical simulation, and it can be concluded that the
analytical theory describes the main characteristics of the flow very well.
DOI: 10.1115/1.3089532
Keywords: Laval nozzle, vortex flow, natural gas, Condi-Cyclone
Introduction
The behavior of a vortex flow through a Laval nozzle was studied in connection with the purification of natural gas. Water normally present in natural gas needs to be extracted because it can
have damaging effects on the installations and pipelines further
downstream. The steps of dehydration and hydrocarbon removal
are normally achieved in multiple process stages, for example, in
turbo-expanders or with added chemicals, notably glycol or
methanol 1. In addition to requiring large vessels and complex
operation, processes such as glycol regeneration can also produce
toxic atmospheric emissions, including benzene, toluene, ethylene, and xylene 1.
To solve these problems a radical new means to purify natural
gas has been invented: the Condi-Cyclone. By creating a vortex in
a Laval nozzle, the gas will be cooled, and water droplets will
form and will be centrifuged out of the gas. This paper will
present an analytical theory describing the basic process of a vortex flow through a Laval nozzle.
Section 2 will provide a brief explanation of the Condi-Cyclone
and its lay-out. Section 3 will present the actual analytical theory.
To verify the predictions of the analytical theory, both Euler numerical simulation and full scale tests with the Condi-Cyclone
have been performed. The results of the full scale tests will be
presented in Sec. 4; the results of the numerical simulation in Sec.
5. The different results will be compared with the analytical
theory in Sec. 6, and finally the conclusions are presented in Sec.
7.
Condi-Cyclone Lay-Out
condensate, and small water droplets will form. These water droplets will be centrifuged to the wall of the Condi-Cyclone by the
vortex and can then be extracted from the flow. A schematic layout of the Condi-Cyclone is given in Fig. 1.
After the natural gas enters the inlet of the Condi-Cyclone, a
static center-body directs the gas flow past a static array of blades.
This static array of blades deflects the flow and thereby generates
the vortex. This part of the Condi-Cyclone is also referred to as
the subsonic vortex generator SVG. The next step is an acceleration of the flow by means of a Laval nozzle. The flow will
become supersonic, and when the temperature drop is sufficient,
water droplets will start to form. When the water droplets have
grown sufficiently, they are driven toward the nozzle wall by the
centrifugal forces and are extracted from the flow in the recovery
section. When the flow has passed the recovery section, it is
slowed down. The cleaned gas then flows straight into the remaining processes.
Analytical Theory
The purpose of the analytical theory is to reveal the most important phenomena of the flow, to first order accuracy. In the case
of rocket nozzles and the design of supersonic windtunnels, it is
well known that the assumption of quasi-one-dimensional isentropic flow may give a useful first indication of the global flow
properties. For this reason it was attempted to develop a theory for
the Condi-Cyclone flow based on similar assumptions. The validity of such an extremely simplified theory will depend on how
strong viscous effects are or on how much the flow is affected by
local shocks. Comparisons with experiments Sec. 4 are therefore
essential before conclusions can be drawn concerning the validity
of the assumptions. Further insight into details such as the occurrence of local shocks will be provided by Euler numerical simulations Sec. 5. The analytical theory for a vortex flow through a
Laval nozzle will be presented in this section. Among others it
will be shown that the total velocity reaches the local speed of
sound upstream of the throat and that the mass flow is a function
of the vortex strength.
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 2
Tt
M
=1+
2
T
1 2
pt
M
= 1+
2
p
1 2
t
M
= 1+
2
1 2
M
a = RT = RTt 1 +
2
2
at2
M2 =
at2
1 2
M
1+
2
V2
=
a2
1/1
/1
A=
m
Vax
1 2
V
2
Vtan
+ V2ax
1 2
Vtan + V2ax
at2
2
1 dVax
1 d M 2 dVax M 2 dVtan
m
dA
+
=
+
dx
Vax dM 2 Vax dx
Vtan dx
Vax dx
10
1 d M 2
1
+
2
dM Vax
Vax
11
1 2
1
d
M
= t 1 +
2
2
dM 2
/1
1 2
1
M
= 1+
2
2
1 a2
= 2
2 at
12
M2
=
Vax
2Vaxat2
1 2
Vtan + V2ax
at2
2
= 2Vax
at2
a4
13
= VaxA = constant
m
Vax = athroat
Vtanr = constant
I=m
The total Mach number of the flow in the throat section, taking
into account that there is an additional velocity component Vtan,
apparently is larger than unity.
at2 a2 = RT
1 2 1 2
Tt
M =
V
1 = a2
2
2
T
so that
041201-2 / Vol. 131, APRIL 2009
M thr 1
2
=
M thr
if Vtan 0
2
2
Vtan
Vtan
+ V2ax
1 2
M thr
=1+ 2 1+
2
2
a
at
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
1+
1 2
M bl
2
+1/21
20
at sin blrblM bl 1 +
M blabl sin blrbl = m
I=m
1 2
M bl
2
1+
2
M thr
=
2
Vtan
at2
2
1 Vtan
2 at2
16
The latter equation again clearly shows that the Mach number in
the throat is larger than unity if a vortex flow is present. Under
such conditions the transition from sub-to supersonic flow occurs
at a position upstream of the throat.
3.3 Mass Flow Through the Channel. Making use of the
property 14 and Eq. 5, the mass flow through the throat section
and therefore the mass flow in every section may be expressed
as follows:
= VaxA = thrathrAthr
m
17
= tatAthr 1 +
m
1 2
M thr
2
+1/21
18
If, for example, is taken equal to the value for air = 1.4, the
exponent in Eq. 18 will equal to 3. It can thus be seen that
increasing M thr will lead to a reduced mass flow. In other words,
increasing the vortex strength will throttle the Laval nozzle, as a
consequence of the additional Vtan and the resulting increased
Mach number in the throat. This throttling effect is quite large, as
will be seen later from some experimental results used to determine the validity of the above given theory.
Using the average radius rthr of the flow annulus in the throat
see the earlier remarks about the quasi-one-dimensional approximation underlying the developed theory, Eq. 6 yields for the
conditions in the throat:
Vtan,thr =
I
rthr
m
1/2
21
, I, and M bl,
Equations 20 and 21 contain the three variables m
apart from the constants that specify the reservoir- or total- conditions of the gas and the geometry of the channel. The set Eqs.
20 and 21 therefore provides a second relation between mass
flow and angular momentum.
At this point there are two separate relations available between
and I; the first one derived from the throat conditions and the
m
second one given by the inflow conditions. The calculation of the
and I can thus be completed.
values of m
3.5 Computational Scheme to Determine Mass Flow and
19
Vx = Mxax =
Mxat
1+
1 2
M x
2
2
x + V2axx
Vx = Vtan
22
23
Mx2at2
2
x
Vtan
1 2
M x
1+
2
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
Analytical results
Experimental results
Input
Output
bl
deg
pt
bar
Tt
K
m
kg/s
m2 / s
Tend nozzle
K
M end nozzle
20
25
30
40
1.5
1.5
1.5
1.5
300
300
300
300
1.07
0.98
0.89
0.73
25.10
28.91
32.27
38.14
180
175
165
155
1.9
1.9
2.0
2.2
t 1 +
1 2
M x
2
1/1
M 2xat2
m
2
x =
Vtan
1 2
Ax
M x
1+
2
25
Solving Eq. 25 will give the Mach number for each axial coordinate: substitution of the Mach number for each axial coordinate
into Eq. 22 will give the total velocity for each axial coordinate.
Finally substitution of the Mach number into Eq. 24 will give
the axial velocity for each axial coordinate.
3.7 Analytical Calculations and Results. Using the analytical theory calculations have been performed on the CondiCyclone test setup that was used during the experiments see Sec.
4. Section 4 will also describe the geometry of the CondiCyclone. Here, the analytical results for four different blade
angles bl are already given see Table 1. The input parameters
are the total pressure and total temperature at the location of the
blades of the subsonic vortex generator. The output parameters are
the mass flow, vortex strength , temperature at the end of the
nozzle Tend nozzle, and Mach number at the end of the nozzle
M end nozzle.
Experiments
SVG
deg
pt
bar
Tt
K
2%
m
kg/s
20%
m2 / s
20
25
30
40
1.5
1.5
1.5
1.5
300
300
300
300
1.04
0.95
0.887
No result
26
24
31
No result
Fig. 4 General geometry of the channel with axisymmetric mesh. The zero station is located at the end of the center-body;
the throat is located at the 0.5 m station.
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
channel.
A crude check on the magnitude of the swirl was made by
mounting a three-hole pressure tube in the widest part of the channel near the aft end of the central body. This pitot tube could be
rotated in such a way that it was roughly aligned with the flow
direction. From this device the flow angle was derived on the
mean radius of this particular station.
The measured vortex strength was expressed in terms of the
circulation defined as
= 2rVtan = 2rblM blabl sin bl
26
The expected accuracy is not more than 20% due to the uncertainties in the measurement techniques. On account of the quasi-onedimensional assumption underlying the analytical theory, the circulation is related to the angular momentum I by the expression
I=
27
4.2 Experimental Results. Experiments with the CondiCyclone have been performed at the Stork Product Engineering
TM620-Air/Water test facility 2,3. The results of the test for
different blade angle deflections are given in Table 2. With increasing blade deflection angle and thus with increasing vortex
strength, the mass flow through the nozzle decreases. At the highest blade deflection angle 40 deg the test setup started to vibrate
heavily, and the mass flow was severely limited. Because of the
vibrations the test was aborted and therefore no results are available for this blade deflection angle. Either vortex breakdown or
the stalling of the blades of the SVG is thought to be the probable
cause for these vibrations.
Fig. 6 Tangential velocity horizontal axis as a function of radial distance vertical axis
for several stations
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
Output
SVG
deg
M bl
pt
bar
Tt
K
m
kg/s
m2 / s
Tend nozzle
K
M end nozzle
20
25
30
40
0.14
0.13
0.13
0.12
1.5
1.5
1.5
1.5
300
300
300
300
1.12
1.02
0.91
0.69
23.31
27.37
32.68
39.07
186
180
174
168
1.75
1.80
1.85
1.95
in the region of the tail of the center-body. In the nozzle throat the
axial velocity is equal to the local speed of sound, as predicted by
the analytical theory. In the nozzle a vortex structure is found with
a solid body rotation at the axis of symmetry and a 1 over r
behavior toward the nozzle wall see Fig. 6. A severe pressure
drop is present in the region of the solid body rotation.
With increasing blade deflection angle the mass flow decreases,
the vortex strength increases, a higher Mach number at the end of
the nozzle is achieved, and a slightly larger temperature drop occurs see Table 3. With increasing blade deflection angles the
solid body rotation region expands, the supersonic region at the
end of the center-body moves upstream, and the effective cross
section of the channel decreases. At the highest vane angle of 40
deg a flow disturbance at the very end of the center-body was
clear see Fig. 7, showing up as a region of back-flow near the
center-line and possibly involving a small shock. It was therefore
recommended to use in later configurations center-bodies with a
pluglike tail shape.
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 Mass flow and vortex strength as a function of the blade deflection angle
Fig. 9 Temperature and Mach number at the end of the nozzle as a function of the blade deflection angle
Mach number becomes supersonic well before the throat, just upstream of the end of the center-body. Both results also show that
the axial velocity in the throat is equal to the local speed of sound.
The Mach number and temperature at the end of the nozzle station 1.3 m downstream of the center-body as predicted by the
analytical theory and calculated by the numerical simulation are
plotted in Fig. 9.
Conclusions
Epilogue
array of blades, the concept performed very well during the experiments. The Condi-Cyclone concept has been patented by
Stork Product Engineering and has been further developed by
Twister BV to make the technology commercially available to the
wider industry. A Twister demonstration unit has been operational
in Nigeria since early 2001, successfully performing water and
hydrocarbon separation, while confirming the ability of Twister to
run virtually without supporting systems Twister BV. Earlier
units have been in operation in the Netherlands in Zuiderveen en
Barendrecht since 1997 5.
References
1 Knott, T., 2000, Twist in the Tale, Offshore Engineer, http://
www.oilonline.com/news/features/oe/20000701.Twist_in.67.asp.
2 Mertens, S., 1999, Test Report: Subsonic Vortex Generation 20 & 40 Degrees Turbine Ring, Stork Product Engineering, Amsterdam, The Netherlands, Technical Report No. SPE-CON-RP-102 1/0.
3 Mertens, S., 1999, Test report: Subsonic Vortex Generation 25 & 30 Degrees
Turbine Ring, Stork Product Engineering, Amsterdam, The Netherlands,
Technical Report No. SPE-CON-RP-113 1/0.
4 Schoones, M. M. J., and Van Holten, Th., 1999, Euler Numerical Simulation
of Channel Flow With a Superimposed Vortex, Faculty of Aerospace Engineering, Delft University of Technology, The Netherlands, Technical Report
No. FM&P-99.001.
5 Twister, B. V., http://www.twisterbv.com/
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
Mohsen Akbari1
Ph.D. Candidate
Mechatronic System Engineering,
School of Engineering Science,
Simon Fraser University,
Surrey, BC, V3T 0A3, Canada
e-mail: maa59@sfu.ca
David Sinton
Associate Professor
Department of Mechanical Engineering,
University of Victoria,
Victoria, BC, V8W 2Y2, Canada
Majid Bahrami
Assistant Professor
Mem. ASME
Mechatronic System Engineering,
School of Engineering Science,
Simon Fraser University,
Surrey, BC, V3T 0A3, Canada
Introduction
Advances in microfabrication make it possible to build microchannels with micrometer dimensions. Micro- and minichannels
show potential and have been incorporated in a wide variety of
unique, compact, and efficient cooling applications in microelectronic devices 1. These microheat exchangers or heat sinks feature extremely high heat transfer surface area per unit volume
ratios, high heat transfer coefficients, and low thermal resistances.
In biological and life sciences, microchannels are used widely for
analyzing biological materials such as proteins, DNA, cells, embryos, and chemical reagents 2. Various microsystems, such as
microheat sinks, microbiochips, microreactors, and micronozzles
have been developed in recent years 36. Since microchannels
are usually integrated into these microsystems, it is important to
determine the characteristics of the fluid flow in microchannels for
better design of various microflow devices.
In parallel to the recent advances in microfluidic devices, microfabrication techniques have evolved. Some of the important
fabrication techniques include lithography soft and photolithography, lamination, injection molding, hot embossing, laser micromachining, and electrochemical or ultrasonic technologies
2,710. Together with new methods of fabrication, it is possible
to exploit certain fundamental differences between the physical
properties of fluids moving in large channels and those traveling
through micrometer-scale channels 2. In recent years, a large
number of papers have reported pressure drop data for laminar
fully developed flow of liquids in microchannels with various
cross sections 1132. Rectangular cross sections have been extensively studied as they are employed in many applications
1116. However, published results are often inconsistent. Some
authors reported a huge deviation from the conventional theories
and attributed it to an early onset of laminar to turbulent flow
transition or surface phenomena such as surface roughness, electrokinetic forces, viscous heating effects, and microcirculation
near the wall 1226.
Jiang et al. 28 conducted an experimental investigation of
1
Corresponding author.
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received June 11, 2008; final manuscript
received December 5, 2008; published online March 6, 2009. Assoc. Editor: James
A. Liburdy. Paper presented at the ECI International Conference on Heat Transfer
and Fluid Flow in Microscale III, Whistler, BC, Canada, September 2126, 2008.
water flow through different microchannel cross sections including circular, rectangular, trapezoidal, and triangular. The hydraulic
diameter of microchannels varied from 8 m to 42 m. They
collected experimental data with the Reynolds number ranging
from 0.1 to 2, and concluded that there was less influence of the
cross-sectional shape on the microflow in the microchannel and
that the experimental data agreed well with the prediction of the
conventional theory. Baviere et al. 29 performed an experimental study on the water flow through smooth rectangular microchannels. Their channels were made of a silicon engraved substrate anodically bonded to a Pyrex cover. Their results showed
that in smooth microchannels, the friction law is correctly predicted by conventional theories. Judy et al. 15 and Bucci et al.
30 showed that their experimental results were in good agreement with conventional theories in the laminar regime. Also Wu
and Cheng 31, Liu and Garimella 32, and Gao et al. 33
reported good agreement between experimental data with conventional theories.
Low Reynolds numbers characterize many microscale liquid
flows 34. Hence, nonlinear terms in the NavierStokes equation
disappear for fully developed flow, resulting in Poissons equation
2u =
1 dp
,
dz
u = 0 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
Theory
L
A2
421 + 2
31 +
42
f ReA = 13.16;
0
=1
Experimental Setup
www.usa.autodesk.com.
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
Channel
Width, 2W
m
Depth, 2H
m
Length, L
mm
=H/W
LD,max / L
%
PPR-0.13
PPR-0.17
PPR-0.4
PPR-0.6
PPR-0.76
780
581
480
189
134
110
101
192
113
103
50
50
58.8
55.5
50
0.13
0.17
0.4
0.6
0.76
3.47
4.03
5.55
6.75
7.86
Viscous dissipation effect was neglected since the dimensionless number 4 EcL f ReA / ReA = 0.006 is much smaller than
1. The value of Tref = 1 C was used to calculate the Eckert number, Ec, for water as a working fluid 26. Hence, the properties of
the water was assumed to be constant.
Total measured pressure drop during the experiment, Pmeasured
is
Pmeasured = Pc + Pin + PD + PFD + Pex + 2Pb + Pev
6
where Pc is the pressure loss due to the flow in the connecting
tubes, Pin and Pex are the inlet and exit losses, PD is the
developing region loss, PFD is the pressure drop in the fully
developed region, Pb is the pressure drop due to 90 deg bends,
and Pev is the pressure drop corresponds to the electroviscous
effect. Since fully developed pressure drop is the focus of this
study, right hand side losses except PFD must be subtracted from
the measured pressure drop.
4.1 Connecting Tube Pressure Loss, Pc. The connecting
tube pressure drop includes the losses due to all fittings and the
capillary tube from the transducer to the microchannel inlet. We
measured this loss directly at each flow rate when there was no
microchannel at the end of the tubing. The measurements were
carefully conducted with the conditions identical to the case when
a microchannel was added to the end of the connecting tube to
avoid the effects of hydrostatic pressure.
4.2 Developing Region, PD. Since the viscous boundary
layer inherently grows faster in microchannels than in macroscales, the developing region in most cases is negligible. There are
few references that can be found in literature, in which the effect
of inlet region was considered 40,44,45. Phillips 44 showed
that the length of the hydrodynamic developing region, LD, depends on the aspect ratio of rectangular cross-section microchannels; the higher the aspect ratio, the longer the developing length.
Maximum value of LD can be obtained from Eq. 7
LD =
4 Q
1 + 2
f ReAQ
2AA
LD + K
Q2
2A2
www.ni.com.
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.13
= 0.17
0.01
0.01
0.01
0.02
0.02
0.03
0.01
0.02
0.02
0.01
0.02
0.05
= 0.4
0.04
0.03
0.04
0.06
0.06
0.10
= 0.6
= 0.76
0.03
0.04
0.06
0.07
0.08
0.20
0.04
0.06
0.08
0.10
0.11
0.24
2.99595.
In Eq. 8, f ReA was calculated based on the measured pressure drop. Table 2 lists the relative pressure loss of the developing
region with respect to the measured pressure drop. As can be seen,
the values are small, less than 0.3%, and can be neglected for the
range of Reynolds numbers studied in this work 1 ReA 35.
For the higher Reynolds numbers Re 100, the developing pressure drop, PD, was found to be less than 2% of the fully developed pressure drop obtained from Eq. 3.
4.3 Minor Losses, Pmin. Other pressure losses associated
with the measured pressure drop are inlet, exit, and bend losses.
These losses are usually obtained from the traditional relationships used in macroscale 15,44,46,47. Phillips 44 showed that
the minor pressure losses can be obtained from
Pmin = Pin + Pex + 2Pb =
Q2
A
Kc + Ke + 2Kb
2A
At
9
where A and At are the channel and connecting tube crosssectional areas, respectively. Kb is the loss coefficient for the
bend, and Kc and Ke represent the contraction and expansion loss
coefficients due to area changes. Phillips 44 recommended Kb to
be approximately 1.2 for a 90 deg bend. Assuming equal crosssectional areas for the channel and connecting tubes and also
maximum possible values for Kc and Ke 48, relative minor
losses with respect to the measured pressure drop are listed in
Table 3. As can be seen, these losses are negligible compared with
the measured pressure drop. For higher Reynolds numbers Re
100, minor pressure drop, Pmin, was found to be less than 5%
of the fully developed pressure drop obtained from Eq. 3.
4.4 Electroviscous Effect, Pev. When a liquid is forced
through a narrow channel under an applied pressure gradient, the
counterions in the diffusive layer of EDL electric double layer
are moving toward the downstream end and a potential gradient is
induced in the flow 49. This so called streaming potential acts to
bD2 eo
= 0.13
= 0.17
= 0.4
= 0.6
= 0.76
40
60
80
100
120
240
0.04
0.05
0.07
0.08
0.09
0.16
0.04
0.07
0.10
0.07
0.10
0.24
0.12
0.10
0.14
0.21
0.19
0.35
0.09
0.13
0.17
0.22
0.25
0.51
0.11
0.17
0.23
0.28
0.35
0.70
10
Poexp =
4 A2A
PFD
Q L
11
Poexp
Poexp
P
P
L
L
25 A
4 A
Q
Q
2 1/2
12
where values indicate uncertainty associated with each variable
2
2 1/2
+ H
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
Another comparison of the present experimental work with analytical model is illustrated in Fig. 7. Since the Poiseuille number,
f ReA, remains constant for the laminar regime as the Reynolds
number varies, the experimental data for each set were averaged
over the laminar region. The 10% bounds of the model are also
shown in the plot to better demonstrate the agreement between the
data and the model.
As shown in Eq. 3, the Poiseuille number, f ReA, is only a
function of the aspect ratio, which is a geometrical parameter.
This dependency is plotted in Fig. 8. Averaged values over the
studied range of the Reynolds number were used in this plot. It
can be observed that for smaller aspect ratios, the Poiseuille number, f ReA, increases sharply. To better show the trend of the
analytical model, two asymptotes of very narrow and square cross
sections are also included in Fig. 8. As can be seen, when 0,
parallel plate, the Poiseuille number asymptotically approaches
42 / 3.
Uncertainty
0.25% of full scale
0.02 mm
3.6 m
3.6 m
0.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
ducted, and the pressure drops due to the developing region, minor losses, and electroviscous effect were estimated using available models in literature. Channel deformation due to the pressure
inside the microchannel was found to be negligible since the pressure drop changes linearly with the flow rate. Comparing the results with the general theoretical model shows good agreement
between the model and experimental data.
Uncertainty analysis showed that the measurement of channel
dimensions and flow rate is critical in microscales. Here, microchannel cross-section geometry was determined through processing high quality images of the channel cross section.
According to the analytical model and experimental data of
present work, the Poiseuille number, f ReA, was found to be only
a function of microchannel geometry in the range of Reynolds
numbers studied in this work. For rectangular cross sections, the
channel aspect ratio is the only parameter that affects the Poi-
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
W
W
H
H
2 1/2
A4
Acknowledgment
The authors gratefully acknowledge the financial support of the
Natural Sciences and Engineering Research Council of Canada
NSERC. Also they thank Dr. Viatcheslav Berejnov and Mr. Ali
K. Oskooi for their assistance.
Nomenclature
Greek
A
At
cp
Ec
f
H
Ip
Ip
Kb
Kc
Ke
L
LD
L
Po
Q
ReA
W
Pb
Pc
PD
Pex
PFD
Pev
Pin
a
eo
R =
R
i
xi
2 1/2
A1
R
=
R
ai
i
xi
2 1/2
A2
A =
A
W
W
A
H
H
2 1/2
A3
References
1 Bahrami, M., Yovanovich, M. M., and Culham, J. R., 2006, Pressure Drop of
Laminar, Fully Developed Flow in Microchannels of Arbitrary Cross-Section,
ASME J. Fluids Eng., 128, pp. 10361044.
2 Whitesides, G. M., 2006, The Origins and the Future of Microfluidics, Nature London, 442, pp. 368372.
3 Grushka, E., McCormick, R. M., and Kirkland, J. J., 1989, Effect of Temperature Gradients on the Efficiency of Capillary Zone Electrophoresis Separations, Anal. Chem., 61, pp. 241246.
4 Fletcher, P. D. I., Haswell, S. J., Pombo-Villar, E., Warrington, B. H., Watts, P.,
Wong, S., and Zhang, X., 2002, Micro Reactors: Principles and Applications
in Synthesis, Tetrahedron, 58, pp. 47354757.
5 DeWitt, S., 1999, Microreactors for Chemical Synthesis, Curr. Opin. Chem.
Biol., 33, pp. 350356.
6 Miyake, R., Lammerink, S. J., Elwenspoek, M., and Fluitman, H. J., 1993,
Micro Mixer With Fast Diffusion, MEMS 93, Fort Lauderdale, FL, pp.
248253.
7 Weigl, B. H., Bardell, R. L., and Cabrera, C. R., 2003, Lab-on-a-Chip for
Drug Development, Adv. Drug Delivery Rev., 55, pp. 349377.
8 Becker, H., and Locascio, L. E., 2002, Polymer Microfluidic Devices, Talanta, 56, pp. 267287.
9 Chovan, T., and Guttman, A., 2002, Microfabricated Devices in Biotechnology and Biochemical Processing, Trends Biotechnol., 203, pp. 116122.
10 Ehrfeld, W., 2003, Electrochemistry and Microsystems, Electrochim. Acta,
48, pp. 28572868.
11 Bayraktar, T., and Pidugu, S. B., 2006, Characterization of Liquid Flows in
Microfluidic Systems, Int. J. Heat Mass Transfer, 49, pp. 815824.
12 Peng, X. F., Peterson, G. P., and Wang, B. X., 1994, Frictional Flow Characteristics of Water Flowing Through Rectangular Microchannels, Exp. Heat
Transfer, 7, pp. 24964.
13 Xu, B., Ooi, K. T., and Wong, N. T., 2000, Experimental Investigation of
Flow Friction for Liquid Flow in Microchannels, Int. Commun. Heat Mass
Transfer, 27, pp. 11651176.
14 Ren, C. L., and Li, D., 2004, Electroviscous Effects on Pressure-Driven Flow
of Dilute Electrolyte Solutions in Small Microchannels, J. Colloid Interface
Sci., 274, pp. 319330.
15 Judy, J., Maynes, D., and Webb, B. W., 2002, Characterization of Frictional
Pressure Drop for Liquid Flows Through Microchannels, Int. J. Heat Mass
Transfer, 45, pp. 34773489.
16 Tuckerman, D. B., and Peace, R. F. W., 1981, High-Performance Heat Sinking for VLSI, IEEE Electron Device Lett., 2, pp. 126129.
17 Mala, G. M., and Li, D., 1999, Flow Characteristics of Water in Microtubes,
Int. J. Heat Mass Transfer, 20, pp. 142148.
18 Pfhaler, J., Harley, J., and Bau, H., 1990, Liquid Transport in Micron and
Submicron Channels, Sens. Actuators, A, A21A23, pp. 431434.
19 Pfhaler, J., Harley, J., Bau, H., and Zemel, J. N., 1991, Gas and Liquid Flow
in Small Channels Micromechanical Sensors, Actuators, and Systems, ASME
Dyn. Syst. Control Div., 32, pp. 4960.
20 Urbanek, W., Zemel, J. N., and Bau, H. H., 1993, An Investigation of the
Temperature Dependence of Poiseuille Numbers in Microchannel Flow, J.
Micromech. Microeng., 3, pp. 206208.
21 Qu, W., Mala, G. M., and Li, D., 2000, Pressure-Driven Water Flows in
trapezoidal Silicon Microchannels, Int. J. Heat Mass Transfer, 43, pp. 3925
3936.
22 Papautsky, I., Brazzle, J., Ameel, T., and Frazier, A. B., 1999, Laminar Fluid
Behavior in Microchannels Using Micropolar Fluid Theory, Sens. Actuators,
A, 73, pp. 101108.
23 Ren, C. L., and Li, D., 2005, Improved Understanding of the Effect of Electrical Double Layer on Pressure-Driven Flow in Microchannels, Anal. Chim.
Acta, 531, pp. 1523.
24 Weilin, Q., Mala, H. M., and Li, D., 2000, Pressure-Driven Water Flows in
Trapezoidal Silicon Microchannels, Int. J. Heat Mass Transfer, 43, pp. 353
364.
25 Guo, Z., and Li, Z., 2003, Size Effect on Microscale Single-Phase Flow and
Heat Transfer, Int. J. Heat Mass Transfer, 46, pp. 149159.
26 Morini, G. L., 2005, Viscous Heating in Liquid Flows in Microchannels, Int.
J. Heat Mass Transfer, 48, pp. 36373647.
27 Bahrami, M., Yovanovich, M. M., and Culham, J. R., 2006, Pressure Drop of
Laminar, Fully Developed Flow in Rough Microtubes, ASME J. Fluids Eng.,
128, pp. 632637.
28 Jiang, X. N., Zhou, Z. Y., Huang, X. Y., and Liu, C. Y., 1997, Laminar Flow
Through Microchannels Used for Microscale Cooling Systems, IEEE/CPMT
Electronic Packaging Technology Conference, pp. 119122.
29 Baviere, R., Ayela, F., Le Person, S., and Favre-Marinet, M., 2005, Experimental Characterization of Water Flow Through Smooth Rectangular Microchannels, Phys. Fluids, 17, p. 098105.
30 Bucci, A., Celata, G. P., Cumo, M., Serra, E., and Zummo, G., 2003, Water
Single-Phase Fluid Flow and Heat Transfer in Capillary Tubes, International
Conference on Microchannels and Minichannels, Vol. 1, pp. 319326, ASME
Paper No. 1037.
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
31 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.
32 Liu, D., and Garimella, S., 2004, Investigation of Liquid Flow in Microchannels, J. Thermophys. Heat Transfer, 18, pp. 6572.
33 Gao, P., Le Person, S., and Favre-Marinet, M., 2002, Scale Effects on Hydrodynamics and Heat Transfer in Two-Dimensional Mini and Microchannels,
Int. J. Therm. Sci., 41, pp. 10171027.
34 Squires, T. M., and Quake, S. R., 2005, Micro Fluidics: Fluid Physics at
Nanoliter Scale, Rev. Mod. Phys., 77, pp. 9771026.
35 White, F. M., 1974, Viscous Fluid Flow, McGraw-Hill, New York, Chap. 3.
36 Bahrami, M., Yovanovich, M. M., and Culham, J. R., 2007, A Novel Solution
for Pressure Drop in Singly Connected Microchannels, Int. J. Heat Mass
Transfer, 50, pp. 24922502.
37 Timoshenko, S. P., and Goodier, J. N., 1970, Theory of Elasticity, McGrawHill, New York, Chap. 10.
38 Yovanovich, M. M., 1974, A General Expression for Predicting Conduction
Shape Factors, AIAA Prog. in Astro. and Aeronautics: Thermophysics and
Space Craft Control, Vol. 35, R. G. Hering, ed., MIT Press, Cambridge, MA,
pp. 265291.
39 Muzychka, Y. S., and Yovanovich, M. M., 2002, Laminar Flow Friction and
Heat Transfer in Non-Circular Ducts and Channels Part 1: Hydrodynamic
Problem, Proceedings of Compact Heat Exchangers, A Festschrift on the 60th
Birthday of Ramesh K. Shah, Grenoble, France, pp. 123130.
40 Shah, R. K., and London, A. L., 1978, Laminar Flow Forced Convection in
Ducts, Supplement to Advances in Heat Transfer, Academic, New York.
41 Erickson, D., Sinton, D., and Li, D., 2003, Joule Heating and Heat Transfer in
Polydimethylsiloxane Microfluidic Systems, Lab Chip, 3, pp. 141149.
42 McDonald, J. C., Duffy, D. C., Anderson, J. R., Chiu, D. T., Wu, H., Schueller,
O. J. A., and Whiteside, G., 2000, Fabrication of Microfluidic Systems in
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
S. M. Ghiaasiaan1
e-mail: mghiaasiaan@gatech.edu
G. W. Woodruff School of Mechanical
Engineering,
Georgia Institute of Technology,
Atlanta, GA 30332-0405
The laminar pulsating flow through porous media was numerically studied. Twodimensional flows in systems composed of a number of unit cells of generic porous
structures were simulated using a computational fluid mechanics tool, with sinusoidal
variations in flow with time as the boundary condition. The porous media were periodic
arrays of square cylinders. Detailed numerical data for the porosity ranging from 0.64 to
0.84, with flow pulsation frequencies of 2064 Hz were obtained. Based on these numerical data, the instantaneous as well as the cycle-average permeability and Forchheimer
coefficients, to be used in the standard unsteady volume-averaged momentum conservation equation for flow in porous media, were derived. It was found that the cycle-average
permeability coefficients were nearly the same as those for steady flow, but the cycleaverage Forchheimer coefficients were significantly larger than those for steady flow and
were sensitive to the flow oscillation frequency. Significant phase lags were observed
between the volume-averaged velocity and the pressure waves. The phase difference
between pressure and velocity waves, which is important for pulse tube cryocooling,
depended strongly on porosity and the mean-flow Reynolds number.
DOI: 10.1115/1.3089541
Keywords: porous media, pulsating flow, permeability, Forchheimer coefficient, phase
shift
Introduction
Pulsating and periodic flows in porous media are poorly understood. An important area of application of such flow is in the
Stirling and pulse tube cryocooler systems. In these systems, microporous structures constitute cold and warm heat exchangers,
and more importantly a regenerator, all of which are subject to a
periodic flow of a gas usually helium.
The fluid-solid interactions in these microporous structures, in
particular, the regenerator, are important with respect to the performance of the cryocooler system. The frictional losses, in particular, adversely affect the overall system performance. A number
of experimental investigations dealing with friction in commonlyused cryocooler regenerator material have been reported in the
recent past 16. These investigations have all dealt with the
macroscopic aspects of friction in porous media subject to periodic flow and have shown that, except at very low frequencies, the
closure parameters that represent friction between the fluid and
solid are different in steady and periodic flows. However, no systematic study of pore-level phenomena during periodic or pulsating flows through porous media has been reported in the past. It is
thus the objective of this investigation to study the pore-level
phenomena in a pulsating flow through a generic two-dimensional
porous medium by numerical analysis, and from there examine
the effects of flow pulsations on the closure parameters that are
encountered in the standard volume-averaged momentum equation for porous media. We have chosen the pulsating flow as an
intermediate step toward the more difficult problem of periodic
flow.
Since pore-level analysis in general is not feasible for the usual
design and analysis calculations for flow through porous media,
model equations are often applied instead. The most widely used
1
Corresponding author.
Present address at Boiling and Two-Phase Flow Laboratory, School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907-1288
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received August 8, 2008; final manuscript received January 5, 2009; published online March 9, 2009. Assoc. Editor:
Neelesh A. Patankar.
2
Geometric Configuration
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. 1 Different porous structure geometries showing a unit cell of continuous porous
structures
u
1
+ u u = P + 2u
t
Governing Equations
1
V
udV
Vf
where V and V f denote the averaging control volume and the fluid
volume contained within the averaging control volume, respectively see Fig. 2. The intrinsic volume-average fluid velocity is
defined as
u f =
1
Vf
udV
Vf
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
ux f
2ux f
ux f
P f
ux f
=
+ ux f
+
x
x2
t
x
Kxx
2Bxxux f ux f
where Kxx and Bxx are the x components of the diagonal tensors K
P f
ux + Bxx,stuxux
x
Kxx,st
L2
+ Bxx,stL
Kxx,st ReL
where
uL
10
P f L
x ux2
11
ReL =
and
u f
u f
+ u f u f = P f + 2u f
t
K
2B u f u f
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
f,n
f,n
f,n
f,n
ui+1
ui1
uif,n uif,n1
Pi1
Pi+1
=
+ uif,n
2x
t
2x
f,n
f,n
ui+1
2uif,n + ui1
u f,n 2Bxx,instt
2
x
Kxx,instt i
uif,nuif,n + Ot,x2
12
where subscript n and superscript i are time step and grid point
indices, respectively. Note that u in the above equation implies ux,
and the subscript x has been left out for convenience.
The calculation procedure for extracting the permeability and
Forchheimer coefficients from the microscopic-level solutions for
pulsating flow is as follows. We note that at the limit of u f 0,
the Forchheimer term in Eq. 7 or Eq. 12 can be neglected.
Thus, for each of the configurations in Fig. 1, the following steps
are followed:
1 Obtain uif,n and Pif,n for i = 1 6 and during n equallyspaced time snapshots.
2 Solve the following equation for the time snapshots covering two pulsation cycles for the low-flow simulations, and
thereby obtain Kxx,instt.
f,n
f,n
ui+1
ui1
uif,n uif,n1
=
+ uif,n
2x
t
f,n
f,n
f,n
f,n
ui+1
2uif,n + ui1
Pi1
Pi+1
+
2
x
2x
u f,n
Kxx,instt i
13
to+1/f
Kxx,insttdt
14
Bxx,insttdt
15
a = 0.4
16
Steady Flow
For CFD simulations, a convergence criterion of 106 was applied to the residuals of the continuity and momentum equations.
Grid independence was examined by using two grid systems with
20 40 and 40 80 per unit cell. The two simulations showed
nearly identical velocity and pressure distributions for the pairs of
simulations with the maximum difference between the volumeaveraged pressure gradients within 1%. For calculation efficiency,
a grid system with 20 40 per unit cell was chosen for all the
subsequent simulations.
Figure 4 represents the typical steady-state streamline patterns.
The streamline patterns are similar for all the unit cells. The flow
details in the unit cells 4, 5, and 6 are virtually identical, indicating a fully developed state. It is seen that for the low Reynolds
number flow, a pair of small vortices appears between the adjacent
square rods. For the high Reynolds number flow, the position of
the vortex center in the cavity is shifted slightly downstream.
The pressure drop characteristics along the flow direction are
presented in terms of the nondimensional pressure gradient in
Fig. 5. The nondimensional pressure gradient undergoes a large
change in the entrance region and reaches an almost constant
value along the flow direction. At very low Reynolds numbers
ReL 1, it is observed that the hydrodynamically fully developed regime is reached just after the second unit cell. However, in
case of high Reynolds numbers, the pressure gradient still varies
slightly from the third to fourth unit cells, and the change in
pressure gradient from one unit cell to the next is decreased as the
fully developed condition is reached further downstream. Similar
to
Bxx,avg = f
to+1/f
to
Fig. 5 Nondimensional pressure gradient along the flow direction with Reynolds number for the steady 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
L2
Kxx,st
L2
Kxx,avg
Bxx,stL
Bxx,avgL
0.64
0.6975
0.75
0.75a
0.75b
0.7975
0.84
153
105
75
75
75
55
41
159
113
79
80
80
59
47
0.071
0.058
0.049
0.049
0.049
0.041
0.036
0.121
0.108
0.096
0.078
0.088
0.085
0.077
20 Hz.
64 Hz,
Fig. 6 Nondimensional pressure gradient as a function of Reynolds number for steady flow
Pulsating Flow
Fig. 7 The Forchheimer term as a function of Reynolds number for steady 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
Fig. 9 Convergence check for = 0.75 with 20 40 nodes per unit cell for
a velocity and b pressure
tion, and the waves are flattened with decreasing porosity. The
results also show that the cycle-average instantaneous velocities at
the measured points increase as the porosity decreases, and as the
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
17
P = P,pt6 P,pt1/5
18
VP = VP,pt1 + VP,pt6/2
19
Fig. 14 Variation in the instantaneous permeability coefficients for different porosities f = 40 Hz calculated from Rem,L
= 0.11
Fig. 15 Variation in the instantaneous Forchheimer coefficients for different porosities f = 40 Hz calculated from Rem,L
= 0.11 and 560
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
Nomenclature
a normalized amplitude of velocity
A constant
B Forchheimer coefficient 1/m
Bxx,avg cycle-average Forchheimer coefficient m2
Bxx,inst instantaneous Forchheimer coefficient along
the x direction for pulsating flow m2
Bxx,st Forchheimer coefficient along the x direction
for steady flow m2
B Forchheimer coefficient tensor 1/m
D side of square rod m
f frequency Hz
I identity tensor
K permeability coefficient m2
Kxx,avg cycle-average permeability coefficient m2
Kxx,inst instantaneous permeability coefficient along the
x direction for pulsating flow m2
Kxx,st permeability coefficient along the x direction
for steady flow m2
K permeability tensor m2
L length of unit cell m
P static pressure N / m2
P f intrinsic volume-average fluid static pressure
over V f N / m2
ReL Reynolds number based on a unit cell length L
in steady flow
Rem,l Reynolds number based on mean velocity and
unit cell length pulsating flow
t time s
T period of pulsation s
u, U velocity m/s
u superficial volume-average fluid velocity m/s
u f intrinsic volume-average fluid velocity m/s
V averaging control volume m3
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
porosity of medium
molecular viscosity kg/ m s
kinematic viscosity m2 / s
phase angle deg
density kg/ m3
phase shift of velocity deg
phase shift of pressure deg
phase shift between velocity and pressure deg
nondimensional pressure gradient
Subscripts
avg
e
exit
i
in
inst
m
st
x, xx
cycle-average value
effective
exit
grid point index
inlet
instantaneous value
mean value
steady state
x direction
Superscripts
f fluid
n time step index
= tensor
References
1 Ju, Y. L., Wang, C., and Zhou, Y., 1998, Numerical Simulation and Experimental Verification of the Oscillating Flow in Pulse Tube Refrigerator, Cryogenics, 38, pp. 169176.
2 Hsu, C. T., Fu, H. L., and Cheng, P., 1999, On Pressure Velocity Correlation
of Steady and Oscillating Flows in Regenerators Made of Wire-Screens,
ASME J. Fluids Eng., 121, pp. 5256.
3 Nam, K., and Jeong, S., 2002, Experimental Study on the Regenerator Under
Actual Operating Conditions, Adv. Cryog. Eng., 47, pp. 977984.
4 Jeong, S., Nam, K., and Jung, J., 2002, Regenerator Characterization Under
Oscillating Flow and Pulsating Pressure, Cryocoolers, 12, pp. 531537.
5 Cha, J. S., Ghiaaisiaan, S. M., Desai, P. V., Harvey, J. P., and Kirkconnell, C.
S., 2006, Multi-Dimensional Flow Effects in Pulse Tube Refrigerators,
Cryogenics, 46, pp. 658665.
6 Clearman, W. M., Cha, J. S., Ghiaasiaan, S. M., and Kirkconnell, C. S., 2008,
Anisotropic Steady-Flow Hydrodynamic Parameters of Microporous Media
Applied to Pulse Tube and Stirling Cryocooler Regenerator, Cryogenics, 48,
pp. 112121.
7 Whitaker, S., 1973, The Transport Equations for Multi-Phase Systems,
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
Aerodynamic Analysis of a
Vehicle Tanker
Ramon Miralbes Buil
Luis Castejon Herrer
Research Group in Vehicles and
Road Safety (VEHIVIAL),
CPS,
University of Zaragoza,
Zaragoza 50012, Spain
The aim of this article is the presentation of a series of aerodynamic improvements for
semitrailer tankers, which reduce the aerodynamic resistance of these vehicles, and,
consequently, result in a positive impact on fuel consumption, which is substantially
reduced (up to 11%). To make the analysis the computational fluid dynamics (CFD)
methodology, using FLUENT, has been used since it allows simulating some geometries
and modifications of the geometry without making physical prototypes that considerably
increase the time and the economical resources needed. Three improvements are studied:
the aerodynamic front, the undercarriage skirt, and the final box adaptor. First they are
studied in isolation, so that the independent contribution of each improvement can be
appreciated, while helping in the selection of the most convenient one. With the aerodynamic front the drag coefficient has a reduction of 6.13%, with the underskirt 9.6%, and
with the boat tail 7.72%. Finally, all the improvements are jointly examined, resulting in
a decrease of up to 23% in aerodynamic drag coefficient. DOI: 10.1115/1.3077135
Keywords: aerodynamic, undercarriage skirts, tank, vehicle tanker, nozzle, boat-tails,
drag, penetration, consumption, vehicle, articulated
Introduction
expenditure, to savings of up to 12,000 l in fuel, and to the lowering of carbon dioxide emissions by almost 30 tons 3.
It should also be kept in mind that aerodynamic improvements
not only influence consumption and the power required by vehicles; they also increase their stability and decrease the effects of
aerodynamics on other vehicles traveling on the same road, as
well as diminishing splashing and improving the vehicles resistance to lateral winds, etc.
Therefore, an aerodynamic analysis of tankers was considered
essential in order to verify how certain elements improve the aerodynamics of the whole and thus contribute to the vehicles overall
performance.
New aerodynamic elements applicable to vehicle tankers and
which considerably reduce the vehicles aerodynamic resistance
were also developed.
Before undertaking a detailed study of the aerodynamic improvements implemented, the vehicles energy budget has to be
calculated using Clarks equation 4, as follows:
P = 2 Cx At V3 + M G V + M G V sin
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
area of the box, although the use of one type of boat-tail or the
other depends on the type of box, the maximum size of the vehicle, etc.
2.2 Aerodynamic Fairings. These aerodynamic elements are
fixed on the upper front area of the semitrailer see Fig. 3 and
redirect the air flow that appears in that area toward the box,
keeping it out of the king-pin cavity and thus involving an aerodynamic improvement of up to 3% 7.
2.3 Aerodynamic Undercarriage Skirt. This part prevents
the entry of air in the lower part of the vehicle, both in the tractor
area and the semitrailer area, decreasing air flow in these areas,
which due to the wheels and several parts of the vehicles lower
area would otherwise generate a turbulent flow that would adversely affect the vehicles aerodynamics. A decrease of up to 7%
8 can be achieved this way see Fig. 4.
Fig. 2 Boat-tails
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
namic accessories that may come with it: spoiler, lateral air flow
streamliners, and bumper spoiler. These elements were not used,
as we only intended to affect the semitrailer. The reason for not
using them is that since these elements are cab accessories, not all
cabs employ them, and also that, due to the fact that they have
been studied and optimized in depth, these elements have many
different configurations and shapes.
The air inlet of the vehicles radiator was also left out of the
modeling, because taking this factor into account would have required a very elaborate model of the vehicle, which would, in
turn, have involved numerical complications that would have
made calculation unfeasible. Besides, the contribution of this part
to the global model is relatively small and is not under study here
9,10,12.
As far as the tanker is concerned, the turn signals, sleeves, and
tubes of the lower area of the tanker, etc., were removed.
The several elements, which comprise the tankers lower area,
were not removed, since this is one of the areas under study,
despite the fact that this entails an increase both in the complexity
of the model and in the number of elements that comprise it.
In relation to the tire area, the wheel pitch and the empty spaces
within it, along with its dimensions, were maintained in order to
correctly model the contribution of these elements to the vehicle,
since the tire area significantly affects the vehicles overall aerodynamics and is one of the areas subject to study 9,10,12. It may
be worth pointing out that a symmetrical model was used to reduce the necessary computational cost. Consequently, only half
the vehicle was modeled.
3.2 Discretization Through Volumetric Elements. In order
to perform modeling using volumetric elements, the vehicles
CAD was taken as the starting point and introduced into a volume,
which is to simulate the air on the outside of the vehicle. The
following volumein keeping with the recommendations of other
authors 9,10,12 and those of the programwas considered for
this purpose 12.
In theory, there should be a distance of 2000 mm in the air
volume of a vehicle similar to the one selected between the front
part of the vehicle and the air inlet, and a distance of 10,000 mm
between the end of the vehicle and the air outlet in order to perfectly collect the slipstream generated by the vehicle. Regarding
both lateral and superior distances, these should be at least 4000
mm in relation to the outermost tangent surface.
The dimensions of this volume are such that the air circulating
in the areas furthest from the vehicle is not altered by the latter,
which implies that the imposition of contour conditions on the
outer area of the volume will not be altered by the vehicle, as this
would involve a modification of the models air flows and, consequently, an error in the results.
APRIL 2009, Vol. 131 / 041204-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
The next step was to extract the vehicles volume from the total
volume, with which we could obtain a negative image of the
vehicle with regards to how discretization will be carried out
To do so, a tetrahedral mesh with lineal Lagrange elements was
chosen, because it allows moderately simple vehicle discretization
and an optimum degree of approximation in relation to computational cost.
First of all, discretization of the outer surfaces was performed
both on the inside and on the outside. The size of these surface
elements depends on the area in which discretization is carried out
9,10,12, on its magnitude and on previous experience, although
it was attempted to reduce the number of elements to the minimum, with the proviso that they be as equilateral as possible in
order to ensure convergence and reduce time and calculation requirements. A more elaborate discretization was performed in
those areas where greater speed gradients were expected, allowing
for improved simulation of air flow in those areas: the wheel
pitch, the front of the vehicle, the area between the cab, and the
tank and the rear part of the vehicle, where detachment of the
outer layer will occur.
After the discretization of the surfaces was performed, boundary layers 13 were imposed. This procedure allows us to per-
form fine meshing of the zones immediately surrounding the surfaces of the vehicle in order to faithfully model the vehicles outer
layer. To this end, an initial element size of 6 mm and a boundary
layer of three elements with a mesh growth factor of 1.2 9,10,12
was chosen. The reason for building such a mesh was to try to
correctly approximate the outer layer, despite the fact that this
entails a substantial increase in the number of elements and an
increase in computational cost; however, this is justifiable and was
necessary for the correct modeling of the outer layer of fluid.
Following this, discretization of the remaining air volume was
carried out; in order to do so, meshing of the outer area with a
larger size element 1000 mm was performed, since attached flow
will be preponderant in this area and there will be no flow alterations or speed gradients. In order to perform discretization, a size
function was used, allowing control of the growth size of the
elements from the vehicle to the outside area. To achieve this, a
maximum element size of 1000 mm and a 1.2 growth factorthe
maximum recommended by the programwas imposed. In this
way, discretization with the most equilateral surfaces possible can
be achieved see Figs. 7 and 8.
Fig. 8 Mesh of the baseline equipment trailer middle zone and tractor
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.3 Aerodynamic Calculation. In order to perform aerodynamic calculation, several tasks have to be carried out prior to the
calculation proper:
3.3.1 Definition of Solution Options. A segregated solution
modelrecommended for this type of calculationswas chosen
in order to lower computational requirements. This model uses a
double precision calculation recommended for aerodynamic calculation of vehicles 14,15 for a stationary steady model based
on node calculation and is recommended by the program in order
to achieve better results with this type of calculation.
3.3.2 Definition of the Turbulence Model. A turbulence model
based on k- type realizable Reynolds average Navier-Stokes
RANS equations 1416 and a wall treatment with nonequilibrated flow were chosen 17. The choice of this turbulence model
was mainly due to the fact that nowadays it is the most commonly
used one for aerodynamic calculations and because it also offers
an optimum degree of approximation in regard to computational
requirements.
3.3.3 Definition of Contour Conditions. Since there are several
areas, different contour conditions were considered for each of
them. An air flow at a constant speed 30 m/s in a tangential
direction 7 with very low turbulence 0.1% turbulence intensity
was introduced for the entry boundary, which means that any flow
that there may be will be laminar. Regarding the upper, lateral,
and rear areas, these were simulated as being open to the outside,
with a pressure of 1 atm and an air flow of 30 m/s in the relevant
direction.
The ground was modeled as a smooth surface traveling at the
same speed as the fluid 30 m/s 18, and the wheels were considered smooth and rotational surfaces, coherent with the speed of
the vehicle 17.
Symmetry was imposed on the interior wall and the vehicles
walls were considered fixed and smooth walls without friction.
This was due to the fact that the contribution of friction between
the fluid and the surface in relation to the global model represents
less than 10% of the resistant contribution and that the implementation of this contribution would entail a substantial increase in
computational cost and might even have lead to divergent results.
Nonetheless, it would be useful to establish, very precisely and by
means of experiments, the friction coefficient of each of the models superficial finish types in order to establish their contribution
to the aerodynamic resistance, which would not be substantially
modified by variations in the vehicles geometry, since the latter
depends mainly on the choice of superficial finishes.
Thus, simulation of the movement of the vehicle is possible, but
from the relative viewpoint of the tanker.
3.3.4 Selection of Options for Solving the Equations. A
13 type resolution algorithm was chosen to solve the
SIMPLE-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 Total resistant forces and aerodynamic resistance coefficients of all the cases considered
flow, turbulence energy maps, and particle flow maps draft lines
were analyzed in order to perform this analysis, as well as an
analysis of the flow of particles in motion, which was not included
here due to its being a moving image.
It must be pointed out that in order to keep the article relatively
short, only the more representative images of the vehicle were
included.
As shown in Fig. 10 and predicted earlier, the low speed areas
are the lower part of the vehicle, the king-pin cavity area and the
slipstream of the tank; i.e., these are the areas with greater turbulence and higher pressures. This figure shows the longitudinal
speed of the vehicle. It shows the wake in the rear part of the
vehicle and indicates a high aerodynamical resistance. It also
shows, in the bottom part of the vehicle past the tractor, a zone
with low velocity and the same in the wheel zone of the tanker.
This indicates high aerodynamical losses too. Analyzing Figs. 11
and 12 that show the k turbulence parameter and then the zones
with high aerodynamical losses, these match with the zones previously mentioned and with Fig. 9, so these zones should be
modified to improve the aerodynamics. These zones also have
whirls that indicate turbulence see Fig. 10 and 13 and cause
aerodynamical losses, usually because the geometry of the vehicle
is discontinuous, and the attached flow changes to a turbulence
one. In Fig. 10 there is a zero velocity zone in the gap between the
tractor and the trailer and it matches with a high losses zone see
Figs. 11 and 15.
We also notice in Figs. 10 and 13 a stagnation zone in the front
area that is one of the zones with high aerodynamical losses of
the tractors, the detachment of the outer layer above the cab in
the vehicles rear part and the lower part of the tractor, as well as
coupling of the flow in the upper part of the tank all was
predicted in Fig. 9.
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 Longitudinal speed of the original configuration at the vehicles plane of symmetry
Fig. 11
Fig. 12
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. 13 Flow lines and associated speeds of the original configuration at the vehicles plane of
symmetry for the cab and rear areas
previously established flow hypotheses were correct, and consequently that modification of the king-pin cavities and the wheels
would be necessary in order to reduce the vehicles aerodynamic
losses. Regarding the slipstream area of the vehicle, no initial
modifications had to be undertaken, since the aerodynamic profile
and the level of turbulence are correct and flow detachment cannot
be avoided in that area; the configuration with the box in the rear
will require modification, since it generates greater flow detachment and more turbulence, which ought to be prevented.
Fig. 14 Longitudinal speed of the original configuration at a plane parallel to the plane of
symmetry and in the tire area
Fig. 15 Longitudinal speed of the original configuration at a plane parallel to the plane of
symmetry tangential to the tank
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. 16
traps in the lower front part of our vehicle is also impossible, and
therefore, this improvement cannot be implemented either.
The main aerodynamic elements subject to study will be undercarriage skirts. An aerodynamic modification based on the boattail will also be designed for the vehicle with a rear box configuration in order to prevent the detachment of the outer layer of this
zone, which we will call aerodynamic box adaptor; on the
other hand, an aerodynamic improvement for the vehicles front
areain order to prevent the generation of vorticity and the detachment of the outer layer of the king-pin cavitywill be designed on the basis of a modification of the aerodynamic fairings
and the lateral air flow streamliners, but will be applied to the tank
and will be named aerodynamic front.
4.1 Undercarriage Skirts. At present, all newly built semitrailers have to implement lateral protection measures in compliance with EEC guideline 98/297. This directive regulates the use
of lateral protectors. The use of these protective elements, in their
horizontal spar configurationthe most widely used one
nowadaysleads to aerodynamic penalization, since they favor
detachment of the outer layer due to their discontinuities with the
tank.
Furthermore, the cavities of the wheel area allow air to penetrate this area and therefore generate an aerodynamic penalization, which increases in proportion to the increase in the air flow
running through that area; therefore, the use of aerodynamic undercarriage skirts will be implemented to reduce the air flow laterally penetrating the area.
Consequently, an aerodynamic improvement was suggested in
order to decrease the aerodynamic penalization caused by these
elements.
The aerodynamic improvement in question is based on the reduction of free edges and geometries capable of generating a detachment of the outer layer in the area between the tractor and the
wheels of the vehicle, since in this area there are many elements
that generate vorticity and aerodynamic losses.
To this end, we suggest a continuous surface going from the
front area of the semitrailer to the rear area, comprising of three
clearly differentiated areas, and additionally, made out of different
materials.
1 The area between the tractor and the parking holder A:
this area is to be made of plastic materials, since it will not
be subjected to important mechanical stress. Its function is
to avoid air penetrating the lower area of the tank from the
lateral areas, mainly issuing from the tractor.
2 The area between the parking holder and the wheels B:
this part was designed according to EEC directive 98/297
and allows the use of said area as a lateral protector in
compliance with previous specifications.
3 The area of the semitrailers wheels C: this part covers
both the cavities located above the vehicles wheelguards
and those between the semitrailers wheels, as these cavities
generate vorticity and aerodynamic losses.
We can see the aerodynamic improvements in Fig. 17.
When we analyze the results obtained we will see that k turbulence is reduced in the lower part of the vehicle see Figs. 16 and
18. This indicates a decrease in turbulence in the area, and therefore a decrease in aerodynamic losses.
The trajectories of particles Fig. 19 follow the direction of the
flow and very small amounts of fluid penetrate the vehicles lower
area, as might have been expected, while the curved part in area A
performs its task, preventing air flow in that area Fig. 20. This
indicates that the turbulence flow first has decreases and then
aerodynamical losses. This is because the underskirt avoids abrupt
geometry changes and then the whirls that are generated are lower
see Figs. 18 and 20 versus Fig. 11.
On the other hand, we may see that the final curved shape of
area C of the undercarriage skirt redirects flow and reduces slipstream in that area see Fig. 20.
Regarding the forces acting on the vehicle see Table 1 and the
aerodynamic resistance coefficient, we see that the vehicles aerodynamic coefficient is 0.56, i.e., 9% less, which means that the
aerodynamic benefit of this part is considerable and has no impact
in terms of weight increase approximately 77 kg
4.2 The Aerodynamic Front. The king-pin area is a zone in
which detachment of the outer layer and high vorticity occur, as
shown in Fig. 11. Consequently, it is an area with high aerodynamic losses, and therefore we suggested an aerodynamic improvement in order to reduce the king-pin cavity and redirect the
air in that area to prevent whirlwinds and detachments of the outer
layer.
The aerodynamic front was designed so as not to interfere with
the vehicles maneuvers or with access to the area when it comes
to connecting it to the tractor. With this in mind, geometry such as
Fig. 17 Aerodynamic undercarriage skirt in a configuration with the box in the lower part versus linebase one
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
Fig. 20 Detail of areas A right and C left of the aerodynamic undercarriage skirt
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. 22 Vehicle tanker with an aerodynamic front versus linebase one. Top view.
Fig. 23 k turbulence parameter on the plane of symmetry and on a vertical plane 2 m above
the ground in the king-pin area for a vehicle with an aerodynamic 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. 25 k turbulence parameter on a horizontal plane 2 m above the ground and at the plane of
symmetry for a rear box configuration without improvements
detachment of the outer layer of the said area and minimize overall aerodynamic losses were used.
Since these elements are not subject to significant mechanical
stress, they can also be made out of plastic, which minimizes the
total weight contribution 57 kg in 3 mm PVC
Fig. 26 k turbulence parameter at the plane of symmetry for a rear box configuration without improvements
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. 28 Detail of the flow lines at the plane of symmetry in a vehicle with final box adaptor
Fig. 29
4.4 Global Contribution of all the Aerodynamic Improvements Shown. To complete the analysis, an aerodynamic calculation of both box configurations including all the aerodynamic improvements at our disposal was carried out in order to verify the
overall contribution of these improvements to the vehicle.
To this end, all previous aerodynamic improvements were
adapted with the purpose of jointly integrating them with the vehicle and achieving a perfect union of the different parts in order
to avoid detachment of the outer layer and make the vehicle even
more aerodynamic.
Consequently, these were the geometries chosen.
k turbulence parameter at the plane of symmetry in a vehicle with final box adaptor
Fig. 30 Front area of a vehicle with all the aerodynamic improvements, regardless of box configuration and equipped with a
rear box in the rear area
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. 31 Speed at the plane of symmetry in a vehicle with lower box configuration and all the
improvements
If we examine the results obtained for the lower box configuration, we find a 0.52 aerodynamic coefficient, i.e., 15.85% less,
due both to the contribution of the undercarriage skirt and that of
the aerodynamic front, as well as to the joint integration of both
elements. This can be appreciated in Figs. 3133, in which we
may notice the reduction of the k vorticity parameter in the undercarriage skirt area and the king-pin cavity. We also notice the
movement of particles, which tends to produce less vortices and,
consequently, to generate less turbulence. Thus the aerodynamic
improvements are expected to prove effective see Table 1.
The implementation of all the aerodynamic improvements for
the rear box configuration yields a 0.505 aerodynamic coefficient,
i.e., 23% less in relation to the original rear box configuration, as
might have been expected see Fig. 34. This aerodynamic coefficient is noticeably inferior to that of the vehicle with the lower
box configuration including all the aerodynamic improvements,
due to the fact that the configuration used for the rear box is more
aerodynamic, although the improvement obtained amounted only
to 3%.
As far as turbulence is concerned, and as a result of the aerodynamic improvements implemented in these areas, we notice a
decrease in the value of the k turbulence parameter in the slipstream, the wheel area and the king-pin cavity Fig. 33 versus Fig.
14.
We may also notice, in Fig. 35 versus Fig. 13, that there are less
vortices in the flow lines and that flow is redirected.
Figures 33 and 36 versus Fig. 13 shows that the turbulence in
the bottom zone and in the wheel zone of the trailer has been
reduced, because the underskirt decreases the geometry changes
and redirects the air flow properly see Figs. 32 and 35 versus Fig.
13.
Economical Analysis
Fig. 32 Flow lines on a vertical plane 1 m above the ground in a vehicle with lower box
configuration and all the improvements
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. 33 k turbulence parameter at the plane of symmetry in a vehicle with a lower box and the
improvements
ers. For the cost analysis a price of 2.1 $/gal 11/09/2008 has
been used and a vehicle consumption of 16 gal/100 miles and
100.000 miles/year for a standard vehicle. This indicates an annual consumption of 16,000 gal and a cost of 33,600 $/year. As
for the improvements, the repercussion of them on the consumer
Fig. 34 Speed at the plane of symmetry in a vehicle with rear box configuration and all the
improvements
Fig. 35 Flow lines on a vertical plane 1 m above the ground in a vehicle with a rear box and all
the improvements
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. 36 k turbulence parameter on a vertical plane 1 m above the ground in a vehicle with a
rear box and improvements
Table 2 Economical analysis of the improvements
the fuel consumption is cut down and the annual saving can be
$2100, so the return on investment ROI is approximately 0.5
years and the economic gain is clearly large. all economic results
are in Table 2. There are also environmental benefits with the
CO2 emissions being reduced by 2.73 Tm during a year.
Conclusions
References
1 Morelli, A., 2000, A New Aerodynamic Approach to Advanced Automobile
Basic Shapes, SAE Technical Paper No. 2000-01-0491, p. 104.
2 Wood, R. M., and Bauer, S. X. S., 2003, Simple and Low Cost Aerodynamic
Drag Reduction Devices for Tractor-Trailer Trucks, SAE Paper No. 2003-013377.
3 Bachman, L. J., Erb, A., and Bynum, C. L., 2005, Effect of Single Wide Tires
and Trailer Aerodynamics on Fuel Economy and NOx Emissions of Class 8
Line-Haul Tractor , SAE Paper No. 05CV45.
4 Rutledge, W. H., Polansky, G. F., and Clark, E.L., 1988, Aerodynamic Design
and Performance of a Bent-Axis Geometry Vehicle, J. Spacecr. Rockets,
254, pp. 257262.
5 Wood, R. M., 2004, Advanced Aerodynamic Trailer Technology.
6 2007, FreightWings web page, http://www.freightwing.com
7 2007, Nosecones web page, http://www.nosecone.com
8 Graham, S., and Bigatel, P., 2001, Freight Wing Trailer Aerodynamics, Technical Report No. JA-132.
9 Kleber, A., 2005, Simulation of Air Flow Around an OPEL Astra Vehicle
With Fluent, International Technical Development Center, Adam, Opel AG,
Journal Articles by FLUENT Software users, www.fluent.com/solutions/articles/
ja132.pdf
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
20
21
22
23
24
25
26
Dr. Ramon Miralbes Buil obtained a degree in Industrial Engineering from Zaragoza University in 2005 (Spain), with specialization
in Vehicle and Machinery Design. From this time on, he continued his Ph.D. studies in the area of New Materials for Road Transport
in this university. In February 2008 he finished his Ph.D. studies with the thesis: Flow-Termo-Mechanical Analysis of Cryogenic
Tanker Vehicles and is working on several research projects involving vehicle strength and rigidity performance simulation by the
FEM. He also works as an Associate Professor in the area of Design and Graphic Expression at the University of Zaragoza.
Dr. Luis Castejon Herrer graduated with a degree in Mechanical Engineering (University of Zaragoza, Spain, 1992) and a Ph.D. in
Mechanical Engineering (University of Zaragoza, Spain, 1998). Dr. Luis Castejon Herrer is a full Professor in Transportation Engineering at the University of Zaragoza. Since 1991, he has been involved in several projects in the areas of Vehicle Dynamics and New
Materials for vehicles, buses, railways, and ground transportation. He has written over 120 papers in technical journals and international congresses and 9 books, published in Europe. He has been a researcher in the European BRITE-EURAM Project MULTEXCOMP and he has collaborated in other European projects such as SEABUS and AFICOSS. He has participated in several R + D
Spanish projects. He developed his Ph.D. thesis on the simulation of roll over of metallic and composite buses by means of FEM,
obtaining the Doctoral Special Award of the University of Zaragoza.
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
R. W. Hewson1
School of Mechanical Engineering,
University of Leeds,
Woodhouse Lane,
Leeds, LS2 9JT UK
e-mail: r.w.hewson@leeds.ac.uk
Introduction
Ca 0,102
where h0 is the residual film thickness, r is the constant meniscus radius, and Ca is the capillary number given by Ca= U /
is the fluid viscosity, U is the substrate speed, and is the interfacial surface tension.
Ruschak 7 used the finite element method FEM to obtain an
empirical expression governing the residual film thickness to gap
ratio for higher capillary numbers than those for which the
Bretherton expression is valid.
h0
= 0.54Ca,
h =
Ca 102,101
where h= is the initial fluid filled gap. One limitation of Ruschaks equation is that although it predicts the residual film thickness accurately, the nonconstant meniscus radius means that Eq.
2 is not able to calculate the local pressure at the point of film
formation, an important parameter, and one that is often used as a
pressure boundary condition in a range of different coating and
lubrication problems. This is noted by Thompson et al. 8 who
compared the finite element method with the lubrication approximation for reservoir fed rigid roll coating.
Coyne and Elrod 9 derived a set of four ordinary differential
equations ODEs to describe the process of film formation by
solving the NavierStokes equations along the free surface. The
model assumes a quadratic velocity profile to describe the flow
perpendicular to the free surface. One limitation of the model is
the stiffness of the set of equations around the free surface stagnation point occurring at 3H. Despite this, the method has been
used by some investigators as a free surface boundary condition
1
Corresponding author.
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received July 24, 2008; final manuscript
received January 15, 2009; published online March 11, 2009. Assoc. Editor: Malcolm
J. Andrews.
Model 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
Free surface
Phase 1
Film
Phase 2
(r, )
U
Fig. 1 Illustration of the fluid film forming process
4 = 0
Moving Substrate
=0
=1
at = 0
at =
Flux = 1
Fig. 2 Definition of the local coordinate system r , , the local
tangent to the free surface = 0, and the moving substrate
=
1 =
at =
1 2
=0
r2 2
at = 0
For the case of a free surface, a moving surface, and zero net flux,
the following equation is obtained:
1
=1
r
2 = r
Stationary surface, = , = 1,
Free surface,
= 0, = 0,
1 2
r 2 2 = 0.
Free surface,
= 0, = 0,
1 2
r 2 2 = 0.
1
r
(a)
= 0.
Moving surface, = , = 0,
(b)
10
Free surface,
= 0, = 0,
1 2
r 2 2 = 0.
1
r
= 1.
Moving surface, = , = 1,
1
r
= 1.
(c)
Fig. 3 Graphical representation of the addition of 1 and 2: note the free surface stagnation point predicted by the
addition of the two stream functions 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
dx
= cos
ds
(a)
(b)
Fig. 4 COMSOL MULTIPHYSICS finite element implementation of the
free surface problem
four first order ODEs with respect to distance along the free
surface s, one of which is Eq. 10 p / r = dp / ds, the other
three are as follows:
d
= Cap
ds
11
dh
= sin
ds
12
13
80
20
70
60
15
50
10
40
30
20
0
10
0
10
20
(a)
100
80
60
40
20
10
30
(b)
25
20
15
10
x
3
6
5
2.5
1.5
1
0.5
0
0
0.5
2
3
(c)
12
10
6
x
(d)
Fig. 5 Free surface meniscus shape, pressure distribution and constant radius at = 90 deg, as obtained by the analytical
model Eqs. 1013
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
ODE
Bretherton
Ruschak
CE
FEM
10
3
10
10
(a)
10
Capillary number, Ca
10
0.9
0.8
10
(d/ds)/(C1)
2.1 FEM. The FEM was used to compare with the model
results. The simulation was undertaken in COMSOL MULTIPHYSICS
using the arbitrary Lagrangian Eulerian ALE method. The fluid
was modeled as Stokes flow with the free surface curvature obtained by calculating the rate of change with distance along the
free surface of the angle the free surface makes with the moving
substrate an implementation that can be achieved in COMSOL by
defining the free surface angle as an additional variable along the
boundary. The simulation used to obtain the results had
27,733DOF and solution times were of the order of 10 min; the
mesh and an example solution are shown in Fig. 4. The velocity
inlet boundary condition on the left hand side of the domain is
quadratic, consistent with the lubrication approximation. The flow
rate across this boundary was varied until the angle the free surface made with the symmetry boundary condition along the top
of the domain was 90 deg. The main results extracted from the
FEM simulations where the final film thickness and the free surface profile. To facilitate the analysis the final film thickness was
solved for with the gap h= kept fixed, this is accounted for in
Sec. 3 were the results have been appropriately scaled.
0.7
0.6
0.5
ODE
CE
FEM
0.4
3
10
(b)
10
10
Capillary number, Ca
10
Fig. 6 Results of the model, compared with Bretherton equation, Ruschak equation, and Coyne and Elrod CE results; for
a gap, h=0 at = 90 deg and b scaled meniscus radius of
curvature d / ds / h= 1
the possibility that the analysis may be able to predict the film
profile at higher capillary numberssomething already shown by
examining the free surface profile data.
Conclusions
Acknowledgment
It is a pleasure to acknowledge Yeaw Chu Lee for his interesting discussions on Stokes flow.
References
1 Greener, J., and Middleman, S., 1979, Theoretical and Experimental Studies
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
3
4
5
6
7
8
9
10
11
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
Wodzimierz Wrblewski
e-mail: wlodzimierz.wroblewski@polsl.pl
Sawomir Dykas
Institute of Power Engineering and
Turbomachinery,
Silesian University of Technology,
Gliwice, Poland
Andrzej Gardzilewicz
Institute of Fluid-Flow Machinery,
Polish Academy of Science,
Gdansk, Poland
e-mail: gar@imp.gda.pl
Michal Kolovratnik
Czech Technical University in Prague (CTU),
Prague, Czech Republic
e-mail: kolovrat@fsid.cvut.cz
Introduction
necessary to take into account the size change of the particles and
the solution concentration in the droplets. The condensation model
requires knowledge of many physical and chemical properties of
the salt solution in steam. It is often very hard to find enough
information in the references, as well as from the feedwater analysis in the power plant.
Two processes have to be considered during steam expansion in
the turbine: homogeneous and heterogeneous condensations.
Steam impurities, necessary to initiate heterogeneous condensation, occur in the form of suspension of particles, which are insoluble or soluble in water. The chemical and physical properties
of a suspension influence the size of the droplets nucleated on the
particles and the nucleation process.
It is possible to model various types of condensation processes,
i.e., homogeneous and heterogeneous ones on soluble and insoluble particles. The necessary step in the nonequilibrium wet
steam flow modeling is the validation of the numerical models
against the experimental results. The development of numerical
methods for wet steam flow calculations has to allow for real gas
properties, which, in fact, is the fundamental one for correct calculations of the condensation process. The correct calculation of
the steam flow for low pressures close to the saturation line and
for high pressures needs to implement the real gas equation of
state into numerical algorithm. The relation between thermal parameters of state is nonlinear within these ranges and does not
follow the ideal equation of state. Moreover, in the case of wet
steam, a two-phase flow model has to be considered.
The developed steam condensing flow models have been validated against experimental data for simple 2D test cases. Usually,
they are the Laval nozzles in which 1D flow form dominates. All
available experimental investigations relate to the homogeneous
condensation, whereas in real turbine the heterogeneous condensation is predominant. Therefore, validation of the numerical algorithm for modeling the flow in real turbine stages by means of
such test cases is not easy; it has been discussed in Refs. 3,4.
The comparison of the numerical modeling and measurements
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
of the flow in real turbine of the large output is not often published. The crucial value of such verification comes down to the
assessment of flow conditions for real LP steam turbine stages by
means of the CFD code.
This paper presents the results of experimental and numerical
investigations of 3D flows through the last two stages of the LP
part of a 360 MW turbine. The problems of condensation in a LP
steam turbine has been presented and discussed also, e.g., by
Hesler et al. 5, Bakhtar and Heaton 6, and Petr et al. 7.
Measurements of thermodynamic parameters, the size, and concentration of water droplets together with CFD results have been
described below. The presented work was performed by research
teams of the Institute of Fluid-Flow Machinery from Gdansk PL,
the Czech Technical University CTU of Prague, and the Silesian
University of Technology PL within the framework of the Research Project No. 3T10B04726 sponsored by the State Committee for Scientific Research in Poland.
4
2 drhom
= l r3Jhom + 3nhomrhom
3
dt
rhet = r3p +
3 y het
4 lnhet
1/3
h = hv1 y + h1y
1/2
mv3/2
v2
4r2
exp
3kTv
l
C= 1+2
1 hv hl hv hl 1
2
+ 1 RTv
RTv
2
lfpv fps pv ps
where
4
= v/1 y
The value of the nonequilibrium wetness fraction y is calculated
as the sum of the homogeneous and heterogeneous wetness obtained from the conservation equations for the liquid phase 1 and
2.
041301-2 / Vol. 131, APRIL 2009
p
BT
= ZT, = AT +
RT
s = sv1 y + s1y
b
b
fp = ln p + c ln
2
2
c
b
c
1
b
1+
b = ATRT
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. 1 Location of probes in the blade system of a 360 MW turbine LP part of Belchatow
Power Plant: meridional section a, cross section b, and photo of measurement stand
c
c = ATRT2 + 4pBTRT1/2
Equation 7 has an ideal gas form since A 1 and B 0 .
The further behavior of the critical droplets can be described by
suitable droplet growth law. The size of the droplets for the vapor
under low pressure is much smaller than mean free path of the
vapor molecules. Therefore, the growth of the droplets should be
governed by considering molecular and macroscopic transport
processes HertzKnudsen model. Difficulties to choose the condensation and accommodation coefficients make the application
of the HertzKnudsen model very difficult to calculate. This problem can be avoided by using Gyarmathys droplet growth model,
which takes into account diffusion of vapor molecules through the
surrounding vapor, as well as heat and mass transfer, and the
influence of the capillarity:
r r Ts Tv
v
dr 1
2
=
dt l 1 + 3.18 Kn
r hv hl
1000ns
nwM ww
aw = 1
aM,
i
i
l
M l M l max
i=1
aw = 1
aM
i
i
l max,
i=1
M l M l max
10
= oT + BM l
11
Numerical 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
Fig. 2 Plate probe for thermal and flow measurements a including the optical system b
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
bly be explained by the heterogeneous nature of the steam condensation process at the presence of other chemical compounds.
The presence of chemical inclusions was confirmed by pH conductivity tests and chemical analyses of the first condensate and
condenser water samples. This is presented in Table 1. Noteworthy are high numbers of Na and Cl ions, which, according to
investigations by Steltz et al. 19 and Svoboda et al. 20, accelerate the condensation process.
area, detected by the optical system of the probe, is situated behind the third stage, below the saturation line, which could possi kg s1 case 1
Table 1 Calculated mass flow rates m
Adiabatic
Homogeneous
condensation
Homo-/heterogeneous
condensation
Heterogeneous
condensation
103.5
102.2
102.1
101.8
4.2 Measurement of Liquid Phase Characteristics. Measuring the liquid-phase structure in the LP part was a special part of
360 MW turbine tests. The optical measurements require a special
measuring device. Therefore, the CTU team was invited, who has
long-time experience with this problem. Measurements of the
droplet size distribution and steam wetness were carried out by an
extinction probe, developed at CTU 7. The probe measures the
multispectral extinction of light, i.e., attenuation I / Io, caused by
the light scattering on droplets, at the series of light wavelengths.
The universal measuring system with an extinction probe was
used for measuring the multispectral light extinction. The probe
measures the transmittance I / Ioi, caused by the light scattering
on droplets, of the light beam transmitted through the wet steam at
the series of light wavelengths ii = 1 , . . .. The layout of the measuring system is shown in Fig. 5a, while the measuring head of
the extinction probe with the 25 mm diameter is presented in Fig.
5b.
The distribution function FD of the droplet sizes, together
with the wetness fraction y, is then found by means of computer
analysis using Mies theory of light scattering on droplets of the
measured transmittance data I / Io. In principal, the problem leads
to the solution of the system of the first kind Fredholm integral
equations.
The structure of fine droplets at the measurement site was
evaluated from the measurement data I / Io by means of the CTU
inversion method using the well-known exponential law:
1
Io
ln
l
I
Dmax
ED/ i,mFDD2dD
12
Dmin
Fig. 5 Measuring system of CTU extinction probe a and head of extinction probe b
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
Numerical Results
The presented experimental data were the basis for the CFD
calculations and the validation tests. A numerical analysis was
performed for the geometry of the last two stages of the five-stage
LP cylinder of a 360 MW steam turbine 3000 rpm. The analyzed
last two stages operate in the conditions under saturation line.
Two different measurements were performed in the power plant
for the loads 300 MW case 1 and 360 MW case 2. In case 1,
only flow parameters were measured. In case 2, the liquid phase
parameters were included in the experimental study.
Fig. 6 Variation in transmittance measured at the blade midspan behind the L-0, L-1, and L-2 stages
D32 =
Dmax
FDD3dD
Dmin
13
Dmax
FDD2dD
Dmin
n=
Dmax
14
FDdD
Dmin
l
1+
6 v
1
Dmax
3
FDD dD
Dmin
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
Fig. 7 Static pressure and flow angle distributions: a between the last two stages and b at the outlet
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 Wetness fraction at the midspan section: a homogeneous, b homoheterogeneous, and c heterogeneous condensations
Fig. 9 Wetness fraction distributions: a between the last two stages and b at the outlet
Table 2 Results of chemical analyses of condensate samples
Notation
Unit
Condensate behind the LP part
First condensate
pH
Na+
K+
NH+4
Cl
SO4
SiO2
Fe
88.4
6.86.9
g / kg
23
180
g / kg
5
40
g / kg
60
g / kg
2
130
g / kg
12
90
g / kg
23
900
g / kg
78
1200
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 Flow angle and static pressure distribution: a between the last two stages and b at the outlet
Fig. 11 Wetness fraction contours at midspan section: global wetness fraction contours homogeneous and heterogeneous a and wetness fraction
contours generated only due to the heterogeneous condensation b
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. 12 Wetness fraction and mean droplet Sauter diameter: a between the last two stages and b at the outlet
are in the rotor row of the first stage as can be deduced from the
wetness contours in Figs. 11a and 11b. The contours are shown
for the global wetness fraction obtained due to the homogeneous
and heterogeneous condensation processes and for the wetness
fraction generated only on the soluble particles.
The spanwise distributions of the wetness fractions and droplet
radii in the gap between stages and at the stage outlet are shown
together with experimental data in Fig. 12. The relatively good
agreement in the wetness fraction distributions is observed below
the midsection at the outlet. The calculated wetness fraction between the stages is about 30% smaller than that estimated from
the experimental data. This difference is high in the midpart of the
blade and decreases in the tip sections. The discrepancies at the
outlet are smaller than those between the stages. This can be explained by fact that at the outlet the equilibrium in the wet steam
was achieved. The simplifications in condensation models and the
accuracy of the real-scale measurements bear on the discrepancies
in the flow parameters between the calculations and experiment.
In Fig. 12, the spanwise distributions of the Sauter mean diameters are presented. From calculations, the volume averaged diameter is obtained; therefore, the Sauter diameter was calculated on
the basis of droplet size spectrum known from experiment.
The calculated droplet radii show a more even distribution and
are confined within the same range in the both sections.
041301-10 / Vol. 131, APRIL 2009
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
Acknowledgment
The authors would like to thank the Polish Ministry of Science
and Information Technology for the financial support of the Research Project No. PB 1431/T10/2004/26 and the management of
the Belchatow Power Plant for the consent for onsite measurements.
Nomenclature
aw
D
D32
E
F
h
Io
I
J
k
l
L
m
m
Ml
Mw
n
p
R
r
r
s
T
u,v,w
y
Z
x,y,z
0
hom
het
v
total parameter
homogeneous
heterogeneous
vapor phase
Subscripts
l liquid phase
s saturation value
w water
References
1 Keller, H., 1974, Erosionkorrosion an Nassdampfturbinen, VGB Kraftwerktechnik, 54, pp. 292295.
2 Jonas, O., 1985, Steam Turbine Corrosion, Mat. Perf., 242, pp. 918.
3 Bakhtar, F., White, A. J., and Mashmoushy, H., 2005, Theoretical Treatments
of Two-Dimensional, Two-Phase Flows of Steam and Comparison With Cascade Measurements, Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. Sci., 219,
pp. 13351355.
4 Wrblewski, W., Dykas, S., and Gepert-Niedobecka, A., 2007, Condensing
Steam Flow ComputationsValidation Problems, Seventh European Conference on Turbomachinery, Conference Proceedings, Athens, K. D. Papailiou, F.
Martelli, and M. Manna, eds., pp. 841848.
5 Hesler, S., Shuster, L., and McCloskey, T., 1998, Optical Probe of Measurement of Steam Wetness Friction in LP Turbine, IPJC II, ASME, New York,
PWR-Vol. 33.
6 Bakhtar, F., and Heaton, A. V., 2005, Effects of Wake Chopping on Droplet
Sizes in Steam Turbines, J. Mech. Eng. Sci., 219C12, pp. 13571367.
7 Petr, V., Kolovratnk, M., and Hanzal, V., 2003, Instrumentation and Tests on
Droplet Nucleation in LP Steam Turbines, Power Plant Chemistry, 57, pp.
389395.
8 White, A. J., Young, J. B., and Walters, P. T., 1996, Experimental Validation
of Condensing Flow Theory for a Stationary Cascade of Steam Turbine
Blade, Philos. Trans. R. Soc. London, Ser. A, 354, pp. 5988.
9 Wagner, W., Cooper, J. R., Dittmann, A., Kijima, J., Kretzschmarr, H.-J.,
Kruse, A., Mare, R., Oguchi, K., Sato, H., Stcker, I., ifner, O., Takaishi, Y.,
Tanishita, I., Trbenbach, J., and Willkommen, Th., 2000, The IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and
Steam, Trans. ASME: J. Eng. Gas Turbines Power, 1221, pp. 150181.
10 Frenkel, J., 1946, Kinetic Theory of Liquids, Oxford University Press, New
York.
11 Gyarmathy, G., 1960, Grundlagen einer Theorie der Nassdampfturbine,
Ph.D. thesis, ETH, Zrich.
12 Kantrowitz, A., 1951, Nucleation in Very Rapid Vapor Expansions, J. Chem.
Phys., 19, pp. 10971100.
13 Gorbunov, B., and Hamilton, R., 1997, Water Nucleation on Aerosol Particles
Containing Both Soluble and Insoluble Substances, J. Aerosol Sci., 282, pp.
239248.
14 Dykas, S., 2001, Numerical Calculation of the Steam Condensing Flow,
Task Quarterly, Scientific Bulletin of Academic Computer Centre in Gdansk,
54, pp. 495519.
15 Wrblewski, W., Dykas, S., and Gepert, A., 2006, Modeling of Steam Flow
With Heterogeneous Condensation, Monografia 95, Wydawnictwo Politechniki Slskiej, Gliwice, in Polish.
16 Lagun, V. P., Simoyu, L. L., and Frumin, Yu. Z., 1975, Full-Scale Tests of the
Exhaust Hoods of a High-Capacity Steam Turbine, Teploenergetika Moscow,
Russ. Fed., 222, pp. 3135.
17 Gardzilewicz, A., and Marcinkowski, S., 1995, Diagnosis of LP Steam Turbine Prospects of Measuring Techniques, Proceedings of the ASME Joint
Power Generation Conference, Vol. 3, pp. 349358.
18 Gardzilewicz, A., Marcinkowski, S., Karwacki, J., and Kurant, B., 2006, Results of Measurements of Condensing Steam Flow in LP Part of Turbine 18 K
in Belchatow Power Plant, IF-FM Report Nos. 5681 and 6157.
19 Steltz, W., Lindsay, P., and Lee, W., 1983, The Verification of Concentrated
Impurities in Low Pressure Steam Turbine, ASME Winter Meeting, WA.
20 Svododa, R., Phlug, H.-D., and Warnecke, T., 2003,Investigation into the
Composition of the Early Condensate in Steam Turbines, Power Plant Chemistry GmbH, Neulusshei.
21 Petr, V., and Kolovratnik, M., 2001, Heterogeneous Effects in the Droplet
Nucleation Process in LP Steam Turbines, Proceedings of the Fourth European Conference on Turbomachinery, Florence, Italy, pp. 783792.
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. Kawaji1
Department of Chemical Engineering and Applied
Chemistry,
University of Toronto,
200 College Street,
Toronto, ON, M5S 3E5, Canada
K. Mori
Department of Mechanical Engineering,
Osaka Electro-Communication University,
18-8 Hatsu-cho,
Neyagawa, Osaka, Japan 572-8530
D. Bolintineanu
Department of Chemical Engineering and Applied
Chemistry,
University of Toronto,
200 College Street,
Toronto, ON, M5S 3E5, Canada
Introduction
cated that there are differences in the flow patterns between microchannels of D 100 m and minichannels of D 1 mm.
Kawaji and co-workers 7,8 investigated the two-phase flow of
nitrogen gas and water through horizontal circular microchannels
of 50 250 m diameter. The gas and liquid were mixed in a
reducer inlet where the diameter changed from 750 m down to
the microchannel diameter. They reported significant differences
in the flow pattern maps and void fraction data from those previously described for minichannels with diameters above 1 mm.
In particular, they reported the existence of flow patterns unique to
microchannels liquid-ring flow and serpentinelike gas core flow.
Furthermore, the time-average void fraction data deviated from
Armands correlation 3, which is normally applicable to minichannels with diameters 1 mm. This was in contrast with the
void fraction data of Serizawa et al. 6 for a 20 m circular
microchannel, which could be well described by Armands type
correlation.
Chung et al. 9 later found that the two-phase flow characteristics in a 96 m square channel and a 100 m circular channel
connected to the same reducer inlet used by Kawahara et al. 7
were quite similar and independent of the channel cross section
for a mixture of nitrogen gas and water. Furthermore, Kawahara et
al. 10 found only a moderate effect of fluid properties such as
liquid viscosity and surface tension on the two-phase flow characteristics in a microchannel.
Besides the reducer-type inlet, several different inlet designs
and gas-liquid mixing methods have been used to inject a twophase mixture into a microchannel. A simple T-junction has also
been popular particularly for microfluidic applications as used by
Thorsen et al. 11, Okushima et al. 12, Gnther et al. 13, Xu et
al. 14, Garstecki et al. 15, and Ide et al. 16, among others. In
many cases, the T-junction diameter was the same as or close to
the microchannel diameter, and its use was intended to produce
monodisperse microbubbles. On the other hand, a much larger
1/8 in. diameter T-junction was used by Xiong and Chung 17
upstream of a rectangular microchannel. The gas was injected into
the liquid stream through a 200 m diameter needle, which was
inserted from the branch of the T-junction so that the needle tip
would be positioned in the middle of the T-junction.
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 recent experiments, Yue et al. 18 used a Y-junction connected to square microchannels with 200 m and 400 m hydraulic diameters to absorb CO2 gas into a water stream. Another
inlet geometry used in the past involved a flow focusing method
in a cross junction by Gan-Calvo and Gordillo 19, Cubaud
and Ho 20, and Garstecki et al. 21, where the gas flowed
straight through the main line into the microchannel and liquid
was injected from the sides through two branches at equal flow
rates. This way, the gas stream could be pinched off to generate
discrete small gas bubbles, short plugs, or long slugs depending
on the gas and liquid injection rates. Another variation of the
T-junction inlet was used by Waelchli and von Rohr 22, who
injected gas and liquid streams from the opposite ends of the main
channel. The two phases mixed at the T-junction where several
50 m diameter round pins were placed to serve as a static mixer,
and then the mixture flowed into the branch line connected to a
microchannel.
As reviewed above, different types of inlet and gas-liquid mixing methods have been used in adiabatic two-phase flow experiments in microchannels, but no systematic comparisons of the
two-phase flow characteristics have been carried out. Some of the
results in the literature have shown significant differences, but the
diameters of the microchannels were not the same, so the effects
of the inlet geometry and mixing method cannot be readily isolated from those of the microchannel diameter. Hence, the objective of this work was to investigate the effect of the inlet geometry
using the same microchannel of 100 m diameter and identify
the key reasons for differences in the two-phase flow characteristics observed in the microchannel. The same experimental apparatus and approach used by Kawahara et al. 7 and Chung and
Kawaji 8 have been used again with the same microchannel but
different inlet designs and gas-liquid mixing configurations.
Since there is no phase change heat transfer involved in this
work, the results are not directly applicable to heat exchangers.
Most relevant applications would be in microreactors in which gas
and liquid streams are contacted in a microchannel to effect gas
absorption with or without a chemical reaction. In the design of
such microreactors, the gas-liquid inlet and mixing method can be
chosen from a T-junction, Y-junction, co-axial reducer, cross junction, and others. Thus, it would be useful to understand better how
the two-phase flow characteristics are affected by the inlet geometry and gas-liquid mixing method.
Fig. 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
Uncertainty range
0.1% to 1%
1% to 10%
2.1% to 3%
3% to 12%
0.05
1% to 4%
2% to 7%
2% to 7%
fects. The window for flow visualization was located as far back
as possible 46.1 mm from the gas-liquid mixing point on the
microchannel to achieve a well-developed flow.
Three pressure transducers PX303 series from Omega Engineering
Stamford,
CT
covering
different
ranges
3450 8.6 kPa, 690 1.7 kPa, and 210 0.52 kPa were used
to measure the pressure drop between the microchannel inlet and
outlet. The liquid was discharged freely from the test section, so
that the pressure was atmospheric and temperature ambient at the
microchannel outlet.
The liquid flow rate was determined by using a load cell
50.97 0.255 g, model ELJ-0.5N from Entran Devices Fairfield, NJ to weigh the discharged liquid collected over time. The
gas flow rate was measured using three mass flow meters model
179A-23012 from MKS Instruments Andover, MA, and models
8172-0411 and 8172-0451 from Matheson Gas Products Montgomeryville, PA to encompass a wide range of conditions
1 0.01 sccm, 10 0.1 sccm, and 50 0.5 sccm. The analog
signals for all of the pressure, temperature, and mass flow rate
readings were sampled at 10 Hz by a 16-bit data acquisition system model NI 6035E from National Instruments Austin, TX.
Background illumination to see the flow patterns was generated
by a cold lamp model DCR III from Schott-Fostec Auburn,
NY located behind the test section. The view field was magnified by attaching a 20 microscope objective lens model 378802-2 from Mitsutoyo America Corporation Aurora, IL to the
video camera. A monochrome CCD camera model TM-1040
from Pulnix America Sunnyvale, CA with a resolution of
1024H 1024V pixels, a maximum frame rate of 30 fps, and
a shutter speed of 1/16,000 s was used for image capturing. The
problem of optical distortion was dealt with earlier by Kawahara
et al. 7. They elected to not correct for optical distortion so that
the images of the gas-liquid interfacial structure near the channel
wall could be enlarged and better examined. Likewise, optical
distortion is not addressed here for the same reason.
The previous experiments covered a range of 0.01 jL
5.77 m / s for the superficial liquid velocity and 0.02 jG,avg
73 m / s for the average superficial gas velocity. The new experiments covered slightly smaller ranges: 0.027 jL 2.05 m / s
and 0.1 jG,avg 50 m / s. The values of dimensionless parameters calculated for the previous experiments are the Bond number
0.000083 Bo 0.0095, the superficial Weber number for liquid 0.00057 WeLS 116 and gas 0.0000020 WeGS 12,
the superficial Reynolds number for liquid 2 ReLS 2300 and
gas 0.4 ReGS 751, and the capillary number 0.00013 CaL
0.070. These values indicate the relative magnitudes of the
forces
involved
inertia surface tension viscous force
gravity.
2.1 Measurement Uncertainties and Data Analyses. The
measurement uncertainties estimated for the present experiments
are summarized in Table 1. A propagation of error analysis 23
was performed to estimate the uncertainties in the derived quantities such as superficial velocities, void fraction, and pressure
gradients.
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. 5 Two-phase flow pattern maps for a micro-T-junction inlet a ConFig. 1 and b ConFig. 2, and c a reducer inlet used
by Kawahara et al. 7 an open circle indicates the occasional
appearance of a ring-film flow pattern
3 Slug flow: the presence of gas slugs that are longer than the
viewing window
4 Semi-annular flow: long gas core surrounded by a liquid
film and interrupted by infrequent liquid slugs
5 Ring flow: 10% or more of the images for a given run have
a ring-shaped liquid film surrounding the gas core
The two-phase flow pattern maps were then constructed based
on the above definitions for the two T-junction inlet configurations
as shown in Fig. 5. At low superficial gas velocities jG 1 m / s,
bubbly flow, plug flow, and slug flow patterns appeared at progressively lower liquid flow rates. At moderate to high superficial
gas velocities jG 2 m / s, slug and semi-annular flow patterns
were observed. An open circle in Figs. 5a and 5b indicates
the occasional appearance of a ring-film flow pattern.
The difference in the gas and liquid injection configurations for
the micro-T-junction inlet had only a small effect on the flow
pattern map as clear from a comparison of Fig. 5a and 5b.
However, these results for the micro-T-junction inlet were significantly different from the map previously obtained for a reducer
inlet case by Kawahara et al. 7 and shown in Fig. 5c. When the
inlet section consisted of a series of reducer sections from
750 m to 100 m Fig. 2, the predominant flow pattern ob041302-4 / Vol. 131, APRIL 2009
served was a slug flow with a long gas core surrounded by a thin,
thick or ring-shaped film. Depending on the frequency of the appearance of different liquid film shapes, the flow patterns were
defined as slug-ring, ring-slug, and multiple flows Fig. 5c.
Although a semi-annular flow pattern was observed in both
micro-T-junction and reducer inlet experiments, a strong difference existed between them in the appearance of short gas bubbles
and plugs. This is due to the fact that small gas bubbles are readily
created at the T-junction inlet as reported in previous studies. The
mechanism of small gas bubble generation at the micro-T-junction
can be readily illustrated as shown in Fig. 6. Such a mechanism
had been previously investigated for the ConFig. 1 case liquid
injected into the main channel while gas is injected into the
branch, Fig. 6a by Shibata et al. 24. Short gas plugs were
observed to be generated as a result of interfacial instability or
pinch-off at the corner of the T-junction as shown in Figs.
6b6e.
Although no flow visualization was conducted for ConFig. 2, a
similar interfacial instability or pinching off of a gas bubble or
plug at the micro-T-junction can be conjectured. Thus, regardless
of which phase is injected into the main channel or the branch, the
use of a T-junction to mix and inject the gas-liquid mixture into a
microchannel can result in two-phase flow patterns that are very
different from those observed in the reducer inlet case. The larger
inlet cross section D = 750 m for the reducer inlet used by
Kawahara et al. 7 would allow larger gas bubbles to form upstream of the microchannel test section. When these gas bubbles
flow into the microchannel with a much smaller diameter D
= 100 m, the gas slug length would increase by more than 50
times, as illustrated in Fig. 7. Thus, for example, if the diameter of
the gas bubble formed in the inlet section were 2 mm, then it
would form a 100 mm long gas slug in the microchannel. This can
explain why Kawahara et al. 7 obtained mostly long gas core
flows with a reducer inlet rather than short gas plug flows.
The present flow pattern results obtained with a T-junction inlet
are now compared with the previous results obtained with other
types of inlet as shown in Fig. 8. Only the bubbly-to-plug or slug
and slug-to-annular or semi-annular flow pattern transition
boundaries are indicated. Although the transition boundaries are
represented by thin lines, the flow patterns did not change as
sharply as suggested by those lines. Nevertheless, the present
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. 7 Formation of a long gas slug from a bubble in a reducer inlet not to scale
void fraction values from about 300 images. The effect of optical
distortion on the above-mentioned method of measuring the void
fraction, particularly for long gas core flows with ring-shaped
films and thick liquid films, was previously addressed by Kawahara et al. 7.
The time-averaged void fraction data for the T-junction inlet
are plotted against the volumetric quality =JG / JG + JL in Figs.
9 and 10 for Configurations 1 and 2, respectively. The homogeneous flow model = and the Armand correlation 3 for narrow channels recommended by Ali et al. 25 are also shown in
the figures, respectively. Ali et al. 25 recommended using a correlation of the form = 0.8 originally developed by Armand 3
for narrow rectangular channels with DH 1 mm.
Both Configurations 1 and 2 data show a linear relationship
between the void fraction and volumetric quality, consistent with
the homogeneous flow model and Armands correlation, but unlike the nonlinear relationship obtained by Kawahara et al. 7 for
the reducer inlet case and given by
=
0.030.5
1 0.970.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
short gas and liquid plugs as shown in Fig. 4 would lead to a slip
ratio close to unity, which is representative of a homogeneous
flow. Since such a plug flow pattern was not observed in the
reducer inlet case, the void fraction deviated strongly from the
homogeneous flow model in the experiments of Kawahara et al.
7. The present data show that the inlet geometry and the method
of gas-liquid mixing can drastically change the time-average void
fraction for the same microchannel at given superficial gas and
liquid velocities.
3.3 Friction Pressure Drop. The two-phase friction pressure
drop was evaluated from the measured total pressure drop by subtracting the acceleration pressure drop. The data were correlated
using the LockhartMartinelli correlation 26, which was developed based on a separated flow assumption. This correlation uses
a two-phase friction multiplier L2 defined as follows:
P f
Z
= L2
TP
P f
Z
2
L
P f /ZL
P f /ZG
L2 = 1 +
C 1
+
X X2
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
References
Conclusion
Acknowledgment
This work was supported by a Discovery research grant from
the Natural Sciences and Engineering Research Council of
Canada. K. Mori would like to thank Osaka ElectroCommunication University for a research leave grant, which allowed a one-year stay at the University of Toronto.
Nomenclature
Bo
C
Ca
D
g
j
L
m
P
Re
u
We
X
Z
Greek Symbols
P
L2
Subscripts
f
G
L
M
S superficial
TP two-phase
Bond number
coefficient in Eq. 4
capillary number
tube inner diameter m
gravitational acceleration m / s2
superficial velocity m/s
length m
mass flow rate kg/s
pressure Pa
Reynolds number
axial velocity m/s
Weber number
LockhartMartinelli parameter
flow direction coordinate m
volumetric quality
void fraction
pressure drop
dynamic viscosity Pa s
two-phase friction multiplier
density kg/ m3
surface tension N/m
friction
gas phase
liquid phase
mixing section
1 Serizawa, A., 2006, Gas-Liquid Two-phase Flow in Microchannels, Multiphase Flow Handbook, C. T. Crowe, ed., CRC Press, Boca Raton, FL.
2 Thome, J. R., 2006, State-of-the-Art Overview of Boiling and Two-Phase
Flows in Microchannels, Heat Transfer Eng., 279, pp. 419.
3 Armand, A. A., 1946, The Resistance during the Movement of a Two-phase
System in Horizontal Pipes, Izv. Vses. Teplotekh. Inst. AERE-Lib/Trans.
828, 1, pp. 1623.
4 Stanley, R. S., Barron, R. F., and Ameel, T. A., 1997, Two-Phase Flow in
Microchannels, Micro-Electro-Mechanical Systems (MEMS), ASME, DSCVol. 62/HTD-Vol. 354, New York, pp. 143152.
5 Serizawa, A., and Feng, Z. P., 2001, Two-Phase Flow in Micro-Channels,
Proceedings of the Fourth International Conference on Multiphase Flow, New
Orleans, LA, May 27Jun. 1.
6 Serizawa, A., Feng, Z., and Kawara, Z., 2002, Two-Phase Flow in Microchannels, Exp. Therm. Fluid Sci., 2667, pp. 703714.
7 Kawahara, A., Chung, P. M.-Y., and Kawaji, M., 2002, Investigation of TwoPhase Flow Pattern, Void Fraction and Pressure Drop in a Microchannel, Int.
J. Multiphase Flow, 289, pp. 14111435.
8 Chung, P. M. Y., and Kawaji, M., 2004, The Effect of Channel Diameter on
Adiabatic Two-Phase Flow Characteristics in Microchannels, Int. J. Multiphase Flow, 3078, pp. 735761.
9 Chung, P. M.-Y., Kawaji, M., Kawahara, A., and Shibata, Y., 2004, TwoPhase Flow Through Square and Circular MicrochannelsEffect of Channel
Geometry, ASME J. Fluids Eng., 1264, pp. 546552.
10 Kawahara, A., Kawaji, M., Chung, P. M.-Y., Sadatomi, M., and Okayama, K.,
2005, Effects of Channel Diameter and Liquid Property on Void Fraction in
Adiabatic Two-Phase Flow Through Microchannels, Heat Transfer Eng.,
263, pp. 1319.
11 Thorsen, T., Roberts, R. W., Arnold, F. H., and Quake, S. R., 2001, Dynamic
Pattern Formation in a Vesicle-Generating Microfluidic Device, Phys. Rev.
Lett., 86, pp. 41634166.
12 Okushima, S., Nishisako, T., Torii, T., and Higuchi, T., 2004, Controlled
Production of Monodisperse Double Emulsions by Two-Step Droplet Breakup
in Microfluidic Devices, Langmuir, 20, pp. 99059908.
13 Gnther, A., Khan, S. A., Thalmann, M., Trachsel, F., and Jensen, K. F., 2004,
Transport and Reaction in Microscale Segmented GasLiquid Flow, Lab
Chip, 4, pp. 278286.
14 Xu, J. H., Li, S. W., Wang, Y. J., and Luo, G. S., 2006, Controllable GasLiquid Two Phase Flow Patterns and Monodispersive Microbubbles in a Microfluidic T-Junction Device, Appl. Phys. Lett., 88, pp. 133506.
15 Garstecki, P., Fuerstman, M. J., Stone, H. A., and Whitesides, G. M., 2006,
Formation of Droplets and Bubbles in a Microfluidic T-JunctionScaling
and Mechanism of Break-Up, Lab Chip, 6, pp. 437446.
16 Ide, H., Kimura, R., and Kawaji, M., 2007, Optical Measurement of Void
Fraction and Bubble Size Distributions in a Microchannel, Heat Transfer
Eng., 288, pp. 713719.
17 Xiong, R., and Chung, J. N., 2007, An Experimental Study of the Size Effect
on Adiabatic Gas-Liquid Two-phase Flow Patterns and Void Fraction in Microchannels, Phys. Fluids, 19, p. 033301.
18 Yue, J., Chen, G. W., and Yuan, Q., 2004, Pressure Drops of Single and
Two-Phase Flows Through T-Type Microchannel Mixers, Chem. Eng. J.,
102, pp. 1124.
19 Gan-Calvo, A. M., and Gordillo, J. M., 2001, Perfectly Monodisperse Microbubbling by Capillary Flow Focusing, Phys. Rev. Lett., 87, p. 274501.
20 Cubaud, T., and Ho, C. M., 2004, Transport of Bubbles in Square Microchannels, Phys. Fluids, 16, pp. 45754585.
21 Garstecki, P., Stone, H. A., and Whitesides, G. M., 2005, Mechanism for
Flow-Rate Controlled Breakup in Confined Geometries: A Route to Monodisperse Emulsions, Phys. Rev. Lett., 94, p. 164501.
22 Waelchli, S., and von Rohr, P. R., 2006, Two-Phase Flow Characteristics in
GasLiquid Microreactors, Int. J. Multiph. Flow, 327, pp. 791806.
23 Kline, S. J., and McClintock, F. A., 1953, Describing Uncertainties in Single
Sample Experiments, Mech. Eng. Am. Soc. Mech. Eng., 75, pp. 38.
24 Shibata, Y., Ikeda, K., Hanri, T., Chung, P. M. Y., and Kawaji, M., 2003,
Study on Two-Phase Flow Through a T-Junction, Proceedings of the Fourth
ASME-JSME Joint Fluids Engineering Conference, Honolulu, Hawaii, Jul.
611, Paper No. FDSM2003-45385.
25 Ali, M. I., Sadatomi, M., and Kawaji, M., 1993, Two-Phase Flow in Narrow
Channels Between Two Flat Plates, Can. J. Chem. Eng., 715, pp. 657666.
26 Lockhart, R. W., and Martinelli, R. C., 1949, Proposed Correlation of Data
for Isothermal Two-Phase, Two-Component Flow in Pipes, Chem. Eng. Prog.,
451, pp. 3948.
27 Chisholm, D., and Laird, A. D. K., 1958, Two-Phase Flow in Rough Tubes,
Trans. ASME, 802, pp. 276286.
28 Mishima, K., and Hibiki, T., 1996, Some Characteristics of Air-Water TwoPhase Flow in Small Diameter Vertical Tubes, Int. J. Multiphase Flow, 224,
pp. 703712.
29 Kawaji, M., 1999, Fluid Mechanics Aspects of Two-phase Flow, Handbook
of Phase Change: Boiling and Condensation, Taylor & Francis, London, Chap.
9, pp. 205259.
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
Thien X. Dinh
e-mail: thien@cfd.ritsumei.ac.jp
Yoshifumi Ogami
Department of Mechanical Engineering,
Ritsumeikan University,
1-1-1 Nojihigashi,
Kusatsu, Shiga 525-8577, Japan
Introduction
structure of the device is compatible to planar layer-by-layer geometry that is feasibly produced by lithography-based microfabrication techniques. The working principle and some primitive
results are presented. All velocity components are examined.
However, the axial component is mostly studied, since it is the
most important flow characteristic.
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
t
u
+ u u = p + u
t
c pT
+ u c pT = T
t
1
fN
r,t = t 1
r
R
2 2
Because 0 of the order of tens of microns is small as compared to the depth of the pump chamber 300 m, it is not
necessary to use the deformation mesh in the simulation. The local
rate of deformation of the diaphragm is converted to local velocity
applied on the inlet surface by taking derivative of Eq. 5 respective with time as
vr,t = 2 f0 1
r
R
cos2 ft
In order to facilitate the discussion, a Cartesian coordinate system is designated with the origin located at the center of the
nozzle outlet. The x-, y-, and z-axes correspond to the axial, horizontal, and vertical directions, respectively, as shown in Fig. 1.
The corresponding velocity components along the x-, y-, and
z-directions are denoted as u, v, and w.
Figures 3a and 3b shows the instantaneous vector velocity
field at the throat for the discharge and suction phases, respecTransactions 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 Flow behavior at the throat in the discharge phase a and suction phase b
tively. The vectors are plotted with constant length for simplification. The magnitude of the vector velocity is represented by the
color bar. The highest and lowest values are represented by red
and blue, respectively. The working principle of the device is
clearly observed. In the discharge phase, there is a ring-shaped
circulatory flow its cross sections are indicated by characters A
and B in the front and side views, as shown in Fig. 3a around
the throat. It intercepts the flow from the throat to the backchannels. The flow from the housing chamber goes straight to the
drive-channel. In the suction phase, the flows from the backchannels are comparable to that from the drive-channel at the
throat. Then, the former slows down the fluid from the drivechannel to the housing chamber. Moreover, the secondary flow
from the wall of the drive-channel indicated by the arrows in Fig.
3b also contributes to slowing down the flow from this channel
to the housing chamber. Consequently, only a small portion of the
flow in the drive-channel, which was created in the pump phase, is
sucked back into the housing chamber. The other major portion
still flows in the drive-channel. We can observe this effect by the
division of the vector velocity field into opposite directions at the
throat, as shown in Fig. 3b.
Figure 4 shows the variation in the spatial-average axial velocity at the nozzle outlet, UJ, against the dimensionless time t f.
We observe that the flow attains the steady state after 4 cycles of
vibration of the membrane. It is also observed that UJ completes 1
cycle during the dimensionless time t f = 1. In other words, the
flow in the main chamber oscillates at the same frequency as the
PZT membrane. It should be noted that the PZT vibrates at a very
high frequency. Therefore, it is expected that the flow is almost in
the steady state for the applications time scales.
Fig. 4 Variation in the spatial-average axial velocity at the outlet of the nozzle against time. The maximum amplitude of the
PZT membrane is 10 m
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 Axial velocity of the jet against cross-distance for different amplitude vibrations of the membrane and different axial
distance, x. The left-hand side of each figure shows a plot of
ux , 0 , z / ucx against z / z1/2x and the right-hand side shows
the variation in ux , y , 0 / ucx against y / y1/2x.
Fig. 7 Axial velocity of the jet against cross-distance for different amplitude vibrations of the membrane at x = 1.5 mm. The
left-hand side of each figure shows a plot of ux , 0 , z / ucx
against z / z1/2x and the right-hand side shows the variation in
ux , y , 0 / ucx against y / y1/2x. The approximate curve is plotted by Eq. 10.
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
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
The two remaining components of the flow are small as compared to the axial component.
Furthermore, the characteristics of the centerline axial velocity,
as well as the half-widths of the jet, should be addressed in a
general way. The device has been investigated only for one specific dimension. Since the geometry, throat, and nozzle may contribute significantly in optimizing the operation of the device,
these factors should be carefully considered in the future.
References
1 Woolley, A. T., Hadley, D., Landre, P., deMello, A. J., Mathies, R. A., and
Northrup, M. A., 1996, Functional Integration of PCR Amplification and
Capillary Electrophoresis in a Microfabricated DNA Analysis Device, Anal.
Chem., 68, pp. 40814087.
2 Taylor, M. T., Nguyen, P., Ching, J., and Petersen, K. E., 2003, Simulation of
4
5
6
7
8
Microfluidic Pumping in a Genomic DNA Blood-Processing Cassette, J. Micromech. Microeng., 13, pp. 201208.
Zhang, L., Koo, J. M., Jiang, L., Asheghi, M., Goodson, K. E., Santiago, J. G.,
and Kenny, T. W., 2002, Measurements and Modeling of Two-Phase Flow in
Microchannels With Nearly Constant Heat Flux Boundary Conditions, J. Microelectromech. Syst., 11, pp. 1219.
Dau, V. T., Dao, V. D., Shiozawa, T., Kumagai, H., and Sugiyama, S., 2005,
A Dual Axis Gas Gyroscope Utilizing Low Doped Silicon Thermistor, 18th
IEEE MEMS, pp. 626629.
Zhou, J., Yan, G., Zhu, Y., Xiao, Z., and Fan, J., 2005, Design and Fabrication
of a Microfluidic Angular Rate Sensor, 18th IEEE International Conference
on Micro Electro Mechanical Systems MEMS, pp. 363366.
Saito, K., Komatsu, G., and Kato, K., 2001, An Analysis of Fluid Behavior on
Gas Angular Rate Sensor, JSME Int. J., Ser. C, 67, pp. 135140.
Bassi, A., Gennar, F., and Symonds, P. S., 2003, Anomalous Elastic-Plastic
Respond to Short Pulse Loading of Circular Plates, Int. J. Impact Eng., 28,
pp. 6591.
Pope, S. B., 2000, Turbulent Flows, Cambridge University, Cambridge.
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
Aijie Han
Yu Qiao1
e-mail: yqiao@ucsd.edu
Department of Structural Engineering,
University of California, San Diego,
La Jolla, CA 92093-0085
By applying a quasihydrostatic pressure, water or electrolyte solution can be compressed into a surface treated MSU-H mesoporous silica. Based on the pressure-volume curves, thermodynamic
and kinetic characteristics of the pressurized flow are analyzed.
For pure water based system, continuum theory explains the testing data quite well but fails to capture the rate effect. For electrolyte solution based system, the classic interface theory breaks
down, probably due to the unique ion behaviors in the
nanoenvironment. DOI: 10.1115/1.3089542
Introduction
Liquid behavior in nanoenvironment has been an important research topic for more than a decade 1,2, which shed much light
on fundamental mechanisms governing solid-liquid interactions at
small length and time scales. The studies in this area have been
immensely useful to developing better techniques for drug delivery, micro-/nanotransportation, micro-/nanofabrication, biosensing, etc.
Over the years, energy dissipation associated with nanofluidic
behaviors has drawn considerable attention 3,4. The initial idea
was based on the assumption that, as a liquid flows in nanopores
and the internal friction takes place across the large pore surface,
the energy dissipation capacity can be ultrahigh. Such a system
may have great applicability in advanced dampers, bumpers, and
armors. For instance, in an experiment performed by Li 5, bending moments were applied on a hydrophilic nanoporous silica
block soaked by water. As water flowed from compressive parts to
tensile parts, there was a detectable damping effect. However, the
energy dissipation efficiency was lower than the expected level.
Recently, an improved concept has been investigated by a few
research teams 611, including the authors of this paper, who
replaced the hydrophilic nanoporous solids by hydrophobic ones.
Due to the capillary effect, as a hydrophobic nanoporous material
is immersed in water, the nanopores remain empty. Once a high
pressure is applied on the liquid phase, pressure induced infiltration PII would occur, accompanied by conversion of a significant amount of external work to excess solid-liquid interfacial
energy. In a number of nanoporous systems, the excess solidliquid interfacial energy cannot be converted back to the hydrostatic pressure as the external loading is removed; that is, the
sorption isotherm curves are hysteretic, which results in ultrahigh
energy dissipation efficiency larger than a few J/g, much better
than many conventional protective materials 12. A similar design concept can be extended to nonpolar liquids using lyophobic
1
Corresponding author.
Contributed by the Fluids Engineering Division of ASME for publication in the
JOURNAL OF FLUIDS ENGINEERING. Manuscript received February 23, 2007; final manuscript received March 11, 2008; published online March 6, 2009. Assoc. Editor: Ali
Beskok.
Experimental Procedure
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
Pressure (MPa)
Sodium Chloride
Solution Based System
Pure Water
Based System
the liquid molecules would be blocked. When the pressure is lowered, since gas phase nucleation may be energetically unfavorable,
the confined liquid does not come out of the nanopores 17; that
is, a decrease in solid-liquid interfacial energy associated with the
expansion of a nanobubble can be lower than the energy required
for evaporating molecules from liquid phase to vapor phase,
which causes a significant hysteresis in sorption isotherm. The
enclosed area of the loading-unloading cycle indicates the dissipated energy Table 1.
According to Fig. 3, the infiltration volume of the silica surface
treated for 3 h, which can be measured by the width of infiltration
In pure water
v = 1 mm/ min In saturated NaCl solution
v = 20 mm/ min in pure water
12
24
5.1
6.8
5.2
5.6
7.0
5.6
3.1
3.6
3.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
Conclusion
To summarize, through a pressure induced infiltration experiment on a surface treated MSU-H nanoporous silica, the liquid
motion in nanopores are examined. For the system under investigation, the classic capillary theory can be employed to relate the
infiltration pressure to the nanopore size. However, the continuum
theory does not explain the low sensitivity of sorption isotherm to
loading rate. It also fails to capture the large variation in infiltration pressure caused by addition of sodium chloride.
Acknowledgment
This work was supported by the Army Research Office under
Grant No. W911NF-05-1-0288, for which the authors are grateful
to Dr. D.M. Stepp.
References
1 Sansom, M. S. P., and Biggin, P. C., 2001, Biophysics: Water at the Nanoscale, Nature London, 414, pp. 156158.
2 Majumder, M., Chopra, N., Andrews, R., and Hinds, B. J., 2005, Experimental Observation of Enhanced Liquid Flow Through Aligned Carbon Nanotube
Membranes, Nature London, 438, pp. 4446.
3 Borman, V. D., Grekhov, A. M., and Troyan, V. I., 2000, Investigation of the
Percolation Transition in a Nonwetting Liquid-Nanoporous Medium System,
J. Exp. Theor. Phys., 91, pp. 170181.
4 Fadeev, A. Y., and Eroshenko, V. A., 1997, Study of Penetration of Water Into
Hydrophobized Porous Silicas, J. Colloid Interface Sci., 187, pp. 275282.
5 Li, J. C. M., 2000, Damping of Water Infiltrated Nanoporous Glass, J. Alloys
Compd., 310, pp. 2429.
6 Martin, T., Lefevre, B., Brunel, D., Galarneau, A., Di Renzo, F., Fajula, F.,
Gobin, P. F., Quinson, J. F., and Vigier, G., 2002, Dissipative Water Intrusion
in Hydrophobic MCM-41 Type Materials, Chem. Commun. Cambridge,
2002, pp. 2425.
7 Borman, V. D., Belogorlov, A. A., Grekhov, A. M., Lisichkin, G. V., Tronin, V.
N., and Troyan, V. I., 2005, The Percolation Transition in Filling a Nanoporous Body by a Nonwetting Liquid, J. Exp. Theor. Phys., 100, pp. 385
397.
8 Surani, F. B., Kong, X., and Qiao, Y., 2005, Two-Staged Sorption Isotherm of
a Nanoporous Energy Absorption System, Appl. Phys. Lett., 87, p. 251906.
9 Surani, F. B., Kong, X., Panchal, D. B., and Qiao, Y., 2005, Energy Absorption of a Nanoporous System Subjected to Dynamic Loadings, Appl. Phys.
Lett., 87, p. 163111.
10 Kong, X., and Qiao, Y., 2005, Improvement of Recoverability of a Nanoporous Energy Absorption System by Using Chemical Admixture, Appl.
Phys. Lett., 86, p. 151919.
11 Kong, X., Surani, F. B., and Qiao, Y., 2005, Effects of Addition of Ethanol on
the Infiltration Pressure of a Mesoporous Silica, J. Mater. Res., 20, pp. 1042
1045.
12 Lu, G., and Yu, T., 2003, Energy Absorption of Structures and Materials, CRC,
Boca Raton, FL.
13 Punyamurtula, V. K., and Qiao, Y., 2007, Hysteresis of Sorption Isotherm of
a Multiwall Carbon Nanotube in Paraxylene, Mater. Res. Innovations, 11, pp.
3739.
14 Punyamurtula, V. K., Han, A., and Qiao, Y., 2007, Damping Properties of
Nanoporous Carbon-Cyclohexane Mixtures, Adv. Eng. Mater., 9, pp. 209
211.
15 Lim, M. H., and Stein, A., 1999, Comparative Studies of Grafting and Direct
Syntheses of Inorganic-Organic Hybrid Mesoporous Materials, Chem. Mater.,
11, pp. 32853295.
16 Qiao, Y., Cao, G., and Chen, X., 2007, Effects of Gas Molecules on Nanofluidic Behaviors, J. Am. Chem. Soc., 129, pp. 23552359.
17 Han, A., Kong, X., and Qiao, Y., 2006, Pressure Induced Infiltration in Nanopores, J. Appl. Phys., 100, p. 014308.
18 Fay, J. A., 1994, Introduction to Fluid Mechanics, MIT, Cambridge, MA.
19 Hartland, S., 2004, Surface and Interface Tension, CRC, Boca Raton, FL.
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