4.2 Seal surface hole-pattern periodic geometry. The surface geometry of these dimensions is used in the computational example. D/H=15.625, h/H=15.625,.
COMPUTATIONS OF GAS ANNULAR DAMPER SEAL FLOWS
By GOCHA CHOCHUA
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2002
Copyright 2002 by Gocha Chochua
To my wife, Tatyana
ACKNOWLEDGMENTS I want to express my sincere gratitude to my advisor Dr. Wei Shyy, for both his guidance and tremendous support during my Ph.D. study. I also want to thank Dr. Bruce Carroll, Dr. William Lear, Dr. Renwei Mei and Dr. Corin Segal for their assistance while serving on my committee. I am grateful for the support I received from Dr. Jeff Moore of Dresser-Rand Co.
iv
TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................. iv LIST OF TABLES ............................................................................................................ vii LIST OF FIGURES..........................................................................................................viii CHAPTERS 1
INTRODUCTION......................................................................................................... 1 1.1 1.2
Gas Annular Damper Seals in Turbomachinery Applications ............................. 1 Goals of the Present Research.............................................................................. 8
2 GOVERNING EQUATIONS AND NUMERICAL TECHNIQUES......................... 10 2.1 Governing Equations.......................................................................................... 10 2.2 Transformation to Body-Fitted Coordinates ...................................................... 13 2.3 Discretized Form of Equations........................................................................... 16 2.4 Pressure Based Solver for Navier-Stokes Equations.......................................... 18 2.5 Convection Scheme............................................................................................ 20 2.6 Multiblock Algorithm ........................................................................................ 22 2.7 Turbulence Modeling ......................................................................................... 23 2.7.1 Transport Equations for k and ε ................................................................ 23 2.7.2 Nonequilibrium k-ε Model........................................................................ 25 2.7.3 Wall Functions and Low Reynolds Number Models ................................ 26 2.7.4 Performance of the Wall Functions and Low Reynolds Number Models 29 2.8 Periodic Boundary Conditions for Gas Annular Seals Analysis........................ 34 2.9 CFD Grid Generation for Regularly Arrayed Surface Geometries.................... 40 3 DETAILED COMPUTATIONAL STUDY OF SEAL FLOWS................................ 42 3.1 Numerical Evaluation for Experimental Studies of Honeycomb Geometries ... 42 3.1.1 Experimental Setup ................................................................................... 42 3.1.2 Computational Results .............................................................................. 44 3.1.3 Loss Mechanisms ...................................................................................... 49 3.2 Assessment of the Periodic Treatment Based on the Scallop Seal Geometry ... 54 3.2.1 Problem Definition.................................................................................... 54 3.2.2 Periodic Boundary Conditions .................................................................. 55 3.2.3 Computational Results .............................................................................. 56 3.2.4 Loss Mechanisms ...................................................................................... 62 3.2.5 Extrapolation of Periodic Results to the Full Geometry ........................... 64
v
4 ROUGH WALL BOUNDARY MODEL ................................................................... 69 4.1 Overview of Flows over Rough Surfaces .......................................................... 69 4.2 Problem Definition............................................................................................. 72 4.3 Rough Model Development ............................................................................... 73 4.4 Computational Evaluation.................................................................................. 80 4.4.1 Computational Grids and Boundary Conditions ....................................... 80 4.4.2 Experimental Validation of the Detailed CFD Solutions.......................... 81 4.4.3 Flow Features in the Hole-Pattern Damper Seals ..................................... 83 4.4.4 Roughness Lengths Evaluation ................................................................. 85 4.4.5 Rough Wall Model Computations ............................................................ 90 5
CONCLUSIONS AND PERSPECTIVES .................................................................. 92
REFERENCES.................................................................................................................. 96 BIOGRAPHICAL SKETCH .......................................................................................... 103
vi
LIST OF TABLES Table
Page
2.1 Constants and functions for k-ε models .......................................................................... 25 2.2 Summary of the computed test cases. ............................................................................. 29 2.3 Periodic boundary conditions.......................................................................................... 37 3.1 Summary of the computed test cases setup and integral results...................................... 58 4.1 Computational setup........................................................................................................ 82
vii
LIST OF FIGURES Figure
Page
1.1 Gas annular seal in a centrifugal compressor. ................................................................. 2 1.2 Schematic of a traditional grooved seal, where the cavities are axisymmetric. .............. 2 2.1 Collocated grid and notation 2-D example. a) Physical plane; b) Transformed plane... 16 2.2 Schematic of the linear extrapolation employed by the second-order upwind convection scheme. ........................................................................................................ 21 2.3 Zones in the turbulent boundary layer for a typical flow over a smooth flat plate.......... 27 2.4 Channel domain dimensions. Mesh with near-wall refinement H/dy1=100 (left) and uniform mesh H/dy1=20 (right). ..................................................................................... 30 2.5 Friction factor vs. Reynolds number in a linear scale. The y+ values shown in the figure correspond to the first grid point next to the solid wall (at the distance dy1 from the wall)................................................................................................................. 31 2.6 Friction factor as in Figure 2.5 but in a logarithmic scale............................................... 32 2.7 Pressure and friction factor plots along the channel length. Re=38,100. ....................... 33 2.8 Friction factor along the channel length for long and short domains. L/H=200. ........... 35 2.9 Schematic of the periodic boundary conditions for compressible flows......................... 36 2.10 Friction factor in the full-length and periodic seal geometries...................................... 39 2.11 Mass flux, turbulent kinetic energy and eddy dissipation profiles at the inlet of the full channel, middle of the full channel, and inlet values computed at the periodic inlet................................................................................................................................. 39 3.1 Geometries of the seal with two honeycomb plates. Dimensions correspond to the test case 17............................................................................................................................ 43 3.2 The mesh system generated for the honeycomb geometry for test case 20 in the experiment. In this illustration, a total of 288,192 grid points with 5 grid blocks are employed. ....................................................................................................................... 44
viii
3.3 Comparison of friction factors obtained by the experiment and numerical computations for the honeycomb test cases 17 and 20 (pin=6.9 bar) and the smooth wall test case 2 (pin=12.4 bar). H=0.38 mm. ................................................................. 46 3.4 Streamline plots, in the symmetry plane of the honeycomb cell. Test case 17. Cell width = 0.79 mm, h = 3.05 mm, H = 0.38 mm and pin = 6.9 bar. .................................. 47 3.5 Streamline plots, in the symmetry plane of the honeycomb cell. Test case 20. Cell width = 0.79 mm, h = 3.81 mm, H=0.38 mm and pin= 6.9 bar. ..................................... 48 3.6 Representative pressure contours on the cell walls. fp=0.077. h=3.05 mm, Re=6,302. 50 3.7 Contours of near-wall shear stress absolute values at the cell walls. fτ=0.005. h=3.05 mm, Re=6,302................................................................................................................ 51 3.8 Pressure contours in Y-Z plane through the plane of symmetry. Flow is from left to right. h=3.05 mm, Re=6,302. ......................................................................................... 52 3.9 Shear stress contours in Y-Z plane through the plane of symmetry. Flow is from left to right. h=3.05 mm, Re=6,302. ........................................................................................ 52 3.10 Typical pressure plot at the middle plane between two honeycomb plates. Flow is in positive Z direction. Periodic domain with 4 cells is shown twice, using advantage of sidewise symmetry..................................................................................................... 53 3.11 Temperature contours in Y-Z plane through the plane of symmetry. Flow is from left to right. h=3.05 mm, Re=6,302............................................................................... 53 3.12 Scallop pattern on the stator surface.............................................................................. 54 3.13 Full-length geometry in the axial (Z) direction (left) and the periodic domain with 4 near-outlet cavities. Bold line shows the middle surface providing data for the periodic boundary conditions......................................................................................... 56 3.14 Full-length computational grid and boundary conditions. ............................................ 57 3.15 Streamlines in the half seal clearance for the scallop-stator gas annular seals with pressure ratio 0.4. Top: preswirl=0; Bottom: preswirl=0.27. Plots are constructed by using three computational domains in the sidewise direction, which have periodic boundary conditions. ...................................................................................................... 59 3.16 Averaged axial Mach number and swirl ratio in the seal for the cases with pressure ratio 0.4. ......................................................................................................................... 60 3.17 Averaged axial Mach number and swirl ratio in the seal for the cases with pressure ratio 0.6. ......................................................................................................................... 60 3.18 Typical streamline pattern in a scallop cavity. .............................................................. 61
ix
3.19 Vector field in the middle of the cavity. Axial (main) flow is in the positive Z direction. Rotor wall motion is in the positive X direction. Velocity vectors are at scale on both figures....................................................................................................... 61 3.20 Pressure contours on the scallop walls. fp=0.013. Pressure ratio is 0.4....................... 63 3.21 Shear stress contours on the scallop walls. fτ=0.0037. Pressure ratio is 0.4. .............. 63 3.22 Pressure plot at the scallop seal rotor surface, periodic results. .................................... 64 3.23 Example of 3 geometrically similar points.................................................................... 65 3.24 Static pressure from the full geometry computations of test case 1 (solid) compared to self-extrapolated results based on data from 3 checkerboard cells. Symbols represent cells, which are employed in the extrapolation. ............................................. 66 3.25 Static pressure from the full geometry computations of test case 1 (solid) compared to extrapolation from the periodic computations. Top: 4-cell periodic computations; bottom: 8-cell periodic computations. ........................................................................... 67 4.1 Schematic of the flow and definitions............................................................................. 70 4.2 Seal surface hole-pattern periodic geometry. The surface geometry of these dimensions is used in the computational example. D/H=15.625, h/H=15.625, Dm/H=16.52.................................................................................................................... 72 4.3 Horizontal (top), vertical (middle) velocities and the dispersive flux (bottom) distributions over the hole-pattern roughness (y=0). Re=105........................................ 74 4.4 Shear stress distribution between the rotor and seal stator, normalized by the sum of the shear stresses on both surfaces. Re=105. Rotor speed, ω=0. The total stress is composed of the Reynolds stresses and stresses induced by the dispersive flux. Turbulent kinetic energy has a smooth maximum at some distance above the seal surface. ........................................................................................................................... 77 4.5 Computational mesh and boundary conditions. Top: Full-cavity 3-D computational grid. h/D=1, D/H=15.625. Bottom: 2-D approximation with rough wall boundary conditions representing the cavity effect........................................................................ 81 4.6 Dependence of the hole-pattern Reynolds number Re=ρU(2H)/µ (nondimensional seal leakage) and swirl ratio from the hole depth........................................................... 82 4.7 Streamlines colored by absolute velocity levels for two hole depths: h/D=0.2 (left) and h/D=0.48 (right). ..................................................................................................... 83 4.8 Velocity magnitude (left) and turbulent kinetic energy, k, (right) contours, obtained from the full geometry computations and normalized by the seal surface frictional
x
velocity, u*, and the streamline pattern resulted in the cavity. The flow is from left to right. The cavity is cut along its axis of symmetry........................................................ 84 4.9 Pressure distribution on the rotor surface. Cavities create discrete pressure drops near the downstream edges. ................................................................................................... 85 4.10 Horizontal velocity, turbulent kinetic energy and dissipation rate of the turbulent kinetic energy (all area averaged) for three Reynolds numbers. All values are normalized by the rough wall frictional velocity u* and friction length y0=0.00335 mm. The coordinate origin is shifted downwards by the zero-plane displacement, d/y0=10.3. ....................................................................................................................... 86 4.11 Pressure (top) and the shear stress component in the axial (x) direction (bottom) contours on the seal stator surface. The high-pressure region appears near the downstream edge of the cavity, as the separated flow reattaches. The pressure (form) drag contributes to 97.6 % of the total stress on the seal surface................................... 88 4.12 Plot of the area-averaged velocity computational results, normalized by the seal friction velocity and the theoretical logarithmic profile with d/H=0.17, y0/H=0.0165, Re=105 and h/D=1. Right figure is a blown up view of the top portion of the left figure. The logarithmic profile is extrapolated down to the apparent wall origin. ....... 89 4.13 Comparison of spatially averaged mean velocity, turbulent kinetic energy and eddy dissipation rate between the rough wall model and full cavity computations for Re=105. All values are nondimensionalized by the stator surface frictional velocity u* and the friction length, y0=0.00335 mm. The zero length displacement d/y0= 10.3. 91
xi
NOMENCALTURE A
coefficients in the discretized equation (2.36)
Aw
effective horizontal area
a,b,c
model coefficients
C
direct damping defined by Eq. (1.1)
C
convective coefficients
Ceff
effective damping, C-k/ω.
Cµ, Cε1, Cε2
turbulence closure coefficient
c
cross-coupled damping defined by Eq. (1.1)
D
hole diameter
D
near-wall damping function for ε
Dm
mean element separation distance,
d
zero plane displacement or error in origin.
E
near-wall function used by low Reynolds number models
E, F, G
Variables and fluxes in compact form of Eq. (2.22)
Eij
mean rate of deformation tensor defined by Eq. (2.52)
e
rotor displacement amplitude
F
general flux at the control volume interface
Fx, Fy
reaction forces on a rotor
f
friction factor defined by Eq. (1.3).
f11, f12, etc.
metrics coefficients defined by Eq. (2.20)
xii
Aw n
fp
friction factor due to pressure loses defined by Eq. (3.2)
fµ
closure coefficient defined by Eq. (2.48)
fτ
friction factor due to loses by shear stresses defined by Eq. (3.3)
H
seal clearance
H
total enthalpy per unit mass, h + 12 ui u j + k
h
roughness height or cavity depth
h
enthalpy per unit mass
J
Jacobian defined by Eq. (2.21)
K
direct stiffness defined by Eq. (1.1)
k
cross-coupled stiffness defined by Eq. (1.1)
k
kinetic energy of fluctuations per unit mass
L
seal length
l0
turbulent length scale
lz
roughness element transverse (z) dimension
M
direct added mass defined by Eq. (1.1)
M
Mach number
Mz
Mach number in the Z direction, based on W velocity
m!
mass flow rate per unit width
N(ui)
Navier-Stokes operator
n
number of roughness elements in a horizontal area A
n,m
empirical coefficients defined by Eq. (1.3)
ni
wall normal unit vector components
PrL
laminar Prandtl number
Prt
turbulent Prandtl number
p
pressure
xiii
q
heat flux
q11, q12, etc
transformation metrics defined by Eq. (2.28)
R
gas constant
Re
Reynolds number,
Ret
turbulent Reynolds number, ρk2/µε
r
position vector in vector notation
S
source term
T
absolute temperature
Tt
total temperature
Txy
shear stress component in the XY plane due to the viscous and Reynolds stresses defined in Eq. (4.9)
Ti
turbulence intensity
t
time
t0
turbulent time scale
ti
near-wall shear stress unit vector components
U
mean flow velocity in the seal in X direction.
U, V, W
Contravariant velocities defined by Eq. (2.29)
u, v, w
Cartesian Favre averaged velocity components
u*
friction velocity, τ w ρ
V
Velocity vector magnitude,
W
mean flow velocity in the seal in Z direction.
X, Y
rotor displacement around the centered position
X, Y, Z
axis of the rectangular Cartesian coordinate system
x, y, z
rectangular Cartesian coordinates
ρU ( 2 H ) µ
u 2 + v 2 + w2
xiv
Y
displaced height, y+d
Yp
distance from a near-wall mesh point, P, to the wall
y0
roughness length
y*
lower limit of the inertial sublayer
y
distance from the wall normalized by the clearance, y/H
Greek Symbols
β
mesh stretching parameter
∆
incremental change
δ
boundary layer thickness
δij
Dirac delta function
ε
turbulent kinetic energy dissipation rate
ε"
effective turbulent kinetic energy dissipation rate, ε - D
Φ
dissipation function
φ
any scalar or vector component or general dependent variable
Γ
effective diffusivity
γ
ratio of specific heats
κ
von Karman constant
λ
thermal diffusivity
λ
roughness density, hl z Dm2
µ
molecular viscosity
µt
eddy viscosity defined by Eq. (2.13)
ν
kinematic molecular viscosity, µ/ρ
Π
generation (production) of k
xv
ρ
gas density
σk, σε
turbulence closure coefficients
τij
shear stress tensor
τw
wall shear stress
ω
angular velocity
ξ, η, ζ
transformed coordinates
Subscripts
E,W,N,S,T,B neighboring cell-centered points for the node P as shown in Figure 2.1 e,w,n,s,t,b
neighboring points at the interfaces for the node P as shown in Figure 2.1
i, j, k, l
cyclical indices
i, j, k,ii, jj, kk indexes of structured grid in
denotes inlet boundary
mid
denotes middle of the domain
nbr
any neighboring point at the interface for the node P
out
denotes outlet boundary
P
cell centered point as shown in Figure 2.1
ref
reference values
ξ, η, ζ
partial differentiation with respect to ξ, η or ζ.
Superscripts
u, v, w
Cartesian velocity components
′
fluctuating component above the Reynolds time averaged mean
′
correction to the intermediate step solution in the pressure based algorithm
″
fluctuating component above the mass-weighted time averaged mean xvi
′′′
fluctuating component above the area averaged mean
`
fluctuating component above the mass-weighted area averaged mean
+
sublayer-scaled value
*
intermediate step solution
Overbars
.
first derivative with respect to time
..
second derivative with respect to time
-
Reynolds temporal average
~
mass-weighted, Favre, temporal average or nondimensional variable
Other
area average
{}
mass-weighted area average
xvii
Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONS OF GAS ANNULAR DAMPER SEAL FLOWS By Gocha Chochua May 2002 Chairman: Dr. Wei Shyy Major Department: Aerospace Engineering, Mechanics and Engineering Science Gas annular damper seals represent a new design approach in turbomachinery applications offering high sealing effectiveness and enhanced rotordynamic stability. Computational and modeling capabilities are developed for analyzing compressible, turbulent flows in the present type of seals. Due to repeated geometric patterns of the seals, it is desirable to develop suitable strategies to treat a small portion of the flow domain using periodic boundary conditions and extrapolate the information to the entire seal. The detailed computational analyses conducted include (i) development of an original periodic treatment for compressible turbulent flows around repeated geometries, (ii) 3-D investigation of velocity and pressure characteristics, and loss mechanisms under given rotor speed and inlet swirl ratio, and (iii) assessment of the extrapolation technique of the full-length results from computations based on a limited periodic domain. Based on the insight gained from the computations in detailed geometries, a simplified engineering model utilizing the roughness treatment is proposed to significantly reduce the computing cost. To develop the concept, the k-ε two-equation
xviii
turbulence closure with extended wall functions is developed. The approach is based on the concept of spatial averaging in the plane parallel to the boundary of interest. A main aspect of the model is that the mean momentum transport, resulting from the wall-normal velocity component, is important for evaluating the near-boundary shear stresses. The present roughness model helps assess the effect of geometry, size and pattern of the seal cavity on the overall flow structures. The computational techniques presented can help address fluid physics and design practices for gas annular seal flows.
xix
Equation Section 1
CHAPTER 1 1 INTRODUCTION 1.1
Gas Annular Damper Seals in Turbomachinery Applications
Gas annular seals are commonly adopted for turbomachinery applications to reduce leakage, thus improving efficiency. A schematic of the seal in a centrifugal compressor is shown in Figure 1.1. Turbo centrifugal compressors are used to elevate gas pressure for example for the transmission of gas in a pipeline. The centrifugal compressor shown has two stages placed in back-to-back scheme to balance thrust on the rotor. After passing through the first stage with 3 impellers the gas is taken out of the compressor. After being cooled it passes trough the second stage. The seal stator surrounds the rotor without contacting it. Its main purpose is to withstand the differential pressure between two stages with minimal leakage. The simplest seals have smooth seal surface. Widely used labyrinth and grooved seals, shown in Figure 1.2, have axisymmetric cavity shapes. They combine good leakage performance and low cost. Seals have small clearances H between the rotor and seal, being order of 0.3% of the seal diameter. Series of grooves or labyrinths serve to reduce the internal leakage. However, operating at high pressures and tight clearances, such seals may develop significant force imbalance to cause rotordynamic instabilities. Von Pragenau [1] suggests that the circumferential velocity in the seal (swirl) in the axisymmetric geometries is the source of the instabilities. The author proposes a damper
1
2 seal configuration, employing a smooth rotor and deliberately roughened stator, which in addition to leakage reduces the swirl.
Gas annular seal stator. Differential pressure between two stages.
Rotor
First stage, then gas is cooled
Second stage
Figure 1.1. Gas annular seal in a centrifugal compressor [2].
Figure 1.2. Schematic of a traditional grooved seal, where the cavities are axisymmetric.
3 A damper seal with honeycomb surface geometry is applied to improve rotordynamic performance in HPOTP (High Pressure Oxygen Turbopump) of the Space Shuttle main engine (see Scharrer [3]). Numerous other damper seal surface geometries are proposed and analyzed, as will be discussed later. The customary used model to define reaction forces (Fx, Fy) for small motion about a centered position (X, Y) uses the following linearized force-displacement equation: F K − X = FY − k
k X C c X! + + M K Y − c C Y!
X!! !! Y
(1.1)
where K, C stands for the direct and k, c for the cross-coupled stiffness and damping respectively. M is direct added mass coefficient. As related to rotordynamics, seal analysis has the objective of defining these rotordynamic coefficients. For gas seals, emphasized in this paper, M and c are typically neglected. The roughened seal stator in addition to the leakage reduces the circumferential velocity in the seal (swirl) and hence the cross-coupled damping coefficient k. Thus, replacing smooth or labyrinth seal with a damper seal increases onset speed of instability. The effect of k is better seen if Eq. (1.1) is rewritten for radial and tangential force components acting on a rotor that is whirling synchronously at constant angular speed ω and amplitude e: Fn K + cω − M ω 2 − = e Cω − k Ft
(1.2)
Childs and Kim [4] show that the effective damping coefficient, Ceff=C-k/ω, has higher values for damper seals because of reduction in k. Childs et al. [5] have conducted extensive experimental investigations for honeycomb seals, and show that the best sealing and rotordynamic performance for
4 damper seals, with swirling incoming flows and followed by labyrinth seals, is achieved for seals that are longer than around 50 mm. However, Kleyhans [6] shows that for shorter seals, say, 25 mm in length, the rotordynamic stability advantage is diminished and swirl brakes must be used. A typical seal is characterized by small clearances between rotor and seal stator, compared to the diameter, as shown in Figure 1.2. The damper seal stator surface is covered with arrayed cavities, producing flow expansion/contraction in both axial and circumferential directions and thus reducing leakage and swirl. A typical damper seal surface geometry, (see, e.g., Figure 3.1) has thousands of three-dimensional cavities, whose dimensions are order of magnitude larger than the seal clearance. Hence, discrete cavities affect fluid physics in the entire clearance region. On the other hand, the clearance length-to-height ratio is high. These aspects complicate analytical studies of the seals. Traditionally, the analytical study of the flow in damper seals is based on the turbulent bulk flow model of Hirs [7]. The theory adopts the Blasius-type pipe friction model for the friction factor, f: ρU ( 2 H ) τw m f = = n ( Re ) = n 1 ρU 2 µ 2
m
(1.3)
where U is the bulk flow velocity relative to the rotor or stator walls, ρ is the fluid density, τw is the wall shear stress. The Reynolds number, Re, is based on the hydraulic diameter, which is twice the seal clearance, H. Empirical coefficients n and m are to be assigned from an experiment. Well-known values for smooth pipes are n = 0.079 and m = -0.25. Nelson [8] provided the first complete analysis, based on Hirs bulk flow
5 theory, for the leakage and rotordynamic coefficients of Eq. (1.1) of tapered annular gas seals operating in turbulent regime. Nelson’s approach is based on zeroth and first order perturbation equation. For non-prerotation the theory gives good agreement. However, for the prerotated case, the agreement is not as good. Elrod et al. [9] adjust the model for honeycomb seals by using the appropriate entrance region friction factor model. Von Pragenau [1] and Nelson and Nguyen [10] propose a Moody model for the damper seals. The Moody model uses the “relative roughness” parameter h/(2H), where h is the roughness height (cavity depth): 1 6 3 h 10 f = 0.0055 1 + 104 + 2 H Re
(1.4)
DeOtte et al. [11] present friction factor measurements for their hole-pattern flat plate tester. The results show that increasing clearance for the same Reynolds number, which presumably reduces the relative roughness, increases the friction factor, which is unrealistic. Ha and Childs [12] propose a two control volume model, where the additional control volume is used to simulate the cavity cell depth and allows flow to only enter end exit in the radial direction. Kleynhans [13] develops an alternative to the Hirs’ viscous shear stress model. The alternative model treats the viscous shear stress as nondirectionally homogeneous. Kleynhans and Childs [14] propose an acoustic based bulk flow model. A honeycomb cell is modeled as “capacitance accumulator” that caused a dramatic reduction in the effective acoustic velocity of the mean flow, dropping the lowest acoustic frequency into the frequency range of interest of rotordynamics. Overall, comparison of the analytical results with measurements shows that there is still much
6 room for improvement. Childs [15] suggests that “further refinements of pipe friction models are not productive, and research effort should be devoted to a clearer understanding of the fluid mechanics involved in annular seal flow.” Computational fluid dynamics (CFD) is being actively pursued in rotordynamic analysis and design. While the CFD approach requires more computational resources, it allows one to account for the seal specific fluid physics in detail. CFD analysis of rotordynamic coefficients in damper seals, characterized by 3-dimensional surface cavities, represents an extremely computationally expensive problem and has not yet been performed. Several authors have studied smooth and axisymmetric (labyrinth, grooved etc.) seals. For example, Williams et al. [16] perform three-dimensional moving boundary Navier-Stokes computations for time accurate incompressible fluid reaction force calculations in a smooth seal. The original k-ε turbulence model [17] is applied. The key element of the algorithm is solving the Helmholtz pressure equation in order to help enforce the continuity equation. First, steady state solution is obtained for geometry of interest. The steady state solution is then used as an initial condition for the unsteady calculations. The rotor shaft moves with an imposed harmonic motion and unsteady calculations are carried out until the time periodicity is resolved. The force time history is then processed to obtain the rotordynamic coefficients. Moore and Palazzolo [18] apply a pressure-based incompressible algorithm to calculate the flow field and rotordynamic forces in a whirling, grooved liquid seal. The whirling, unsteady problem is transformed into a steady one by solving the threedimensional, eccentric seal flow field in the rotating frame of reference, attached to the whirling rotor. Quasi-steady solutions are obtained for multiple whirl frequency ratios.
7 Predictions are made using both a original and low Reynolds number k-ε model. Authors report good qualitative agreement with the experiment. It should be noted that although the quasi-steady approach is computationally efficient, it is limited to axisymmetric geometries and hence may not be applied to damper seal computations. Numerous experimental studies are performed for seal flows. For example, Ha and Childs [19] and Ha [20] present interesting and comprehensive results for flows between two parallel honeycomb plates. They demonstrate that the friction factor is sensitive to changes in cell geometry, making it necessary to adjust empirical coefficients of the theories for an individual surface geometry. In addition, experimental studies show a strong dependence of the friction factor from the Reynolds number in certain situations, including dramatic drop and rise of the friction factor. This phenomenon, termed “friction factor jump," is observed in 42% of the test cases. It is explained by interaction between shear flow and acoustic resonance in honeycomb cells [21]. Hendricks et al. [22] have reported their Full Navier-Stokes solutions for flows between two flat honeycomb plates, based on the experimental studies by Ha and Childs [19]. Computations are performed in a small periodic domain. However, little quantitative information is provided to allow detailed evaluation of their efforts. More detailed description of the experimental setup and computational results are presented in Section 3.1. Nava [23] experimentally studies friction factor in incompressible flow between two parallel plates of different roughness patterns in combination with a smooth plate. Considered surface geometries include hole-pattern and knurl- (diamond shape) pattern of different scales. Clearances range from 76 to 1270 µm and Reynolds number is up to
8 105, leading to the both laminar and transition regimes. As the clearance between two plates is increased, the friction factor increases until it reaches a plateau. Beyond this point, further increase in clearance lead to reduction in the friction factor. Carter et al. [24] experimentally study leakage and rotordynamic coefficients for smooth, hole-pattern and honeycomb gas annular seals. Six tests are conducted for each seal, representing 3 rotor speeds and 2 pressure ratios. The results show that geometric aspect ratio, more than cavity shape, has a significant effect on the seal performance. For example, in their study, increasing the hole depth results in a noticeable decrease in the seal leakage. The experimental setup and results as well as the effect of the cavity proportions on the seal efficiency are discussed in Chapter 4. Sprowl et al. [26] perform similar studies for the scallop geometry, which is numerically studied in Section 3.2. The results show that the scallop geometry results in higher leakage but lower effective stiffness than a typical hole-pattern seal. This property may be advantageous in some practical applications, when reduced system stiffness is required. 1.2
Goals of the Present Research
In the present study, computational fluid dynamics (CFD) is adopted to utilize the damper seal flow problem. Computational techniques are refined from both numerical and modeling viewpoints. The main goals of the current research are the following: •
Develop a computational technique for damper seal analysis by considering a limited in length and circumference computational domain.
•
Study fluid mechanics involved in the gas annular damper seal flows.
•
Develop and validate a rough (generalized) wall functions boundary model under the k-ε turbulence model for the damper seal analysis.
9 A numerical scheme for handling periodic boundary conditions for compressible, turbulent flows around geometries with repeated patterns is developed. After the technique is developed and validated, detailed numerical computations allow one to study the flow physics in the seals, including effects of inlet swirl and rotor velocity. Several surface geometries and regimes are chosen in order to evaluate the effect of surface geometry on flow physics. Extrapolation techniques of the periodic solutions enable one to analyze the fluid flow in the full-length geometry based on computations in reduced domains. Based on the computational analysis, a rough surface model, which smears out the details of the cavities, is developed and assessed. The proposed rough boundary model is based on simplified treatment of the seal surface as a rough surface. As already mentioned, the major difficulty in this approach is that the seal roughness dimensions are substantially higher than the main flow length scales. Approaches developed in the meteorological literature are considered and adopted in the present study. Results obtained using the generalized wall treatment are then assessed by the computations accounting for detailed cavity geometries. The resulting generalized wall boundary treatment can significantly reduce the computational costs, and allow the damper seals to be analyzed with simplified resources as those for the smooth seals. Finally, Chapter 5 presents conclusions and perspectives.
Equation Section (Next)
CHAPTER 2 2 GOVERNING EQUATIONS AND NUMERICAL TECHNIQUES 2.1
Governing Equations
The present study utilized a computational tool (Thakur and Wright [27]). Basic governing equations, numerical techniques and modifications proposed by this study are shown next. The Navier-Stokes equations governing compressible fluid flow are the continuity, momentum, total enthalpy equations and equation of state ∂ρ ∂ + ( ρ ui ) = 0 ∂t ∂xi
(2.1)
∂ ∂ ∂p ∂τ ( ρ ui ) + ( ρu j ui ) = − + ij ∂t ∂x j ∂xi ∂x j
(2.2)
∂q ∂ ∂ ∂ ( ρ H − p ) + ( ρ u j H ) = − j + ( uiτ ij ) ∂t ∂x j ∂x j ∂x j
(2.3)
p = ρ RT
(2.4)
Time averaged form of the above equations for turbulent flows is obtained by Reynolds time averaging for density and pressure (its mean and fluctuating components are denoted by bar and prime (′) respectively). Favre-averaging (Favre [28], see also Cebeci and Smith [29] or Tannehill et al. [30]) is convenient to use for velocity components and thermal variables for compressible flows (~ denotes Favre averaged mean quantities defined as φ" = ρφ ρ and double prime (″) denotes their fluctuating components): 10
11
ρ = ρ + ρ′
p = p + p′
T = T" + T ′′ q j = q" j + q′′j
ui = u"i + ui′′ τ ij = τ"ij + τ ij′′
(2.5)
After time averaging the governing equations become ∂ρ ∂ + ( ρ u"i ) = 0 ∂t ∂xi
(2.6)
∂ ∂ ∂p ∂ ( ρ u"i ) + ( ρ u"iu" j ) = − + τ"ij + τ ij′′ − ρ ui′′u′′j ∂t ∂x j ∂xi ∂x j
(2.7)
∂ ∂ ∂ ρ H" − p ) + ρ u" j H" ) = − ( ( ∂t ∂x j ∂x j ∂ ∂x j
(
q" + q′′ + ρ u′′h′′ + j j j
)
u" τ" + τ ′′ − ρ u′′u′′ + ∂ u′′τ" + u′′τ ′′ − ρ u′′ 1 u′′u′′ i j i j i i ij i ij ij ∂x j i ij 2
1 where H" = h" + u"i u"i + k" ; 2
(2.8)
1 H ′′ = h′′ + u"i ui′′ + ρ ui′′ui′′ . Term − ρ ui′′u′′j is called 2
Reynolds stress tensor and it usually gives the strongest contribution from the fluctuating components. Terms in equations (2.7) and (2.8) underlined by double lines are not equal to zero, in general, unless ρ ′ = 0 . Instead, the time average of the doubly primed fluctuations multiplied by the density is equal to zero:
ρφ ′′ = 0
(2.9)
The underlined terms are typically neglected. However, for quite high Mach number effect of these terms becomes noticeable and may not be omitted. Further discussion and some proposed models might be found in Shyy et al. [31]. Expressing mean fluxes in terms of the gradients of mean variables yields − q" j = λ
∂T µ ∂h" = ∂x j PrL ∂x j
(2.10)
12 ∂u"i
τ ij = µ
∂x j
+
∂u" j ∂xi
−
2 ∂u"l ⋅ δ ij 3 ∂xl
(2.11)
where PrL, is the laminar Prandtl number and λ is the thermal diffusivity. Employing the eddy viscosity hypothesis, which relates the Reynolds stress and turbulent thermal stresses to the gradient of mean variables, the various correlation terms are modeled as follows: − ρ u′′j h′′ =
µt ∂h"
(2.12)
Prt ∂x j
∂u" ∂u" 2 ∂u" 2 − ρ ui′′u′′j = µt i + j − µt l ⋅ δ ij − ρ k" ⋅ δ ij ∂x 3 j ∂xi 3 ∂xl
µ 1 ui′′τ ij′′ − ρ u′′j ui′′ui′′ = µ + t σk 2
(2.13)
∂k" ∂x j
(2.14)
where µt is turbulent eddy viscosity, Prt is the turbulent Prandtl number and σk is a turbulence model closure coefficient. It should be noted that equations (2.12)-(2.14) are based on the Boussinesq eddy viscosity hypothesis, which assumes isotropy of the Reynolds stress tensor. Thus, with the modeling assumptions discussed above, the governing equations are expressed as ∂ρ ∂ + ( ρ u"i ) = 0 ∂t ∂xi ∂ ∂ ∂p ∂ ( ρ u"i ) + ( ρ u"iu" j ) = − + ∂t ∂x j ∂xi ∂x j
(2.15)
∂u" ∂u" j 2 ∂u"l δ ij − ( µ + µt ) i + ∂x j ∂xi 3 ∂xl
∂ ∂ ∂ µ µ ∂h" + t ρ H" − p ) + ρ u" j H" ) = ( ( + ∂t ∂x j ∂x j PrL Prt ∂x j ∂u" ∂u" j 2 ∂u"l ∂ ∂ + − ⋅ δ ij + u"i ( µ + µt ) i + ∂x ∂x ∂x j j j ∂xi 3 ∂xl
µt µ + σk
∂k ∂x j
(2.16)
(2.17)
13 p = ρ RT"
(2.18)
In the following discussions all flow variables are time-averaged. Time averaging symbols are dropped to simplify the equations. 2.2
Transformation to Body-Fitted Coordinates
For arbitrary shaped geometries, the above governing Navier-Stokes equations are transformed into generalized curvilinear coordinates (ξ, η, ζ), where ξ=ξ(x,y,z),
η=η(x,y,z) and ζ=ζ(x,y,z). The transformation of the physical domain (x,y,z) to the computational domain (ξ,η,ζ) is achieved by the following relations (see Tannehill et al. [30]): ξ x ξ y ξ z f11 1 η x η y η z = J f 21 ζ x ζ y ζ z f 31
f13 f 23 f 33
f12 f 22 f 32
(2.19)
where the metrics are f11 = yη zζ - zη yζ f12 = zη xζ - xη zζ f13 = xηyζ - yη xζ f12 = zξ yζ - yξ zζ f22 = xξ zζ - zξ xζ f23 = yξ xζ - xξ yζ
(2.20)
f31 = yξ zη - zξ yη f32 = zξ xη - xξ zη f33 = xξ yη - yξ xη
and J is the Jacobian determinant of the transformation given by xξ ∂ ( x, y , z ) = yξ J= ∂ (ξ ,η , ζ ) zξ
xη
xζ
yη zη
yζ zζ
(2.21)
Transforming each term in equations (2.15), (2.16) and (2.17) to the generalized body-fitted curvilinear coordinate system (ξ, η, ζ), the following set of equations written in compact form results in ∂W ∂E ∂F ∂G + + + =S ∂t ∂ξ ∂η ∂ζ
(2.22)
14 where ρ ρu W = ρv ρw ρ H
(2.23)
ρU Γ ρUu − ( q11uξ + q12uη + q13uζ ) + f11 p J Γ ρUv − ( q11vξ + q12vη + q13vζ ) + f12 p E= J Γ ρUw − ( q11wξ + q12 wη + q13wζ ) + f13 p J Γh Γ q11hξ + q12 hη + q13hζ ) − k ( q11kξ + q12 kη + q13kζ ) ( ρUH − J J
(2.24)
ρV Γ ρVu − ( q21uξ + q22uη + q23uζ ) + f 21 p J Γ ρVv − ( q21vξ + q22vη + q23vζ ) + f 22 p F= J Γ ρVw − ( q31wξ + q32 wη + q33wζ ) + f 23 p J Γh Γk q21hξ + q22 hη + q23hζ ) − q21kξ + q22 kη + q23kζ ) ( ( ρVH − J J
(2.25)
ρW Γ ρWu − ( q31uξ + q32uη + q33uζ ) + f 31 p J Γ ρWv − ( q31vξ + q32vη + q33vζ ) + f 32 p G= J Γ ρWw − ( q31wξ + q32 wη + q33wζ ) + f 33 p J Γh Γk q31hξ + q32hη + q33hζ ) − q31kξ + q32 kη + q33kζ ) ( ( ρUH − J J
(2.26)
15 0 Su S = Sv S w Φ
(2.27)
where
q11 = f112 + f122 + f132 q12 = q21 = f11 f 21 + f12 f 22 + f13 f 23 q 22 = f 212 + f 222 + f 232 q13 = q31 = f11 f13 + f12 f 32 + f13 f 33
(2.28)
q 33 = f 312 + f 322 + f 332 q23 = q32 = f 31 f 21 + f 32 f 22 + f 33 f 23 and Γ = µ + µ t
Γh =
µ PrL
+
µt Prt
Γk = µ +
µt . σk
Here U, V and W are the contravariant velocities given by U=f11u+f12v+f13w V=f21u+f22v+f23w
(2.29)
W=f31u+f32v+f33w
Momentum source terms (Su, Sv, Sw) absorb body forces. The dissipation function is denoted by Φ. When solving the above equations in curvilinear coordinates, even though the formal coordinate transformation is introduced, in practice, the finite-volume formulation is adopted. To achieve this goal, combined Cartesian and contravariant velocity components are employed. Specifically, while the Cartesian velocity components are the primary variables, the contravariant velocity components are used in establishing the flux estimation at the cell faces, and in enforcing the mass continuity in the pressurecorrection equation. Moreover, the contravariant components are computed by explicitly accounting for the geometric meaning of the metric terms, as discussed in Shyy [32].
16 The same is true when constructing the second-order upwind convection scheme (Thakur and Shyy [33]). 2.3
Discretized Form of Equations
The governing equations are discretized in the structured grid. All variables are arranged in the collocated manner (Rhie and Chow [34], Peric et al. [35], Miller and Schmidt [36]). In terms of the notation shown in Figure 2.1, for a node P enclosed in its cell and surrounded by neighbors N, S, E, W, T, B the finite difference approximation to the governing equations may be obtained by taking the integral of the equations over the control volume with dimensions ∆ξ=1, ∆η=1, ∆ζ=1.
Figure 2.1. Collocated grid and notation 2-D example. a) Physical plane; b) Transformed plane. The discretized form of the continuity equation (without the time dependent term) can be written as
[ ρU ]ew + [ ρV ]ns + [ ρW ]bt = 0
(2.30)
17 where [i ]w = ( i ) e − ( i ) w , etc. For later use, we denote the mass fluxes at each control e
volume face by F and rewrite the continuity equation as
[ F ]ew + [ F ]ns + [ F ]bt = 0
(2.31)
The discretized form of u-momentum equation can be obtained in similar manner:
[ ρUu ]ew + [ ρVu ]ns + [ ρWu ]bt = − [ f11 p ]ew + [ f 21 ]ns + [ f 31 ]bt + e
n
Γ ∂u ∂u ∂u Γ ∂u ∂u ∂u + q12 + q13 + q22 + q23 q11 + q21 + ∂η ∂ζ w J ∂ξ ∂η ∂ζ s J ∂ξ
(2.32)
t
Γ ∂u ∂u ∂u + q32 + q33 q31 ∂η ∂ζ b J ∂ξ
The three terms on the left hand side of equation (2.32) are the convective fluxes at the control volume faces; the first three terms on the right hand side are the pressure fluxes and the last three terms are the diffusion fluxes. The above equations can be put in the generalized form as Γ Γ ρUφ − J ( q11φξ + q12φη + q13φζ ) + ρV φ − J ( q21φξ + q22φη + q23φζ ) + w s e
n
Γ ρW φ − J ( q31φξ + q32φη + q33φζ ) = S b t
(2.33)
where φ is the general dependent variable and S is the source term. The discretized equation relating the variable at a center point P to its neighboring values is obtained as APφ P = AEφ E + AW φW + ASφS + ANφ N + ABφ B + ATφT + S =
∑A
nbr
⋅ φnbr + S
(2.34)
nbr
where pressure terms, cross-derivative portion of the viscous terms, due to nonorthogonal grid effect, and higher-order contributions of the convective fluxes are taken into the source term, S.
18 2.4
Pressure Based Solver for Navier-Stokes Equations
The SIMPLE algorithm (Patankar [37]) is essentially a guess-and-correct procedure for the calculation of pressure fields. This algorithm was originally developed for incompressible flows and has been successfully extended for compressible flows by taking into account density dependence (Shyy and Braaten [38]). The algorithm also has been extended to body-fitted curvilinear coordinates. Pressure based methods can handle flows at all speeds. Clearly, it would be easier in coding procedures if one uses a grid system that places all the flow variables, scalar and vector, at the same physical location (collocated grid). The use of such a collocated grid system has been shown to produce nonphysical oscillations in the pressure field. This situation can permit two different pressure fields to coexist, giving rise to the nonphysical “checkerboard” field (Patankar [37]). Rhie and Chow [34] developed a momentum interpolation method in which cellface velocity, being interpolated from the neighboring cell centers, also contains pressure of the neighboring nodes. Suppose that the velocity field at the intermediate step and pressure field are given by u*, v*, w*, p*. The corrected velocities and pressure can then be obtained by adding the corrections as follows: p=p*+p′ u=u*+u′ v=v*+v′ w=w*+w′
(2.35)
The u-momentum equation for this intermediate velocity field is written as * APu uP* = ∑ Anbr ⋅ unbr + S P* + SD + SC nbr
where S P* = − ( f11 pξ* + f 21 pη* + f 31 pζ* ) , SD, SC are respectively pressure terms, crossderivative diffusion terms and higher order terms of the convective fluxes.
(2.36)
19 Writing the same equation for the correct pressure and velocity field, subtracting (2.36) from it and ignoring all the terms on the right hand side except pressure terms (note that omission of these terms is the main approximations of the SIMPLE algorithm) yields 1 u′ = − u ( f11 pξ′ + f 21 pη′ + f 31 pζ′ )e AP 1 v′ = − v ( f12 pξ′ + f 22 pη′ + f 32 pζ′ )e AP
(2.37)
1 w′ = − w ( f13 pξ′ + f 23 pη′ + f 33 pζ′ )e AP
The corresponding corrections for the contravariant velocities are given by U′=f11u′+f12v′+f13w′ V′=f21u′+f22v′+f23w′
(2.38)
W′=f31u′+f32v′+f33w′
Substituting velocity corrections from (2.37) into (2.38) and dropping derivatives of p′ those are not in the direction of U, V and W respectively yield f112 f122 f132 U = U − u + v + w pξ′ = U * − C u pξ′ AP AP AP *
f2 f2 f2 V = V * − 21u + 22v + 23w pη′ = V * − C v pη′ AP AP AP
(2.39)
f2 f2 f2 W = W * − 31u + 32v + 33w pζ′ = W * − C w pζ′ AP AP AP
Substituting the contravariant velocity components into the continuity equation (2.30) and introducing ( pξ′ )e = pE′ − pP′ , etc., we obtain the pressure correction equation: AP p′P = AE p′E + AW pW′ + AN p′N + AS pS′ + AT pT′ + AB p′B + S P′
where
(2.40)
20 AE = ( ρ C u ) AN = ( ρ C AT = ( ρ C
v
u
AW = ( ρ C u )
e
) )
AS = ( ρ C
n
AB = ( ρ C
t
v
) )
w
s
u
(2.41)
b
AP = AE + AW + AN + AS + AT + AB e
n
t
S p′ = − ρU * − ρV * − ρW * w s b
It should be noted that S p′ , represents mass imbalance in the control volume. Once the pressure correction field is known, the correct pressure, velocity field and scalar quantities can be obtained. The process repeats itself until convergence using the underrelaxation. 2.5
Convection Scheme
The central difference operator is employed for the pressure and diffusive fluxes. The second-order-central difference operator for the first derivative convection terms gives rise to wiggles when the cell-Peclet number (Reynolds number based on the dimensions of the largest grid cell) is greater than some critical value (Shyy [39]). For the one-dimensional Burgers equation, the critical value is 2. An enormous variety of alternative discretization schemes have been proposed for velocities. To overcome this problem, the first-order upwind scheme has been widely used. Although it is a safe choice, the excessive numerical dissipation inherent to this scheme has necessitated the search for more accurate schemes. In the present work the second-order upwind scheme (Thakur and Shyy [33]) has been adopted as a basis. The discretized form of the convection terms in the u-, v-, w-momentum equations can be written from the left hand side of equation (2.32) as
[ ρUφ ]ew + [ ρV φ ]ns + [ ρWφ ]bt ≡ Feφe − Fwφw + Fnφn − Fsφs + Ftφt − Fbφb
(2.42)
21 where φ=u, v or w. The mass fluxes at the control volume interfaces, F’s, are estimated by a linear interpolation between the two neighbors. For example for a uniformly spaced grid we have Fe =
1 1 ( FE + FP ) = ( ρU ) E + ( ρU ) P 2 2
(2.43)
The variable φ at the control volume interface can be determined using an interpolation or extrapolation of the values of φ from one or several grid points. The choice of the geometric shape of the interpolation/extrapolation profile dictates the type and the local order of accuracy of the convection scheme.
Figure 2.2. Schematic of the linear extrapolation employed by the second-order upwind convection scheme. Adapted from Thakur and Wright [27]. Figure 2.2 shows schematic of the linear extrapolation employed by the secondorder upwind convection scheme depending on the flux direction. This can be summarized as follows: 1 1 3 3 Feφe = φP − φW $ Fe ,0% − φ E − φ EE $ − Fe , 0% 2 2 2 2
where the operator $ a, b% denotes the larger of a and b.
(2.44)
22 In the present implementation the first-order upwind scheme is used as a basis and terms that are different between the two schemes are assigned to the source term SC, which is treated explicitly (Thakur and Shyy [33]). 2.6
Multiblock Algorithm
Complex physical domains occur in many fluid flow problems, making implementation of single block structured grids either inefficient or impossible. A structured grid is implemented in this study (Rai [39], Wright and Shyy [41], Shyy [32], Shyy et al. [31]). Although the structured multiblock approach is not as flexible as the unstructured grid [42], [43] for very complex domains, it is simple to implement, it allows employing well-established algorithms and it provides accurate resolution of nearwall behavior in turbulent flows. Composite structured grids can be classified as either patched grids (also known as abutted grids) or overlaid grids. Neighboring subgrids, in general, can overlap and do not match exactly at the internal boundaries for overlaid grid. In the current multiblock computational method the patched approach is implemented, requiring grid lines to be continuous between blocks. Although it is more restrictive, it offers straightforward information transfer across the internal boundaries with the flux conservation. The primary issue of concern for a multiblock solution technique is the transfer of information in the vicinity of the internal block boundaries. In order to facilitate the exchange of information between the two blocks, two layers of extended control volumes in each internal boundary are constructed. These control volumes extend two layers deep into the neighboring block. Taking value from the neighboring blocks the extended lines are treated as part of the parent block (Wright and Shyy [41]).
23 2.7 2.7.1
Turbulence Modeling
Transport Equations for k and ε
Practical turbulent calculations at high Reynolds numbers are often performed using the Navier-Stokes equations in the averaged form (equations (2.6) - (2.8)). As already mentioned, the time averaging process gives rise to new unknowns, including socalled Reynolds stress tensor terms, − ρ ui′′u′′j . It is a symmetric tensor, having six independent components. Introduction of the new unknowns requires introduction of new governing equations to close the system. One of the earliest proposals for modeling Reynolds stresses is the Bousinesq’s eddy-viscosity concept, which assumes the tensor to be a scalar function proportional to the mean flow velocity gradients: ∂u" ∂u" 2 ∂u" 2 − ρ ui′′u′′j = µt i + j − µt l ⋅ δ ij − ρ k ⋅ δ ij ∂x 3 j ∂xi 3 ∂xl
(2.45)
where µt is the eddy viscosity to be modeled. The eddy viscosity can be written as
µt = ρ
l02 t0
(2.46)
where l0 is the turbulent length scale and t0 is the turbulent time scale quantities, which can be varying significantly for different flows. Jones and Lauder [44] introduce a two-equation turbulence model for two turbulent quantities directly related to the turbulent length and time scales. In the most widely used k-ε turbulence model, the turbulent length and time scales are defined based on the
turbulent kinetic energy k and dissipation rate ε, defined as follows: k=
1 ui′′ui′′ 2
ρε = µ
∂ui′′ ∂u′′j ∂x j ∂xi
(2.47)
24 The eddy viscosity is given by k2 µt = ρ C µ f µ ε"
(2.48)
where Cµ and fµ are empirical proportionality constants to be defined later and ε" = ε − D is the difference between the dissipation rate and a near-wall function used for near-wall damping in low Reynolds number models. Transport equations for k and ε can be obtained by the following procedure (Wilcox [45]). Let N(ui) denote the Navier-Stokes operator, then forming the following ensemble averages yields
(
)
2µ
∂ui′′ ∂ ( N (ui ) ) = 0 ∂x j ∂x j
trace ui′′N (u j ) + u′′j N (ui ) = 0 (2.49)
Introduction of modeling assumptions yields ∂ ∂ ∂ ρk + ( ρ ui k ) = ∂t ∂xi ∂xi
µt µ + σk
∂ ∂ ∂ µt ρε" + ( ρ uiε" ) = µ + σε ∂t ∂xi ∂xi
∂k + Π − ρε ∂xi
(2.50)
∂ε" ε" ε" 2 ρ + Π − +E C C ε1 ε2 k k ∂xi
(2.51)
where E is a near-wall function used by low Reynolds number models and Π is the generation (production) of k, which is the function of the mean rate of deformation Eij: ∂u ∂u ∂u Π = − ρ ui′′u′′j ⋅ Eij = 2 µt Eij Eij = µt i + j ⋅ i ∂x j ∂xi ∂x j
(2.52)
The original k-ε model, formulated by Launder and Spalding [17], has five closing constants, as shown in Table 2.1. Numerous modifications have been proposed to the original k-ε models. Non-equilibrium model by Shyy et al. [31] and the low Reynolds model by Chien [47] are discussed next.
25 Table 2.1. Constants and functions for k-ε models Model D Cε1 Cµ Original [17] 0 .09 1.44 Nonequilibrium 0 .09 1.15+ Shyy et al. [31] 0.25Π/(ρε) Chochua et al. [46] ≤1.8 2 .09 1.35 Low Reynolds 2νk/y Chien [47] Model Original [17]
fµ 1.0
Nonequilibrium Shyy et al. [31]
−3.4 exp 2 (1 + Re t 50) 1-exp(-.0115y+)
Low Reynolds Chien [47]
*
Cε2 1.92 1.45+ 0.45Π/(ρε) ≤2.375 1.8
σk
σε
1.0
1.3
1.0 .8927
1.3 1.15
fε1 1.0
fε2 1.0
E 0
1.0
1.0
0
1.0
1-.22exp[-(Ret/6)2]
−2 µ (ε" y ) × 2
(
exp −0.5 y
+
)
* Ret=ρk2/µε
2.7.2
Nonequilibrium k-ε Model
Shyy et al. [31] propose nonequilibrium modifications to the original k-ε turbulence model, for flows where the assumption of equilibrium between rates of production and dissipation of the turbulent kinetic energy does not valid. Flow equilibrium does not exist in important fluid flows such as flows with rotation, adverse pressure gradient, recirculation, and large streamline curvature. The main idea of the modifications is adding a time scale in the ε equation by introducing functional forms for Cε1 and Cε2 as function of imbalance between the production and dissipation of k. This model also absorbs the low Reynolds number effect with modified function for Cµ. Table 2.1 summarizes the modeling constants. The constants are tuned up by considering performance of backward facing step computations. Chochua et al. [46] show that jet-incrossflow computations produce strong imbalance between the production and
26 dissipation of k, leading to computational instabilities. Adding cutoff values helps improve computational stability. 2.7.3
Wall Functions and Low Reynolds Number Models
Owing to the presence of the solid wall, the flow behavior and turbulence structure are considerably different from the free turbulent flow and it is usually based on the correlated experimental data in the turbulence modeling. If we consider the Reynolds number based on a distance y away from the wall and friction velocity u* = τ w ρ , i.e. Rey =y+=(ρu*Yp)/µ, it characterizes the importance of inertia forces, which dominate in the flow far away from the wall. In the fully turbulent region wall shear stresses are defined by turbulent stresses. As y is deceased to zero, however, the Reynolds number based on y also decreases to zero, and at some location it will be order of 1. At this distance from the wall and closer the viscous forces will be equal in order of magnitude to the inertia forces or higher. In flows along solid boundaries there is usually a substantial region of inertia-dominated flow far away from the wall and a thin layer within which viscous effects are important. For a fully developed turbulent flow, near a no-slip smooth wall, the normalized near wall tangential velocity, assuming a two-layer structure (viscous sublayer followed by the log layer), can be written as follows: y+≤11.63
u+=y+ u+ =
1
κ
log( y + ) + 5.15 y + > 11.63
where u+ =
ut u*
u* =
The von Karman constant κ is usually taken as 0.41.
τw ρ
(2.53) (2.54)
27 Figure 2.3 shows the agreement between theoretical equations (2.53) and (2.54) in their respective areas of validity and experimental results for a typical flow.
Figure 2.3. Zones in the turbulent boundary layer for a typical flow over a smooth flat plate. Adapted from Tannehill et al. [30]. The original k-ε model is based on assumption of the local equilibrium between the production and dissipation of turbulent kinetic energy:
− ρ ui′′u′′j =
ρε
(2.55)
Eij
Form the assumption of the constant stress layer near the wall, dominated by Reynolds stresses the turbulent kinetic energy can be related to the wall shear stress as
τ w = − ρ ui′′u′′j = µt Eij = ρ Cµ ⇒
k2
ε
Eij = ρ 2Cµ
k2 k2 = ρ 2Cµ τw − ρ ui′′u′′j
(2.56)
τ w = ρ Cµ k
Combining equations (2.54) and (2.56) the expression for the near-wall shear stress yields
ρκ Cµ k 1
τw =
4
1
2
ln ( y + ) + 5.15
ut
(2.57)
28 Patel et al. [48], reviewing modifications to the k-ε model at low Reynolds numbers, shows that k + ( = k u*2 ) starting from k+=0 at the wall becomes maximum at around y+=15. In the interval 60 < y+ < 150, k+ becomes almost constant and takes value of around 3.3 (it corresponds to Cµ =
1 k+
2
=0.09 in equation (2.56)). The normalized
+ −u′v′ Reynolds stresses component −u′v′ = 2 has zero value at the wall and u*
monotonically increases towards the value of 1.0, which is one of the assumptions accumulated in the wall functions. The ratio between production and dissipation steeply µε drops, as y+ becomes less then 15. The dissipation rate ε+ = has a finite nonzero 4 ρ u*
value at the wall. The wall functions are used to produce accurate results without involving relatively fine grid resolution in the wall region, which would be needed in order to capture steep variations in the flow profiles. However they are limited by relatively high values of y+. If this requirement is violated then solutions may contain strong discrepancies. Many suggestions have been made for the extension of the turbulence closure models to enable their use at low Reynolds numbers. Low-Reynolds number models approximate the governing equations all the way to the wall and thus obviate the need to make any assumptions about the nature of turbulence or the velocity profile near solid walls. Wall damping coefficients are commonly used to ensure that viscous stresses take over from turbulent Reynolds stresses at low Reynolds numbers at viscous sub-layers adjacent to the wall. Table 2.1 presents closure coefficients and damping functions for low Reynolds number model by Chien [47].
29 2.7.4 Performance of the Wall Functions and Low Reynolds Number Models
In this subsection effect of the wall treatment on computational results is discussed. 2-D computations of airflow are performed based on the experimental setup for flows between two smooth plates of the test case 2 by Ha and Childs [19]. Figure 2.4 shows the computational domain representing a 2-D channel. The inlet pressure pin=12.4 bar, the channel clearance H=0.38 mm and channel length L=76 mm, making the channel aspect ratio L/H=200. Experiments are conducted for the Reynolds number range, Re, (based on the hydraulic diameter 2H) from 8,090 to 8,330, corresponding to the transitional flow regime in pipes. Table 2.2 shows summary of the experimental setup. The well-known analytical solution for laminar flows in a channel yields f = 24
(2.58)
Re
For turbulent flows, an often-used Blasius-type friction factor formula for fully turbulent flow with smooth walls has the following form (e.g. Childs [49]): f = 0.079(Re) −0.25
Table 2.2. Summary of the computed test cases. Turbulence # Case Inlet BC/grid Model description size/ ω/Inlet swirl Full/12,663 Wall 1 Flow between 0.0/0.0 Functions flat plates. Uniform mesh. Test 2 (Ha and Childs [19]) Full/12,663 Low Re 2 Flow between 0.0/0.0 flat plates. Strong nearwall refinement. Test 2 (Ha and Childs [19])
(2.59)
Reynolds number 8,090 23,000 38,100 53,030 68,330 8,090 23,000 38,100 53,030 68,330
Axial Mach number 0.036 0.103 0.177 0.258 0.356 0.036 0.103 0.177 0.258 0.356
Inlet Pressure, (bar) 12.4
12.4
30 The domain is discretized using a 201×21 mesh. The mesh is uniform along the length of the channel. Figure 2.4 shows a portion in length of the computational domain with two types of transverse clustering.
Figure 2.4. Channel domain dimensions. Mesh with near-wall refinement H/dy1=100 (left) and uniform mesh H/dy1=20 (right). Both uniform and non-uniform meshes are considered in the transverse, y-, direction. For the non-uniform mesh, a grid-clustering law is chosen from the family of general stretching transformations as given by Tannehill [30]: ( 2 y −1)
− β +1 ( β + 1) ( β + 1) ( β − 1) y=H 2 − y ( 1) 2 1 + ( β + 1) ( β − 1)
{
}
(2.60)
where y is the nondimensional distance from the wall of the channel for the uniformly distributed grid. The stretching parameter β is related (approximately) to the nondimensional boundary-layer thickness (δ/h) by
δ β = 1 − H
−1 2
0
