- Email: [email protected]
Streamwise computation of three-dimensional parabolic flows
COMPUTER
METHODS
NORTH-HOLLAND
IN APPLIED
PUBLISHING
STREAMWISE
MECHANICS
COMPUTATION
of Mechanical
Wichita State Unitlersity.
Received
A
new
the
approach
computation
dependent difference mass.
of
velocity
three-dimensional
plane,
of the chosen
the chosen Results
obtained streamlines;
of the calculations.
data for flow through
rectangular
19 November
1981
received
10 March
three-dimensional
in the present
equations
manuscript
parabolic
along a set of streamlines.
variables
cross-stream around
to calculate
by calculating
36 (19X3) 7 I-93
flows
approach
are
the
The
Euler’s
the streamline
Wichita. KS h72Ob’. U.S.A.
1982
flows
is presented.
The dependent
three
velocity
are the streamwise
streamlines.
by applying
PARABOLIC
S. GREYWALL
Engineering,
Revised
computed
ENGINEERING
OF THREE-DIMENSIONAL FLOWS
Mahesh Department
AND
COMPANY
velocity velocity
theorem
are
calculated
approach.
ducts;
is satisfactory.
the agreement
flow field is used in
In contrast,
and the two coordinates.
streamwise
coordinates
The
commonly
components.
momentum
based on the present
variables
is calculated
from
to streamtubes from
are compared
the in the
the finite
constructed
the conservation
of
with the experimental
1. Introduction Viscous flows have been calculated for a long time by separating the flow into a boundary layer and a potential flow region. The boundary layer then is calculated by using integral methods and the potential flow by using inviscid flow equations. This approach does not require large computers, and even can be carried out by use of hand calculations. This approach, however, has its limitations. For example, in calculating flow through nozzles, channels, and diffusers this approach is limited to the region near the entrance. As the boundary layer thickens, the flow cannot be meaningfully separated into a potential flow region and a boundary layer region. Further, this approach cannot describe the corner effects, such as the secondary flows, or adequately describe flows in which the body forces are nonuniform in the cross-stream plane (such as magnetohydrodynamic channel flows). In the past two decades, with the advent of high speed computers, many numerical methods to compute velocity and temperature fields without separating the flow into a boundary layer and a potential flow region have been developed. In this paper, a new approach to calculate three-dimensional compressible viscous flows is presented. The approach is an extension of the author’s previous work [l] to three-dimensional flows. The method, as presented in this paper, is for ‘parabolic’ flows through rectangular ducts, and its possible extensions to other ‘Parabolic’ or ‘boundary layer’ flows types of three-dimensional flows is under investigation. are characterized by the existence of a predominant how direction along which downstream conditions have a negligible influence on the upstream conditions. This assumption allows use of a marching integration procedure. A lucid discussion of the parabolic flow assumption is given by Caretto et al. [2], and Patankar and Spalding [3]. 00457825/83/0000-0000/$03.00
@ 1983 North-Holland
72
M.S. Greywall,
Stream wise computation
of parabolic flows
Methods developed so far to compute three-dimensional viscous flows are, among others, those given in (2-101. These methods differ in their finite difference approximations, and whether the flow is assumed to be parabolic or not. From among these methods, those proposed for calculation of confined flows differ further in the procedure used to calculate the pressure field. These methods are discussed here briefly. Included in the discussion is the method for the time and convective terms, and Dufort-Frankel for the ditiusion terms. The transient viscous flows, it has some features in common with the methods developed later for of Chorin [4] calculates computation of three-dimensional parabolic flows. The method incompressible elliptic flows. Finite difference equations are obtained by using the leapfrog method for the time and convective terms, and Dufort-Frankel for the diffusion terms. The pressure field is calculated by coupling the incompressible continuity equation to the momentum equations via an artificial variable density and an artificial equation of state. Miller [S], on the basis of the work of Chorin, computed steady incompressible three-dimensional elliptic flows through rectangular ducts. In the method of Harlow and Welsh [ 111 the finite difference equations are obtained by using forward in time and center in space differencing. The pressure field is calculated from an elliptic equation obtained from the continuity equation and the divergence of the momentum equations. The method uses a ‘staggered’ grid in which the velocities are defined at the cell boundaries and pressure at the cell centers. A similar type of staggered grid is used in the cross-stream plane, later, in [2,X 6,8] for the computation of three-dimensional flows. In these methods the cross-stream velocities are specified at the cell boundaries and the axial velocity, along with all the other variables. at the cell centers. Methods [2,3,6, 71 calculate three-dimensional confined parabolic flows. In these methods the pressure field in the cross-stream plane is calculated from an elliptic pressure equation obtained from the continuity equation and the cross-stream momentum equations. The axial. along the duct axis, pressure gradient is calculated to insure conservation of mass along the duct axis. In these methods Caretto et al. [2], and Patankar and Spalding 1.31 obtain finite difference equations by using the control volume approach; Emery et al. [6] use DufortFrankel for the axial velocity momentum equation and implicit differencing for the crossstream velocity momentum equations; Roberts and Forrester [7] use one-sided upstream differencing in the flow direction and second order differencing in the other two directions. Pratap and Spalding [8] extend the method of Patankar and Spalding to ‘partially parabolic’ flows. Methods of Nash [‘i’] and Wang [lo] are for external (pressure field specified) threedimensional parabolic flows: Nash uses forward explicit differencing, and Wang uses a Crank-Nicolson type of differencing scheme. All the methods described in the preceding paragraph use the Eulerian system. The flow field is computed by calculating the velocity components u,, u,, and u, along a set of spatial grid points fixed in advance. In contrast, in the present approach, the flow field is computed by calculating the streamwise (not axial) velocity, u, along a set of chosen streamlines. Streamlines are identified by a pair of indices, (i. j), such that the coordinates of the streamline (i, j) in the cross-stream plane are Yi and Z,. Coordinates of the streamlines in the cross-stream plane are calculated as part of the flow computations. The dependent variables in the present approach are Ui,iTthe streamwise velocities, and (Yi, Zj). the coordinates of the streamlines. To calculate the streamwise velocities streamtubes are constructed around the individual streamlines such that the sum total of all the streamtubes fills the duct. One such streamtube. constructed around the streamline (i. j). is shown in Fig. 1. Finite difference equations for ld,~,
M.S. Greywall.
Streamwise
computation
of parabolic flows
An adjacent
---
X
Fig. 1. Streamtube
73
streamline
Streamline
(i,j)
Duct
a xis
X+&X
(i, j) constructed
around
the streamline
(i j).
and hi.,, velocity and enthalpy along the streamline (j, j), are now obtained by applying Euler’s momentum theorem and the first law of thermodynamics directly to the flow in the streamtube between x and x + Ax. After the streamwise velocities and densities have been calculated at x + Ax. the streamline coordinates at x + Ax are cafculated next by using mass conservation. From (Yi, 2’) and (Y’, ZT), coordinates of the streamline (i j) at x and x + Ax, one can calculate the slope of the streamline (i, j), and then using this slope decompose u,., into its components u;i, uz’, and ufi.
2. Finite difference equations 2 1. Background In most numerical approaches to the computation of the viscous flows one starts with the partial differential equations-Navier-Stokes equations, the energy equation, and the continuity equation-obtained from the conservation principles. These equations, either in their primitive forms or in transformed forms, are then differenced to obtain the needed finite difference equations. In contrast, in the present approach the flow in the duct is partitioned into a finite number of streamtubes and the finite difference equations are then obtained by applying the conservation principles directly to the flow through the individual streamtubes. Consider the flow through a rectangular duct with its axis along the x-direction. To define the various streamtubes lines Y2, Y3, . . , , Yiq . . . , YI are drawn parallel to the z-axis (left duct wall), and lines Z,, Z,, . . . , Z”, . . . , Z, parallel to the y-axis (lower duct wall), as shown in Fig. 2. Lines Y1 and 2, (not shown in Fig, 2) are drawn through the middle of the lower wall and the middle of the left wall, respectively. The lower and left walls are also referred to as Z, and Y, lines. Line Yi+l/z is drawn in the middle of Yi and Y,,,, lines Yi_1,2. Zj+j/zr and Zj-l,z are drawn similarly. Various streamtubes, called streams for short, are now defined as follows. Stream (2,2): flow bounded by the left wall, the lower wall, Z2+1,7_,and Y2+1,7. Stream (i. 2) i = 3, . . . , I: flow bounded by the lower wall, Z2+1,2r YI:-l,z, and Yi+i/z. Stream (2, j), j = 3, . . . , J: flow bounded by the left wall, Y2+l,z$ Zj-x,z, and Zj+r/z. Stream (i, j), i = 3. . . . , I, j = 3, . . . , J: flow bounded by Yi_-1,2. Y<+j/z, ZjP,,z, and .Zi+r,2. Streamline (i, j) is defined as the streamline passing through the intersection of Yi and Zi.
M.S. Greywall,
74 z
r2 ly 2+%
’ I
yi_‘/
’ I
’ I
_L_L__L____L___Zj, I I I
I
I /q/j, ‘,// : *;‘j:
I , I
i,j ;;g+ ;_ ‘,/ ,,,; ij /, ,,///
L _: I
yi+‘/t
(i,j) l&$~;’
:
I --
Fig. 2. Definition
--zj 2.
J-S
I _ _ L _=_2+: z2 ,:
i.1
1.1
LI I
Tgg;;“: $ (1.2)::
‘j+l/,
I
., ‘1
I
L ,
1
I
I Y”“’ w;;, q,: , I
_L_L__L____L___~3 I I I
ggg ,(2#2)
of parabolic flows
yi+l
yi T3
I
Stream wise computation
of the various
I
Y
streams.
+.I -
%
Fig. 3. Control
Yi+K
volume for the stream
(i, j).
Other nomenclature is as follows: !Py,! is the mass flow rate through stream (i, j). The mass flow rate !Pc is partitioned further into four parts !PlT, q;,+, V’i,, and q’:; ; mass flow rates through the upper-right, upper-left, lower-left, and lower-right quadrants, respectively (see Fig. 2).
AC is the cross-sectional area of the stream (i, j). To apply conservation principles control volumes are constructed for each stream distance Ax along x. The control volume for the stream (i j) is shown in Fig. 3.
using a
As1, A?: surface areas of the control volume normal to the y- and z-axis respectively. ui,,, Ti,,, hi,j, pi,,, ciJ: velocity, temperature, enthalpy, density, and constant pressure specific heat along the streamline (i, j) at x. U;, etc.: variables along the streamline (i, j) at x + Ax. Ui,j: average of Ui,j and u ;. ZA&, ufj: of iii,j the X, y, and z directions. Yi,Zj: coordinates of the streamline (i, j) in the cross-section plane at x. Y:, Zi+: coordinates of the streamline (i, j) at x + Ax. Si,j s
12
TU r,j.
+_1
t2 = 2lA i,j . ',I I p i+ 1/2,j, p f- 1/2,j, p
s
{j+ ,127 p :j-,121 effective viscosities, including the turbulent contributions, for the momentum transfer across the surfaces of the control volume, as shown in Fig. 3, evaluated using the information at x. pfzIIz,j, etc.: effective viscosities evaluated using the information at x + Ax. K!+112,j, etc.: effective thermal conductivities, defined similarly as the effective viscosities.
M.S. Greywall, Stream wise computation of parabolic flows
V’s and H’s with superscript
Corresponding
+ are defined
in terms of p+, K’,
Y’, and Z’.
In estimating momentum and energy fluxes through a streamtube u and h are assumed, in general, to vary linearly between the neighboring streamlines. Between the walls. however, and the adjacent streamlines (streamlines with i or j equal to 2 in Fig. 2). a nonlinear variation of u and 11 is allowed. Thus, the mass. momentum, and energy flux through the streams adjacent to the walls are expressed using coefficients y, (Y, p, and 6 defined as follows (see Appendix A for further discussion):
The following relations are for the lower quadrants of the streams along the lower duct wall. Formulae for the lower-left quadrants are obtained from these relations by substituting !P;; for ?&
and for the lower-right
quadrants
by substituting
U = $(3Z(j.2+ Ili+l,z) 1
AY = +(r;+, - Y), 6 = $3h,,,
+ hi+t.,)
PT; for ?P,‘,.?,
and
&? = &$
+ iuf+,.J, Z2
Z2 Pi.2 = A Y
AY
I
pudz
=
y:@AYZ,
,
AY -3
ZZ
0
pUUdZ~CXfTPi,~G, I0
pu tu2dz = S; ‘&p,
Ay
f0
puhdz
Coefficients 77, a;, ST, and py, that take into account the left wall, are defined similarly.
= /?‘;‘iQh.
the non-linear
variations
normal
to
76
MS.
Greywall,
Strearnwise
computation
of parabolic flows
In evaluating the various integrals for the near-wall correction factors, p has been set equal to some constant p. For incompressible flows p is, of course, just equal to p. In computing two-dimensional compressible flows, p set equal to the average of p between the wal and the adjacent streamline gives satisfactory answers. For the three-dimensional compressible case the appropriate value to be used for /s is still under consideration. 2.2. Finite difference equations for Idi,, Applying
Euler’s
momentum
Momentum = Pressure
theorem
to stream
flux outlat r+A.r- Momentum forces + Viscous
(i, j) between
x and x + Ax one obtains,
flux ini at x =
forces + fi,j(Body
forces) .
denotes the quantities in the preceding bracket [ ] with U+ by replaced by 24. This symbol is used similarly later on. In calculating pressure and viscous forces, it is assumed that the streamline inclinations with with u.*. Thus, in respect to the duct axis are small; that is, U” and uz are small compared calculating pressure force A;j+ the cross-sectional area of the streamtube (i, j) projected onto the Y-Z-plane, is used. In calculating viscous forces, the velocity gradient is based on An, (see Fig. l), the distance between the neighboring streamlines projected onto the Y-Z-plane. Pressure force on stream (i, j) is estimated as 141
- DP A;,
(2.3)
where DP is the pressure difference between x and x + Ax and is assumed to be the same for all the streams. The assumption that DP is the same for all the streams restricts the analysis to straight ducts without any abrupt change in area along the duct. The viscous force on a surface between x and x + Ax is approximated by the average of the viscous forces evaluated by using the effective viscosities and the velocity gradients at x and x + Ax. One thus obtains for the viscous forces on the stream (i, j) (see Appendix A), Viscous
Forces = i[ V’+ t+l,j(l4T+l., - 44:,)+ -
V:j+(U:j
-
14:
1.~)~
V::+I(UZ,+ Vi,+(Z4lj
-
I - Z4lj) L4,tl- I)]
+ +[
V’.
U+ +
V, U] .
(2.4)
M.S. Greywall,
Streamwise
computation
of parabolic flows
77
For streams along the walls, momentum fluxes are expressed using the non-linear correction factors introduced in Section 2.1. Viscous forces exerted by the lower and the left wall on the adjacent streams are, respectively, set equal to V&U~,~and V;,jU,,z. The calculation of Vi2 and Vi,, is discussed in Appendix A. After the various estimates in (2.1) have been substituted, all the terms involving Ucj are collected on the left-hand side, and the resulting equation is organized in the form
Coefficients A, B, C and D for the various streams are given next. 1 Z, are i times the ly’s defined in Section Note. In equations that follow all the q’s, except !J’_, 2.1, and the V’s are $ times the V’s defined earlier. Coefficients for stream (2.2):
Coefficients for stream (i, 2) i > 2:
Coefficients for stream (2. j), i > 2:
78
M.S. Greywall.
Coefficients
for stream
-7.3. Finite difference
Applying
computation
of
flows
(i. j), i and j > 2.
equations
for /Ii,,
first law of thermodynamics
(i, j) between
to the stream
x and x + Ax, one obtains,
Energy flux outlat x+A.Y - Energy flux inlat X= = Heat added by conduction + Work done by the viscous forces + Q,,(Internal heat generation)+ wi,j(Work done by the body forces).
(2.6)
For the streams (i, j) with i and j > 2, various terms in (2.6) are estimated as follows. The energy flux terms are estimated in the same way as the momentum flux terms and can be obtained by substituting (h + s) for u in (2.2). The heat conducted across the top surface of the control volume (see Fig. 3) is estimated as
Work done by the viscous forces
$[V<~+,(U~j+l
-
14,:I)(i(Ut+l
=f[V{,‘+,(.S~,+[ -
aS:j)]
on the top surface +
Ulj))l
+t[V+S++
+
is taken to be
$[V+Ll++
VU] =
VS] .
(2.X)
Work done by the viscous forces and the heat conducted across the other three surfaces of the control volume are estimated similarly. After the various estimates in (2.6) have been substituted, all the terms involving h:, are collected on the left-hand side and the resulting equation is organized in the form I
+
-A,,jh;+i,j-
I
+
+
C,.jh-1.1 + Bi.,h ,,,+I-A{jh,T,+l-
CfjjlZ,p, =Dt.,.
Coefficients A. B, C and D for the various streams are given next. Note. In equations that follow all the P’s, except W.5, are $ times Section 2.1, and the H’s are $ times the H’s defined earlier. Coefficients
for stream
(2.2):
(2.9)
the
q’s
defined
in
MS.
Greywall,
Streamwise
of parabolic flows
computation
79
02.2 = Q2.2 + W,,, + h,,z( Yf;,t + #of !P;.; + El:,,) + h,,(2!P;.;
+ 3pS !P,; + Pf PT., + 3p5 P’t; - H:,, - H;., - Hi., - Hi,,)
+ hz..J ‘P’tt + PP Y’,: + Hi.3) + s&Y’;.: + s2,42ly;.;
+ 36; !P,; + s; ?K,
+ 6; W;,s + V:.,)
+ 36; Ye;
- vi.2 - v;,, - v;,, - VS.,)
+ s2,3(!I’;.; + 6; T’;;
+ Vi,,) + h ,.,H:,, + h,l H;,z
+ h,,,H:.:; + h,,H::
+ s:,,(-
+ 3c5;*;,:
+ s;,2(2K;
!P;.; - S; P;,, + V:r’,)
+ s; YP_2 + 3s; w;
+ Sl.3(- !P;.; - sp !P;,; + v:.;,
+ v:;
+ vg
+ v;r+, + VC)
.
Coefficients for stream (i, 2), i > 3: A:2 = - !l’;; - Pf W;; + Hf;,.2 , C:,= -!&-/3f!P;;+ Bi,2=2~~~+3Pf~(i17_+Hf=,,2+H:_~+H:f+H:2+, A? = - IJY(?+ + H-!’ c:, = 0 . I.3 .
H:;,
D;;f = Qi,2 r;’ Wi,2
+Si+l,Z(!P~~+61*~~+
Vf+l,Z)+Si~,,Z(W;~+~flVi.2+
Vf.2)
+ Si,2(2T:z + 36: Py; - Vf+ 1,’- Vf3 - V{2 - V:2) + Si,3(q?:
+ V$) + h,lH:2
+ h,lH{i
+ s:+,,J(- !P’1+2 - Sf !+Jq-y+ VfIl,,) + s:_,.2(- !P’;; + s;22w?+I.2+ 36; p- 1.2+ Vf’,+ 1.2+ v:.; + v:; +s;3-
!P.:I;+
sf 9,
+ v::
+ v:;
v;:).
Coefficients for stream (2, j), j > 2: A’.=
c:, = 0 ,
-~V+‘?+ff’+.
B:; z 21$!+
3&+
+ Hi; + H;,;+l + H:; + H:,; , c;~~ = - ~7 Pi,; - !P;,; + H:; A;,j = - ‘I’;.; - /3; ‘Pi,; + Hi,;+ 1 , Di,j = Q1.j + W2.j
,
+ h,j ( ‘JT,: + Hi,) + hz.j(2Pz,Y + 3p7 *y.J! - H:,j - HSI,j+l- H:j - H:j) + hz,j+,( !$‘z.J+ ~7 Pi,f + H;,j+l) + hz.j-l(pyqY.Y + s3.j( qi.7 + Vi.j) + &,j(2WZJ + 363 YP2.ys2,j+*(Tl,f
+
h1,jHi.j + h1,jHi.T + Sl,j(- pZ,y + V:,$)
- &j(2Pl,7 +
S2,j+l(-
+
67
P2.f
Vs,j+,) +
+
+
S2,j-1(6;
+ 3S7!P’;q + V$,J + V;,T+l + ql,f
-
SyP2.J
+
V5,j
%,y
Vi,j
+
+ %.T + Hij) -
VT&+1
+
Ti,i
-
+
Kj
-
vCj)
vii)
vCj)
V:,T+l,>+ S2,j-1(-8~pS.j
- pii
+ v:,T) ’
80
M. S. Greywull,
Coefficients
for stream
Stream
wise cornputution
of parabolic
flows
(i, j), i and j > 3: c’
A:, = - !P:;‘+ H!:,,, ,
= - I)-!’ + H?t r., . 1.1 I.,
+ H:;+, + H:J + H::
B,,, = ‘Py; + Hf:,,,
A-’ = - I$% + H!?r.,+l 7 I., I.,
c:,
3
= - ‘P::; + H;I; ,
4.j = Qt., + W, + h,+,,,(v:;‘+ + h,,(‘J$-
IX+,.,)+
h,- &K;‘+
H:,)
HJ+,,,-H:,+,-H:,-H:j)
+ h,.,+,(T:IT + H:,+l)+
h+I(CIJ
+ HL)
+ s;+,.,(!Plt:)+ Vf+,.,)+ him.,+ +
S,j(
PI!:’
-
Vf+
1.1
-
K,)
V:j+ I - vt] - v:,)
+ S,.,+l(T:tT + Vij+l) + Si.,pI( %17 + v:]) + ST+I,,(- P,f:’ + VfI1.j) + ST-~ I.j(- qt,” + vt:) - s:,pP:y+
v;:,.,
+s:,+,(-P:;;+
+ V:: + V$)
v:;,,+s:.,.
of u and
2.4. Decomposition
+ v::+,
the secondary
,(-El,+
Vii).
flows
Computation of secondary flows. defined as the velocity components normal to the main flow, is discussed in this section. After y+ and Zt have been calculated, the streamwise velocity iii,j, the average of U,,j and Cl:j. then can be decomposed into II;,, Iit+ and U;j, by using the formulae,
*
11i.1
=
ii,,, where
AX Ar
u:,_AY ’
Ax is the integration AY=
Y’-y,,
iii.1 -
and
Ar
uf,_AZ -_C,,j
(2.10)
Ar
step, AZ=ZT-Zj
and
Ar’=AXz+AY’+AZ’.
The calculation of Y’ and ZT as discussed in Appendix A is based on the conservation of mass alone. Velocity components uy and u’, obtained from Y’ and ZT thus calculated, are the secondary flows associated with the flow development in the ducts-may the flow be laminar or turbulent. In turbulent flows through rectangular ducts, besides the primary shear stresses which are in the streamwise direction, there are also shear stresses in the plane normal to the streamwise direction (see, for example, [12]). Force due to these stresses will produce bending of the streamlines in the z-x and z-y-planes, and thus displace Y’ and Zt beyond that displaced by the conservation of mass. For example, force per unit area AT_ due to the cross-stream shear stress r,, distribution, will bend the streamline (a j) according to the formula (2.11)
M.S. Greywall,
Stream wise computation
of parabolic flouts
81
where R,, is the radius of curvature in the z-x-plane. With Rg, (approximated as a straight line for small Ax) known one can calculate displacement of ZT (and similarly of Y’) caused by the cross-stream turbulent shear stresses. Contributions to II.’ and u”, as calculated from (2. lo), by this additional displacement of the streamlines in the cross-stream plane represents, as viewed from the present computational approach. the secondary flows associated with turbulent flow through noncircular ducts. It should be noted that the turbulent shear stresses in the cross-stream plane do not enter into the finite difference equations for the streamwise velocity presented in Section 2.2. These equations are valid for laminar and turbulent flows as long as the streamwise shear stresses, appearing in the equations via V’s, are modelled appropriately.
3. Solution technique A methodology to solve the momentum and energy finite difference equations given in Section 2 is presented in this section. This is the methodology used for the sample flow computations presented in Section 4. Equations (2.5) and (2.9) are of the general form - A :,jX
T+I
.I
-
C:jX:-l,j
+
Bi,,X:,
-
A:,+,Xlj+, - C:jXZj-1= 4.1 .
(3. I )
Equations (3.1) constitute a set of coupled nonlinear algebraic equations. The solution technique discussed is for flow through a rectangular duct, with the duct cross-section specified along the axis. In this case, the pressure drop, DP, that appears in Di,j for the velocity equations, is not known in advance and is to be calculated as a part of the computations. Equations are solved iteratively by using two iteration loops, one placed inside the other. The outer loop iterates on the unknown DP, and the inner loop solves the nonlinear algebraic equations, given a value for DP, iteratively. Assign a value to DP. At the beginning of the outer iteration loop, the value assigned is based on the upstream value of DP. For the subsequent iterations. DP is assigned a value as discussed in Step 7. Step 2. Calculate coefficients A, B, C and D for the velocity equations (2.5) evaluating V’s by using the information at x for the first iteration, and the information at x + Ax obtained in the previous iteration for the subsequent iterations. Step 3. Once A. El, C and D are known (2.5) constitutes a pentadiagonal set of linear algebraic equations. These equations are solved by using alternate direction implicit method. Step 4. Repeat Steps 2 and 3 for the enthalpy equations (2.9). Step -5. With the newly calculated values at x + Ax of h z, and P (obtained by adding DP to the value of P at x), calculate p; using the equation of state. Next, calculate Y: and ZT using u:, and p; with the help of formulae given in Appendix A. Step 6. Check newly calculated u:, and hl, for convergence. If the convergence check is met, proceed to the next step; otherwise repeat the steps starting from Step 2. Step 7. Check Y: and Z;, calculated duct half width and half height for the value of DP assigned in Step 1, against the duct dimensions specified at x + Ax. If the check is within the desired accuracy, the computation of flow from x to x + Ax is complete. If not. DP is adjusted, on the basis of the difference between the specified and the calculated duct dimensions. and Step
1.
82
M.S. Greywall,
Stream wise computation
of parabolic flows
the calculations started again at Step I, with the new value of DP. The formulae used to adjust DP for the succeeding iteration is the one given in [I]. This formula, based on the difference between the calculated and specified duct cross-sections, is applicable to both two and three dimensional flow computations. Spacing between the streamlines in the cross-stream plane changes result, at the beginning of each new integration step, streamtubes accordance with the definitions given in Section 2.1.
from x to x + Ax. As a are defined anew in
4. Sample flow computations The approach, presented in this paper, to compute three-dimensional flows by calculating velocity and enthalpy along the streamlines, was checked by making three sample flow computations. All the computations are for laminar flow through rectangular ducts. Computation of laminar flows is free from the uncertainties associated with the modeling of turbulent shear stresses in three-dimensional flows and, thus, provides a clearer assessment of the numerical method itself. All the computations are for the lower-left quarter of the duct. Each sample flow computation was carried out using one hundred streamlines; this corresponds to, including the grid points along the walls, a 11 x 11 grid in the y-z-plane. Streamtubes around these streamlines were redefined at the beginning of each new integration step, as discussed in Section 3. In the present computations. however, the same set of streamlines was used for the entire length of the duct. The integration step along the duct, Ax, was selected for each sample flow computation so that the locations of the computer output will match the locations of the reported experimental data; the integration step varied from one-tenth to one-twentyfifth of the equivalent duct diameter. The effects of the integration step size and the number of streamlines used are discussed in Appendix B. Velocity distribution at the inlet in all the sample flows was taken to be uniform across the duct cross-section. The following nomenclature is used in the figures in this section. -.-._
u Ul UC W, H D X
YandZ Re Pr Gr Nu,
Dots along the velocity distribution curves represent the streamline which the flow was computed. axial velocity. U at the inlet, centerline velocity, duct half width and half height, respectively, equivalent diameter (= 4 WH/( W + H)), distance along the duct measured from the inlet, distances normal to x measured from the lower-left duct corner, Reynolds number (= U,,D/v, where v is the kinematic viscosity), Prandtl number (= 0.72) Graetz number (= Pr Re/x). Logarithmic mean Nusselt number - Tw,,, =$Gr log T,,,,,, T hulk
T hulk
bulk temperature.
-
Twa,,
>>
’
locations
alon
M.S. Greywall,
Streamwise
computation
83
of parabolic flows
In Figs. 4-7 the results for flow through a square-sectioned duct are shown. Calculations are compared with the experimental data of Goldstein and Kreid [13] (Figs. 4-6) and Beavers et al. [14] (Fig. 7). Agreement between the theoretical and experimental results is good. In Figs. 8-l 1 the results for flow through a rectangular duct of 2 : 1 aspect ratio are shown. Calculations are compared with the experimental data of Sparrow et al. [1.5]. Once again agreement between the theoretical and experimental results is satisfactory. In the two flow computations just described only the finite difference equations for velocity were solved. Density was considered constant and, thus, the energy equation was bypassed. To check the finite difference equations for enthalpy a sample flow with heat transfer to the walls through a rectangular duct of 2 : 1 aspect ratio was computed. The inlet temperature and velocity distributions were taken to be uniform across the duct cross-section. The walls were maintained at a constant temperature of 100 degrees Celsius below the inlet flow temperature.
Fully developed
0.02 >iO s E” ._
6
0.4
0.3
E &
o.2
a 5
0.1
work date Goldstein et al. (13)
o Expt.
1 0
0
2
1 0
2
1
2
u/u,
Fig. 3. Development
0.5
r
of the velocity profile in a square-sectioned
ductdiagonal
Fully developed
&=0.0075
-
plane.
3
Present work 0 Expt. data Goldstein et al. (13)
u/u, Fig. 5. Development
of the velocity profile in a square-sectioned
duct-central
plane.
84
MS.
Greywall. 2.2
Streamwise
computation
of parabolic flows
r
I
-
Present 0 Expt.
work
data
Goldstein
1.0 1
0
I 2
I
I
I
4
I
I
6
8
10
12
&
Fig. 6. Center
line velocity
et al. (13)
J
14
x102
development
in a square-sectioned
duct.
6
5
4 N >” z3
Present
:
Re 2 651 v 735 a 1222 01564 ,5 1858
2
Expt.
work
data
Beavers
et al
1
C
I
1
I 2
I
I
I
3
4
5
*
Fig. 7. Pressure
I 6
I
7
x lo2
drop in a square-sectioned
duct.
I
8
(14)
MS.
Greywall.
Stream wise computation
of parabolic flows
85
The velocity and enthalpy finite difference equations were solved together as discussed in Section 3. Calculations were carried out using constant specific heat, perfect gas behavior, and the Prandtl number equal to 0.72. Results of the computations are shown in Fig. 12, and are compared with the calculations of Wibulwas as reported in [ 16, Table 52, p. 2201. To illustrate the velocity field in the cross-stream plane, in Fig. 13 are plotted the streamline locations in the lower-left quarter of the duct for air flow through a $ meter by : meter rectangular duct (second sample flow computation) at two locations along the duct. The beginnings and the ends of the arrows represent, respectively, the streamline locations at 1 meter and 5 meters from the duct inlet. Bulk velocity through the duct is 0. IS cm/set.
0.0128
r
0.03881
Present work Expt. date Sparrow
Fig. 8. Development
of velocity profile in rectangular duct of 2 : 1 aspect ratio-across
et al. (15)
long side.
l.O-
NII ? t
0.6
-
0
k
_
s
Present work @ Expt. data Sparrow
0
0
0
1 0
1
2 1
2
et al. (15)
2
u/u,
Fig. 9. Development
of velocity profile in rectangular duct of 2 : 1 aspect ratb-across
short
side
86
M.S.
Grevwall,
Stream
of parabolic
wise cornputatiorl
jlows
Present @ Vertical 0 Horizontal
1.0.
I 3
I 2
I 1
0
traverse traverse
I 4 X D Re
x 10
Sparrow
et al. (15)
I
I 6
I 5
work
f=xpt, date
I 8
7
2
Fig. 10. Center line velocity development in rectangular duct of 3 : I aspect ratio.
8.0-
5.0
-
3.0
2.0
-
-
0
1
2
3 D
X
4 x I0
6
Present B
Re 897
Expt.
t
:ig
Sparrow
o A 0 m
1233 1416 1876 1986
6
work
data n+ .-1. . . ,_ _.
7
2
Fig. 11. Pressure drop in rectangular duct of 2 : 1 aspect ratio.
(15)
M.S. Greywall. Stream wise computation of parabolic flows
87
E 2’
5
2t 20
Fig. 12. Variation of logarithmic with constant wall temperature.
0.4-
Fig. 13. Motion
I
I
50
mean
Nusselt
number
**-e-e_,--+
of the streamlines
in the cross-section
-
- Wibulwas
I
100
I
150
with Graetz
--+
work
Present
--
200
number
(16)
J
220
for flow through
a rectangular
duct
--.+
plane-figure
is for the lower-left
quarter
of the duct
5. Concluding remarks A new approach to compute three-dimensional parabolic flows has been presented. The dependent variables in the present approach are the streamline velocity and the streamline coordinates in the cross-stream plane. Strictly speaking, the present approach is not Lagrangian, as one does not follow the individual fluid particles, but rather follows the flow through individual streamtubes. This approach, however, seems closer to the Lagrangian than the Eulerian flow description. The advantage of the present approach is that one needs to solve only one set of finite difference equations---for the streamwise velocity u-rather than three sets of equations for u,, u,, and u,. The disadvantage is that extra computations are required to calculate streamline coordinates.
h4.S. Greywall.
88
Stream wise computation
of parabolic flows
Appendix A In this appendix formulae which supplement the finite difference equations given in Section 2 are presented. The lower-left corner of the duct is taken as the origin of the y-z-plane. A. 1. Mass flow rates Let
Ri.j
=
&jUi.j.
Formulae
for stream
(2,2):
= :(2R,,, + Rs.2+ R,,)
Pi.;
. ;( Yi - YJ +;(Z, - Z,) .
‘z&f = y; . 1:(3R2.2+ Rz.x)( YJ . ;(Z.i - Z2) , ‘Pyt, = yf . +(3R,,, + R?.z) . +(Y3 - Y?)(Z,) . Formulae
for stream
(i, 2), i >2:
Formulae
for stream
(2, j),j > 2:
l~l,T = :(2Rz,j + R3.j + R,.j+l) . i( Y?-
Yz) . $(Z,+I - 4) 3
Ti,t = ~7 *i(3Rz.j + Rz,j+l)( Yz) . i(Zj+l - Zj) 3 lpp4 = ry *$(3Rz., + Rz,j-I)( YL) *$(Z, - Zjm1)3 Ti.7 = i(2Rz.j + Rz,j-1 + R,,j) . $( Y.3- J’z) . $(q - &I> Formulae
for stream IJZ-C+T = WLY+
A.2.
(i j). i and j > 2:
i(2Ri.j + Ri+l,j + Ri,j+l) * $(J’+I - X) . t(4+1-
~$(2Ri,,
+R,j+,+
Ri_I,j).$(X-
$(X - K-1) * %Zj- -T-l)
$(;?Ri,j + Ri_I,, + Ri,j-1) .
IJY~,T=
i(2Ri.j + Ri,j_1 + Ri+l,j) . t( X+I -
Surface and cross-sectional areas
$(Y3 -
Yz)
+
Y,)Ax
,
4)
3
X-I>~~(G+I-Z),
IJFL,T =
A;’ =
.
9
Yi> * k
M.S. Greywall.
Streamwise
computation
89
of parabolic flows
AsJ = $(y,,, - k;-,)Ax, A;'
= &Z, - Z, + Z&Lx,
A;'
= +(Zj+, - Z,-,)Ax,
A+
r.1
=-@‘A: hX2.
A.3. Streamline coordinates Y t, ZT After the velocities and densities have been calculated used to calculate the streamline coordinates at x + Ax. Let RT, = p:iU; and
For a given value of j. by summing i = I, one obtains for j = 2,
at x + Ax, the following
the mass flow rates through
the streamtubes
procedure
is
with i = 2 to
(A-1) and for j > 2, (Zf - ZT_,)Yf
By adding
(A.l)
16P’,; = 77(3Rl,i+ R;i_l)+
16P;,;_, rr-1(3Rl,ip1 + R&)+
I l6 z? Gi’i*
64.2)
to (A.2) for j = 3 to j = J, one gets 64.3)
z:y:=c, I
where c, represents the sum of the right-hand sides. By using the duct aspect ratio, specified at x + Ax, (A.3) can be solved for Y:. With Y: known, (A.l) and (A.2) can now be solved for Y G, Yf, . . . , Y: are calculated similarly. z;, z;. . . . ) Z:. The coordinates A.4. Mass, momentum and energy correction factors For laminar flow, the streamlines adjacent to the duct walls can be chosen near enough to the walls so that the velocity and enthalpy distributions between the walls and the adjacent streamlines can be approximated by linear functions. In this case V& and V:,i, introduced in Section 2.2, are given by the general expressions for V’s given in Section 2.1, and the
M.S. Greywall.
90
correction
factors
y,
Streamwise
computation
and 6 take on the following
a,
of parabolic
flows
values
These are the values for y, cy, and 6 used in the sample flow computations given in Section 4. The finite difference equations for the enthalpy were solved by using the wall temperature as the reference temperature; for this case the p’s become equal to the corresponding cy’s. For turbulent flow the correction factors y, (Y, 6 and p, and V$ and V:j will depend on the turbulence mode1 chosen for the near wall region in three-dimensional flows. Their calculation can be carried out by using three-dimensional extension of the procedure given in [l] for the case of two-dimensional flows. A._5. Viscous forces OHstream (i, j) The viscous forces exerted on the stream (i, j) by the neighboring streams are calculated as follows. The force exerted on the stream (i, j) by the stream (i, j + 1). call it F(top), is given by = Ap
F(top)
(a24
dz) .
(A.4)
where A is the interface area between the streams (i, j) and (i, j + l), and p is the effective viscosity at this interface. By substituting the appropriate symbols for A and ,u defined in Section 2.1 (see Fig. 3) and by approximating the velocity gradient by l&j+1
-
Ui.1
Zj+, -Zj one obtains
’
from (A.4)
F(top)
= AsJp{j+l/z ‘g,+’ I >I .
By using the definition F(top)
=
V:]+I(W,j+l
Similarly one obtains call it F(bottom), F(bottom)
of V:,,,,
(AS) kj)
is rewritten
= V:j(ni,j -
as
.
for the viscous force exerted
The net force on the stream by,
-
(A.5)
I
/+I
(A.6)
by the stream
(i. j - I),
(A.7)
u,,,_,).
(i, j) by the neighboring
(i, j) on the stream
streams
(i, j + 1) and (i. j - 1) is thus given
MS.
Greywall,
Streamwise
computation
of parabolic flows
91
V{j+l(Ui,j+lUi.j)Vi.j(Ui,j - Ui,j-l)e
F(top)-F(bottom)=
64.8)
By adding to (A.8) a similar expression for the net force exerted by the streams (i + 1. j) and (i - 1, j), one obtains the required relation for the viscous forces exerted on the stream (i, j) by the neighboring streams.
Appendix
B
In this appendix the effects of the integration step size along the duct axis, and the number of streamlines employed, are discussed. In Fig. 14 is shown the pressure drop in a squaredsection duct calculated with various integration step sizes and with various numbers of streamlines. The calculations were carried out for the lower-left corner of the duct, and the numbers given in Fig. 14 for the streamlines employed are for the quarter section. The streamlines were arranged in equal numbers of rows and columns (see, for example, Fig. 13). The spacing between the streamlines was non-uniform; the general nature of the spacing used can be seen from the velocity distribution figures given in Section 4. The integration step size, Ax, shown in Fig. 14 is normalized by D, the equivalent duct diameter. All the calculations were started with an initial integration step equal to 1/512th the final size. The initial integration step was then doubled at the end of each integration step until it reached its final value. The pressure drops for the three cases: Ax = O.lD and 100 streamlines; Ax = O.lD and 225 streamlines; Ax = 0.05D and 100 streamlines, were found to be very close (A&& = 6.478,
7_
6-
6-
Number of streamlines
AX/D 0.1 0.1 0.1 0.06
--0
1
2
3
4
6
-
Fig. 14. Effects of the integration sectioned duct.
X D Re
6
x10
0.2
100 7
8
6
2
step size and the number of streamlines
on the pressure drop in a square-
M.S. Greywall.
92
Streamwise
computation
of parabolic
flows
6.444, and 6.483 at x/DRe equal to 0.09 for the three cases, respectively). The difference was too small to be shown in the figure; thus, the three cases are shown as a single curve. The calculations were carried out using FORTRAN IV programming language on IBM 370/3031 computer operating under CMS Batch. The CPU times to carry out the computations, for duct length equal to 10 times the equivalent duct diameter, were as follows: Ax = O.lD;
25 streamlines, 100 streamlines, 225 streamlines,
100 streamlines;
Ax = 0.050, Ax = O.lD, Ax = 0.20.
CPU time = 1.28 mns, CPU time = 4.00 mns, CPU time = 10.60 mns, CPU time = 4.31 mns, CPU time = 4.00 mns, CPU time = 6.00 mns.
For a fixed value of Ax, the CPU time is scaled, as expected, by the number of streamlines used. For a fixed number of streamlines, however, the CPU time does not necessarily decrease when Ax is increased. As Ax was increased, although, the number of integration steps required to reach the end of the duct decreased, the number of iterations needed for each integration step increased.
Acknowledgment We would like to thank Dr. John M. Smith of NASA-Lewis Research Center for his valuable discussions and suggestions during the course of this work. The work was supported by NASA/DOE Grant NSG-3186.
References [l] M.S. Greywall. Streamwise computation of duct flows. Comput. Meths. Appl. Mech. Engrg. 21 (1980) 231-247. [2] L.S. Caretto, R.M. Curr and D.B. Spalding. Two numerical methods for three-dimensional boundary layers, Comput. Meths. Appl. Mech. Engrg. I (1972) 39-57. [3] S.V. Patankar and D.B. Spalding, A calculation procedure for heat, mass, and momentum transfer in three-dimensional parabolic flows, Internat. J. Heat and Mass Transfer 15 (1972) 1787-1806. [4] A.J. Chorin. A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1967) 12-26. [S] J.A. Miller, Laminar incompressible flow in the entrance region of ducts of arbitrary cross section, J. Engrg. Power Trans. ASME 93 (1971) 11>118. [6] A.F. Emery, P.K. Neighbors and F.B. Gessner, The numerical prediction of developing turbulent flow and heat transfer in a square duct, J. Heat Transfer Trans. ASME 102 (1980) 51-57. [7] D.W. Roberts and C.K. Forester. Parabolic procedure for flows in ducts with arbitrary cross sections, AIAA J. 17 (1979) 33-40. [8] V.S. Pratap and D.B. Spalding, Fluid flow and heat transfer in three-dimensional duct flows, Internat. J. Heat and Mass Transfer 19 (1976) 118>1188. [9] J.F. Nash, The calculation of three-dimensional turbulent boundary layers in incompressible flow, J. Fluid Mech. 37 (1969) 625-642. [ 101 K.C. Wang, Three-dimensional boundary layer near the plane of symmetry of a spheroid at incidence, J. Fluid Mech. 43 (1970) 187-209.
M.S. Greywall,
Streamwise
computation
of parabolic flows
93
[ 1 I] F.H. Harlow and J.E. Welsh, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Phys. Fluids 8 (1965) 2182-2189. [12] B.E. Launder and D.B. Spalding. The numerical computation of turbulent flows, Comput. Meths. Appl. Mech. Engrg. 3 (1974) 269-289. [I31 R.J. Goldstein and D.K. Kreid. Measurement of laminar flow development in a square duct using a Laser-Doppler flowmeter. J. Appl. Mech. Trans. ASME 34 (1967) 81.3818. [I41 G.S. Beavers, E.M. Sparrow. and R.A. Magnuson, Experiments on hydrodynamically developing flow in rectangular ducts of arbitrary aspect ratio. Internat. J. Heat and Mass Transfer 13 (1970) 689-702. [ 151 E.M. Sparrow, C.W. Hixon, and G. Shavit. Experiments on laminar flow development in rectangular ducts, J. Basis Engrg. Trans. ASME 89 (1967) 116-124. (161 R.K. Shah and A.L. London, Laminar Flow Forced Convection in Ducts (Academic Press. New York, 1978).
METHODS
NORTH-HOLLAND
IN APPLIED
PUBLISHING
STREAMWISE
MECHANICS
COMPUTATION
of Mechanical
Wichita State Unitlersity.
Received
A
new
the
approach
computation
dependent difference mass.
of
velocity
three-dimensional
plane,
of the chosen
the chosen Results
obtained streamlines;
of the calculations.
data for flow through
rectangular
19 November
1981
received
10 March
three-dimensional
in the present
equations
manuscript
parabolic
along a set of streamlines.
variables
cross-stream around
to calculate
by calculating
36 (19X3) 7 I-93
flows
approach
are
the
The
Euler’s
the streamline
Wichita. KS h72Ob’. U.S.A.
1982
flows
is presented.
The dependent
three
velocity
are the streamwise
streamlines.
by applying
PARABOLIC
S. GREYWALL
Engineering,
Revised
computed
ENGINEERING
OF THREE-DIMENSIONAL FLOWS
Mahesh Department
AND
COMPANY
velocity velocity
theorem
are
calculated
approach.
ducts;
is satisfactory.
the agreement
flow field is used in
In contrast,
and the two coordinates.
streamwise
coordinates
The
commonly
components.
momentum
based on the present
variables
is calculated
from
to streamtubes from
are compared
the in the
the finite
constructed
the conservation
of
with the experimental
1. Introduction Viscous flows have been calculated for a long time by separating the flow into a boundary layer and a potential flow region. The boundary layer then is calculated by using integral methods and the potential flow by using inviscid flow equations. This approach does not require large computers, and even can be carried out by use of hand calculations. This approach, however, has its limitations. For example, in calculating flow through nozzles, channels, and diffusers this approach is limited to the region near the entrance. As the boundary layer thickens, the flow cannot be meaningfully separated into a potential flow region and a boundary layer region. Further, this approach cannot describe the corner effects, such as the secondary flows, or adequately describe flows in which the body forces are nonuniform in the cross-stream plane (such as magnetohydrodynamic channel flows). In the past two decades, with the advent of high speed computers, many numerical methods to compute velocity and temperature fields without separating the flow into a boundary layer and a potential flow region have been developed. In this paper, a new approach to calculate three-dimensional compressible viscous flows is presented. The approach is an extension of the author’s previous work [l] to three-dimensional flows. The method, as presented in this paper, is for ‘parabolic’ flows through rectangular ducts, and its possible extensions to other ‘Parabolic’ or ‘boundary layer’ flows types of three-dimensional flows is under investigation. are characterized by the existence of a predominant how direction along which downstream conditions have a negligible influence on the upstream conditions. This assumption allows use of a marching integration procedure. A lucid discussion of the parabolic flow assumption is given by Caretto et al. [2], and Patankar and Spalding [3]. 00457825/83/0000-0000/$03.00
@ 1983 North-Holland
72
M.S. Greywall,
Stream wise computation
of parabolic flows
Methods developed so far to compute three-dimensional viscous flows are, among others, those given in (2-101. These methods differ in their finite difference approximations, and whether the flow is assumed to be parabolic or not. From among these methods, those proposed for calculation of confined flows differ further in the procedure used to calculate the pressure field. These methods are discussed here briefly. Included in the discussion is the method for the time and convective terms, and Dufort-Frankel for the ditiusion terms. The transient viscous flows, it has some features in common with the methods developed later for of Chorin [4] calculates computation of three-dimensional parabolic flows. The method incompressible elliptic flows. Finite difference equations are obtained by using the leapfrog method for the time and convective terms, and Dufort-Frankel for the diffusion terms. The pressure field is calculated by coupling the incompressible continuity equation to the momentum equations via an artificial variable density and an artificial equation of state. Miller [S], on the basis of the work of Chorin, computed steady incompressible three-dimensional elliptic flows through rectangular ducts. In the method of Harlow and Welsh [ 111 the finite difference equations are obtained by using forward in time and center in space differencing. The pressure field is calculated from an elliptic equation obtained from the continuity equation and the divergence of the momentum equations. The method uses a ‘staggered’ grid in which the velocities are defined at the cell boundaries and pressure at the cell centers. A similar type of staggered grid is used in the cross-stream plane, later, in [2,X 6,8] for the computation of three-dimensional flows. In these methods the cross-stream velocities are specified at the cell boundaries and the axial velocity, along with all the other variables. at the cell centers. Methods [2,3,6, 71 calculate three-dimensional confined parabolic flows. In these methods the pressure field in the cross-stream plane is calculated from an elliptic pressure equation obtained from the continuity equation and the cross-stream momentum equations. The axial. along the duct axis, pressure gradient is calculated to insure conservation of mass along the duct axis. In these methods Caretto et al. [2], and Patankar and Spalding 1.31 obtain finite difference equations by using the control volume approach; Emery et al. [6] use DufortFrankel for the axial velocity momentum equation and implicit differencing for the crossstream velocity momentum equations; Roberts and Forrester [7] use one-sided upstream differencing in the flow direction and second order differencing in the other two directions. Pratap and Spalding [8] extend the method of Patankar and Spalding to ‘partially parabolic’ flows. Methods of Nash [‘i’] and Wang [lo] are for external (pressure field specified) threedimensional parabolic flows: Nash uses forward explicit differencing, and Wang uses a Crank-Nicolson type of differencing scheme. All the methods described in the preceding paragraph use the Eulerian system. The flow field is computed by calculating the velocity components u,, u,, and u, along a set of spatial grid points fixed in advance. In contrast, in the present approach, the flow field is computed by calculating the streamwise (not axial) velocity, u, along a set of chosen streamlines. Streamlines are identified by a pair of indices, (i. j), such that the coordinates of the streamline (i, j) in the cross-stream plane are Yi and Z,. Coordinates of the streamlines in the cross-stream plane are calculated as part of the flow computations. The dependent variables in the present approach are Ui,iTthe streamwise velocities, and (Yi, Zj). the coordinates of the streamlines. To calculate the streamwise velocities streamtubes are constructed around the individual streamlines such that the sum total of all the streamtubes fills the duct. One such streamtube. constructed around the streamline (i. j). is shown in Fig. 1. Finite difference equations for ld,~,
M.S. Greywall.
Streamwise
computation
of parabolic flows
An adjacent
---
X
Fig. 1. Streamtube
73
streamline
Streamline
(i,j)
Duct
a xis
X+&X
(i, j) constructed
around
the streamline
(i j).
and hi.,, velocity and enthalpy along the streamline (j, j), are now obtained by applying Euler’s momentum theorem and the first law of thermodynamics directly to the flow in the streamtube between x and x + Ax. After the streamwise velocities and densities have been calculated at x + Ax. the streamline coordinates at x + Ax are cafculated next by using mass conservation. From (Yi, 2’) and (Y’, ZT), coordinates of the streamline (i j) at x and x + Ax, one can calculate the slope of the streamline (i, j), and then using this slope decompose u,., into its components u;i, uz’, and ufi.
2. Finite difference equations 2 1. Background In most numerical approaches to the computation of the viscous flows one starts with the partial differential equations-Navier-Stokes equations, the energy equation, and the continuity equation-obtained from the conservation principles. These equations, either in their primitive forms or in transformed forms, are then differenced to obtain the needed finite difference equations. In contrast, in the present approach the flow in the duct is partitioned into a finite number of streamtubes and the finite difference equations are then obtained by applying the conservation principles directly to the flow through the individual streamtubes. Consider the flow through a rectangular duct with its axis along the x-direction. To define the various streamtubes lines Y2, Y3, . . , , Yiq . . . , YI are drawn parallel to the z-axis (left duct wall), and lines Z,, Z,, . . . , Z”, . . . , Z, parallel to the y-axis (lower duct wall), as shown in Fig. 2. Lines Y1 and 2, (not shown in Fig, 2) are drawn through the middle of the lower wall and the middle of the left wall, respectively. The lower and left walls are also referred to as Z, and Y, lines. Line Yi+l/z is drawn in the middle of Yi and Y,,,, lines Yi_1,2. Zj+j/zr and Zj-l,z are drawn similarly. Various streamtubes, called streams for short, are now defined as follows. Stream (2,2): flow bounded by the left wall, the lower wall, Z2+1,7_,and Y2+1,7. Stream (i. 2) i = 3, . . . , I: flow bounded by the lower wall, Z2+1,2r YI:-l,z, and Yi+i/z. Stream (2, j), j = 3, . . . , J: flow bounded by the left wall, Y2+l,z$ Zj-x,z, and Zj+r/z. Stream (i, j), i = 3. . . . , I, j = 3, . . . , J: flow bounded by Yi_-1,2. Y<+j/z, ZjP,,z, and .Zi+r,2. Streamline (i, j) is defined as the streamline passing through the intersection of Yi and Zi.
M.S. Greywall,
74 z
r2 ly 2+%
’ I
yi_‘/
’ I
’ I
_L_L__L____L___Zj, I I I
I
I /q/j, ‘,// : *;‘j:
I , I
i,j ;;g+ ;_ ‘,/ ,,,; ij /, ,,///
L _: I
yi+‘/t
(i,j) l&$~;’
:
I --
Fig. 2. Definition
--zj 2.
J-S
I _ _ L _=_2+: z2 ,:
i.1
1.1
LI I
Tgg;;“: $ (1.2)::
‘j+l/,
I
., ‘1
I
L ,
1
I
I Y”“’ w;;, q,: , I
_L_L__L____L___~3 I I I
ggg ,(2#2)
of parabolic flows
yi+l
yi T3
I
Stream wise computation
of the various
I
Y
streams.
+.I -
%
Fig. 3. Control
Yi+K
volume for the stream
(i, j).
Other nomenclature is as follows: !Py,! is the mass flow rate through stream (i, j). The mass flow rate !Pc is partitioned further into four parts !PlT, q;,+, V’i,, and q’:; ; mass flow rates through the upper-right, upper-left, lower-left, and lower-right quadrants, respectively (see Fig. 2).
AC is the cross-sectional area of the stream (i, j). To apply conservation principles control volumes are constructed for each stream distance Ax along x. The control volume for the stream (i j) is shown in Fig. 3.
using a
As1, A?: surface areas of the control volume normal to the y- and z-axis respectively. ui,,, Ti,,, hi,j, pi,,, ciJ: velocity, temperature, enthalpy, density, and constant pressure specific heat along the streamline (i, j) at x. U;, etc.: variables along the streamline (i, j) at x + Ax. Ui,j: average of Ui,j and u ;. ZA&, ufj: of iii,j the X, y, and z directions. Yi,Zj: coordinates of the streamline (i, j) in the cross-section plane at x. Y:, Zi+: coordinates of the streamline (i, j) at x + Ax. Si,j s
12
TU r,j.
+_1
t2 = 2lA i,j . ',I I p i+ 1/2,j, p f- 1/2,j, p
s
{j+ ,127 p :j-,121 effective viscosities, including the turbulent contributions, for the momentum transfer across the surfaces of the control volume, as shown in Fig. 3, evaluated using the information at x. pfzIIz,j, etc.: effective viscosities evaluated using the information at x + Ax. K!+112,j, etc.: effective thermal conductivities, defined similarly as the effective viscosities.
M.S. Greywall, Stream wise computation of parabolic flows
V’s and H’s with superscript
Corresponding
+ are defined
in terms of p+, K’,
Y’, and Z’.
In estimating momentum and energy fluxes through a streamtube u and h are assumed, in general, to vary linearly between the neighboring streamlines. Between the walls. however, and the adjacent streamlines (streamlines with i or j equal to 2 in Fig. 2). a nonlinear variation of u and 11 is allowed. Thus, the mass. momentum, and energy flux through the streams adjacent to the walls are expressed using coefficients y, (Y, p, and 6 defined as follows (see Appendix A for further discussion):
The following relations are for the lower quadrants of the streams along the lower duct wall. Formulae for the lower-left quadrants are obtained from these relations by substituting !P;; for ?&
and for the lower-right
quadrants
by substituting
U = $(3Z(j.2+ Ili+l,z) 1
AY = +(r;+, - Y), 6 = $3h,,,
+ hi+t.,)
PT; for ?P,‘,.?,
and
&? = &$
+ iuf+,.J, Z2
Z2 Pi.2 = A Y
AY
I
pudz
=
y:@AYZ,
,
AY -3
ZZ
0
pUUdZ~CXfTPi,~G, I0
pu tu2dz = S; ‘&p,
Ay
f0
puhdz
Coefficients 77, a;, ST, and py, that take into account the left wall, are defined similarly.
= /?‘;‘iQh.
the non-linear
variations
normal
to
76
MS.
Greywall,
Strearnwise
computation
of parabolic flows
In evaluating the various integrals for the near-wall correction factors, p has been set equal to some constant p. For incompressible flows p is, of course, just equal to p. In computing two-dimensional compressible flows, p set equal to the average of p between the wal and the adjacent streamline gives satisfactory answers. For the three-dimensional compressible case the appropriate value to be used for /s is still under consideration. 2.2. Finite difference equations for Idi,, Applying
Euler’s
momentum
Momentum = Pressure
theorem
to stream
flux outlat r+A.r- Momentum forces + Viscous
(i, j) between
x and x + Ax one obtains,
flux ini at x =
forces + fi,j(Body
forces) .
denotes the quantities in the preceding bracket [ ] with U+ by replaced by 24. This symbol is used similarly later on. In calculating pressure and viscous forces, it is assumed that the streamline inclinations with with u.*. Thus, in respect to the duct axis are small; that is, U” and uz are small compared calculating pressure force A;j+ the cross-sectional area of the streamtube (i, j) projected onto the Y-Z-plane, is used. In calculating viscous forces, the velocity gradient is based on An, (see Fig. l), the distance between the neighboring streamlines projected onto the Y-Z-plane. Pressure force on stream (i, j) is estimated as 141
- DP A;,
(2.3)
where DP is the pressure difference between x and x + Ax and is assumed to be the same for all the streams. The assumption that DP is the same for all the streams restricts the analysis to straight ducts without any abrupt change in area along the duct. The viscous force on a surface between x and x + Ax is approximated by the average of the viscous forces evaluated by using the effective viscosities and the velocity gradients at x and x + Ax. One thus obtains for the viscous forces on the stream (i, j) (see Appendix A), Viscous
Forces = i[ V’+ t+l,j(l4T+l., - 44:,)+ -
V:j+(U:j
-
14:
1.~)~
V::+I(UZ,+ Vi,+(Z4lj
-
I - Z4lj) L4,tl- I)]
+ +[
V’.
U+ +
V, U] .
(2.4)
M.S. Greywall,
Streamwise
computation
of parabolic flows
77
For streams along the walls, momentum fluxes are expressed using the non-linear correction factors introduced in Section 2.1. Viscous forces exerted by the lower and the left wall on the adjacent streams are, respectively, set equal to V&U~,~and V;,jU,,z. The calculation of Vi2 and Vi,, is discussed in Appendix A. After the various estimates in (2.1) have been substituted, all the terms involving Ucj are collected on the left-hand side, and the resulting equation is organized in the form
Coefficients A, B, C and D for the various streams are given next. 1 Z, are i times the ly’s defined in Section Note. In equations that follow all the q’s, except !J’_, 2.1, and the V’s are $ times the V’s defined earlier. Coefficients for stream (2.2):
Coefficients for stream (i, 2) i > 2:
Coefficients for stream (2. j), i > 2:
78
M.S. Greywall.
Coefficients
for stream
-7.3. Finite difference
Applying
computation
of
flows
(i. j), i and j > 2.
equations
for /Ii,,
first law of thermodynamics
(i, j) between
to the stream
x and x + Ax, one obtains,
Energy flux outlat x+A.Y - Energy flux inlat X= = Heat added by conduction + Work done by the viscous forces + Q,,(Internal heat generation)+ wi,j(Work done by the body forces).
(2.6)
For the streams (i, j) with i and j > 2, various terms in (2.6) are estimated as follows. The energy flux terms are estimated in the same way as the momentum flux terms and can be obtained by substituting (h + s) for u in (2.2). The heat conducted across the top surface of the control volume (see Fig. 3) is estimated as
Work done by the viscous forces
$[V<~+,(U~j+l
-
14,:I)(i(Ut+l
=f[V{,‘+,(.S~,+[ -
aS:j)]
on the top surface +
Ulj))l
+t[V+S++
+
is taken to be
$[V+Ll++
VU] =
VS] .
(2.X)
Work done by the viscous forces and the heat conducted across the other three surfaces of the control volume are estimated similarly. After the various estimates in (2.6) have been substituted, all the terms involving h:, are collected on the left-hand side and the resulting equation is organized in the form I
+
-A,,jh;+i,j-
I
+
+
C,.jh-1.1 + Bi.,h ,,,+I-A{jh,T,+l-
CfjjlZ,p, =Dt.,.
Coefficients A. B, C and D for the various streams are given next. Note. In equations that follow all the P’s, except W.5, are $ times Section 2.1, and the H’s are $ times the H’s defined earlier. Coefficients
for stream
(2.2):
(2.9)
the
q’s
defined
in
MS.
Greywall,
Streamwise
of parabolic flows
computation
79
02.2 = Q2.2 + W,,, + h,,z( Yf;,t + #of !P;.; + El:,,) + h,,(2!P;.;
+ 3pS !P,; + Pf PT., + 3p5 P’t; - H:,, - H;., - Hi., - Hi,,)
+ hz..J ‘P’tt + PP Y’,: + Hi.3) + s&Y’;.: + s2,42ly;.;
+ 36; !P,; + s; ?K,
+ 6; W;,s + V:.,)
+ 36; Ye;
- vi.2 - v;,, - v;,, - VS.,)
+ s2,3(!I’;.; + 6; T’;;
+ Vi,,) + h ,.,H:,, + h,l H;,z
+ h,,,H:.:; + h,,H::
+ s:,,(-
+ 3c5;*;,:
+ s;,2(2K;
!P;.; - S; P;,, + V:r’,)
+ s; YP_2 + 3s; w;
+ Sl.3(- !P;.; - sp !P;,; + v:.;,
+ v:;
+ vg
+ v;r+, + VC)
.
Coefficients for stream (i, 2), i > 3: A:2 = - !l’;; - Pf W;; + Hf;,.2 , C:,= -!&-/3f!P;;+ Bi,2=2~~~+3Pf~(i17_+Hf=,,2+H:_~+H:f+H:2+, A? = - IJY(?+ + H-!’ c:, = 0 . I.3 .
H:;,
D;;f = Qi,2 r;’ Wi,2
+Si+l,Z(!P~~+61*~~+
Vf+l,Z)+Si~,,Z(W;~+~flVi.2+
Vf.2)
+ Si,2(2T:z + 36: Py; - Vf+ 1,’- Vf3 - V{2 - V:2) + Si,3(q?:
+ V$) + h,lH:2
+ h,lH{i
+ s:+,,J(- !P’1+2 - Sf !+Jq-y+ VfIl,,) + s:_,.2(- !P’;; + s;22w?+I.2+ 36; p- 1.2+ Vf’,+ 1.2+ v:.; + v:; +s;3-
!P.:I;+
sf 9,
+ v::
+ v:;
v;:).
Coefficients for stream (2, j), j > 2: A’.=
c:, = 0 ,
-~V+‘?+ff’+.
B:; z 21$!+
3&+
+ Hi; + H;,;+l + H:; + H:,; , c;~~ = - ~7 Pi,; - !P;,; + H:; A;,j = - ‘I’;.; - /3; ‘Pi,; + Hi,;+ 1 , Di,j = Q1.j + W2.j
,
+ h,j ( ‘JT,: + Hi,) + hz.j(2Pz,Y + 3p7 *y.J! - H:,j - HSI,j+l- H:j - H:j) + hz,j+,( !$‘z.J+ ~7 Pi,f + H;,j+l) + hz.j-l(pyqY.Y + s3.j( qi.7 + Vi.j) + &,j(2WZJ + 363 YP2.ys2,j+*(Tl,f
+
h1,jHi.j + h1,jHi.T + Sl,j(- pZ,y + V:,$)
- &j(2Pl,7 +
S2,j+l(-
+
67
P2.f
Vs,j+,) +
+
+
S2,j-1(6;
+ 3S7!P’;q + V$,J + V;,T+l + ql,f
-
SyP2.J
+
V5,j
%,y
Vi,j
+
+ %.T + Hij) -
VT&+1
+
Ti,i
-
+
Kj
-
vCj)
vii)
vCj)
V:,T+l,>+ S2,j-1(-8~pS.j
- pii
+ v:,T) ’
80
M. S. Greywull,
Coefficients
for stream
Stream
wise cornputution
of parabolic
flows
(i, j), i and j > 3: c’
A:, = - !P:;‘+ H!:,,, ,
= - I)-!’ + H?t r., . 1.1 I.,
+ H:;+, + H:J + H::
B,,, = ‘Py; + Hf:,,,
A-’ = - I$% + H!?r.,+l 7 I., I.,
c:,
3
= - ‘P::; + H;I; ,
4.j = Qt., + W, + h,+,,,(v:;‘+ + h,,(‘J$-
IX+,.,)+
h,- &K;‘+
H:,)
HJ+,,,-H:,+,-H:,-H:j)
+ h,.,+,(T:IT + H:,+l)+
h+I(CIJ
+ HL)
+ s;+,.,(!Plt:)+ Vf+,.,)+ him.,+ +
S,j(
PI!:’
-
Vf+
1.1
-
K,)
V:j+ I - vt] - v:,)
+ S,.,+l(T:tT + Vij+l) + Si.,pI( %17 + v:]) + ST+I,,(- P,f:’ + VfI1.j) + ST-~ I.j(- qt,” + vt:) - s:,pP:y+
v;:,.,
+s:,+,(-P:;;+
+ V:: + V$)
v:;,,+s:.,.
of u and
2.4. Decomposition
+ v::+,
the secondary
,(-El,+
Vii).
flows
Computation of secondary flows. defined as the velocity components normal to the main flow, is discussed in this section. After y+ and Zt have been calculated, the streamwise velocity iii,j, the average of U,,j and Cl:j. then can be decomposed into II;,, Iit+ and U;j, by using the formulae,
*
11i.1
=
ii,,, where
AX Ar
u:,_AY ’
Ax is the integration AY=
Y’-y,,
iii.1 -
and
Ar
uf,_AZ -_C,,j
(2.10)
Ar
step, AZ=ZT-Zj
and
Ar’=AXz+AY’+AZ’.
The calculation of Y’ and ZT as discussed in Appendix A is based on the conservation of mass alone. Velocity components uy and u’, obtained from Y’ and ZT thus calculated, are the secondary flows associated with the flow development in the ducts-may the flow be laminar or turbulent. In turbulent flows through rectangular ducts, besides the primary shear stresses which are in the streamwise direction, there are also shear stresses in the plane normal to the streamwise direction (see, for example, [12]). Force due to these stresses will produce bending of the streamlines in the z-x and z-y-planes, and thus displace Y’ and Zt beyond that displaced by the conservation of mass. For example, force per unit area AT_ due to the cross-stream shear stress r,, distribution, will bend the streamline (a j) according to the formula (2.11)
M.S. Greywall,
Stream wise computation
of parabolic flouts
81
where R,, is the radius of curvature in the z-x-plane. With Rg, (approximated as a straight line for small Ax) known one can calculate displacement of ZT (and similarly of Y’) caused by the cross-stream turbulent shear stresses. Contributions to II.’ and u”, as calculated from (2. lo), by this additional displacement of the streamlines in the cross-stream plane represents, as viewed from the present computational approach. the secondary flows associated with turbulent flow through noncircular ducts. It should be noted that the turbulent shear stresses in the cross-stream plane do not enter into the finite difference equations for the streamwise velocity presented in Section 2.2. These equations are valid for laminar and turbulent flows as long as the streamwise shear stresses, appearing in the equations via V’s, are modelled appropriately.
3. Solution technique A methodology to solve the momentum and energy finite difference equations given in Section 2 is presented in this section. This is the methodology used for the sample flow computations presented in Section 4. Equations (2.5) and (2.9) are of the general form - A :,jX
T+I
.I
-
C:jX:-l,j
+
Bi,,X:,
-
A:,+,Xlj+, - C:jXZj-1= 4.1 .
(3. I )
Equations (3.1) constitute a set of coupled nonlinear algebraic equations. The solution technique discussed is for flow through a rectangular duct, with the duct cross-section specified along the axis. In this case, the pressure drop, DP, that appears in Di,j for the velocity equations, is not known in advance and is to be calculated as a part of the computations. Equations are solved iteratively by using two iteration loops, one placed inside the other. The outer loop iterates on the unknown DP, and the inner loop solves the nonlinear algebraic equations, given a value for DP, iteratively. Assign a value to DP. At the beginning of the outer iteration loop, the value assigned is based on the upstream value of DP. For the subsequent iterations. DP is assigned a value as discussed in Step 7. Step 2. Calculate coefficients A, B, C and D for the velocity equations (2.5) evaluating V’s by using the information at x for the first iteration, and the information at x + Ax obtained in the previous iteration for the subsequent iterations. Step 3. Once A. El, C and D are known (2.5) constitutes a pentadiagonal set of linear algebraic equations. These equations are solved by using alternate direction implicit method. Step 4. Repeat Steps 2 and 3 for the enthalpy equations (2.9). Step -5. With the newly calculated values at x + Ax of h z, and P (obtained by adding DP to the value of P at x), calculate p; using the equation of state. Next, calculate Y: and ZT using u:, and p; with the help of formulae given in Appendix A. Step 6. Check newly calculated u:, and hl, for convergence. If the convergence check is met, proceed to the next step; otherwise repeat the steps starting from Step 2. Step 7. Check Y: and Z;, calculated duct half width and half height for the value of DP assigned in Step 1, against the duct dimensions specified at x + Ax. If the check is within the desired accuracy, the computation of flow from x to x + Ax is complete. If not. DP is adjusted, on the basis of the difference between the specified and the calculated duct dimensions. and Step
1.
82
M.S. Greywall,
Stream wise computation
of parabolic flows
the calculations started again at Step I, with the new value of DP. The formulae used to adjust DP for the succeeding iteration is the one given in [I]. This formula, based on the difference between the calculated and specified duct cross-sections, is applicable to both two and three dimensional flow computations. Spacing between the streamlines in the cross-stream plane changes result, at the beginning of each new integration step, streamtubes accordance with the definitions given in Section 2.1.
from x to x + Ax. As a are defined anew in
4. Sample flow computations The approach, presented in this paper, to compute three-dimensional flows by calculating velocity and enthalpy along the streamlines, was checked by making three sample flow computations. All the computations are for laminar flow through rectangular ducts. Computation of laminar flows is free from the uncertainties associated with the modeling of turbulent shear stresses in three-dimensional flows and, thus, provides a clearer assessment of the numerical method itself. All the computations are for the lower-left quarter of the duct. Each sample flow computation was carried out using one hundred streamlines; this corresponds to, including the grid points along the walls, a 11 x 11 grid in the y-z-plane. Streamtubes around these streamlines were redefined at the beginning of each new integration step, as discussed in Section 3. In the present computations. however, the same set of streamlines was used for the entire length of the duct. The integration step along the duct, Ax, was selected for each sample flow computation so that the locations of the computer output will match the locations of the reported experimental data; the integration step varied from one-tenth to one-twentyfifth of the equivalent duct diameter. The effects of the integration step size and the number of streamlines used are discussed in Appendix B. Velocity distribution at the inlet in all the sample flows was taken to be uniform across the duct cross-section. The following nomenclature is used in the figures in this section. -.-._
u Ul UC W, H D X
YandZ Re Pr Gr Nu,
Dots along the velocity distribution curves represent the streamline which the flow was computed. axial velocity. U at the inlet, centerline velocity, duct half width and half height, respectively, equivalent diameter (= 4 WH/( W + H)), distance along the duct measured from the inlet, distances normal to x measured from the lower-left duct corner, Reynolds number (= U,,D/v, where v is the kinematic viscosity), Prandtl number (= 0.72) Graetz number (= Pr Re/x). Logarithmic mean Nusselt number - Tw,,, =$Gr log T,,,,,, T hulk
T hulk
bulk temperature.
-
Twa,,
>>
’
locations
alon
M.S. Greywall,
Streamwise
computation
83
of parabolic flows
In Figs. 4-7 the results for flow through a square-sectioned duct are shown. Calculations are compared with the experimental data of Goldstein and Kreid [13] (Figs. 4-6) and Beavers et al. [14] (Fig. 7). Agreement between the theoretical and experimental results is good. In Figs. 8-l 1 the results for flow through a rectangular duct of 2 : 1 aspect ratio are shown. Calculations are compared with the experimental data of Sparrow et al. [1.5]. Once again agreement between the theoretical and experimental results is satisfactory. In the two flow computations just described only the finite difference equations for velocity were solved. Density was considered constant and, thus, the energy equation was bypassed. To check the finite difference equations for enthalpy a sample flow with heat transfer to the walls through a rectangular duct of 2 : 1 aspect ratio was computed. The inlet temperature and velocity distributions were taken to be uniform across the duct cross-section. The walls were maintained at a constant temperature of 100 degrees Celsius below the inlet flow temperature.
Fully developed
0.02 >iO s E” ._
6
0.4
0.3
E &
o.2
a 5
0.1
work date Goldstein et al. (13)
o Expt.
1 0
0
2
1 0
2
1
2
u/u,
Fig. 3. Development
0.5
r
of the velocity profile in a square-sectioned
ductdiagonal
Fully developed
&=0.0075
-
plane.
3
Present work 0 Expt. data Goldstein et al. (13)
u/u, Fig. 5. Development
of the velocity profile in a square-sectioned
duct-central
plane.
84
MS.
Greywall. 2.2
Streamwise
computation
of parabolic flows
r
I
-
Present 0 Expt.
work
data
Goldstein
1.0 1
0
I 2
I
I
I
4
I
I
6
8
10
12
&
Fig. 6. Center
line velocity
et al. (13)
J
14
x102
development
in a square-sectioned
duct.
6
5
4 N >” z3
Present
:
Re 2 651 v 735 a 1222 01564 ,5 1858
2
Expt.
work
data
Beavers
et al
1
C
I
1
I 2
I
I
I
3
4
5
*
Fig. 7. Pressure
I 6
I
7
x lo2
drop in a square-sectioned
duct.
I
8
(14)
MS.
Greywall.
Stream wise computation
of parabolic flows
85
The velocity and enthalpy finite difference equations were solved together as discussed in Section 3. Calculations were carried out using constant specific heat, perfect gas behavior, and the Prandtl number equal to 0.72. Results of the computations are shown in Fig. 12, and are compared with the calculations of Wibulwas as reported in [ 16, Table 52, p. 2201. To illustrate the velocity field in the cross-stream plane, in Fig. 13 are plotted the streamline locations in the lower-left quarter of the duct for air flow through a $ meter by : meter rectangular duct (second sample flow computation) at two locations along the duct. The beginnings and the ends of the arrows represent, respectively, the streamline locations at 1 meter and 5 meters from the duct inlet. Bulk velocity through the duct is 0. IS cm/set.
0.0128
r
0.03881
Present work Expt. date Sparrow
Fig. 8. Development
of velocity profile in rectangular duct of 2 : 1 aspect ratio-across
et al. (15)
long side.
l.O-
NII ? t
0.6
-
0
k
_
s
Present work @ Expt. data Sparrow
0
0
0
1 0
1
2 1
2
et al. (15)
2
u/u,
Fig. 9. Development
of velocity profile in rectangular duct of 2 : 1 aspect ratb-across
short
side
86
M.S.
Grevwall,
Stream
of parabolic
wise cornputatiorl
jlows
Present @ Vertical 0 Horizontal
1.0.
I 3
I 2
I 1
0
traverse traverse
I 4 X D Re
x 10
Sparrow
et al. (15)
I
I 6
I 5
work
f=xpt, date
I 8
7
2
Fig. 10. Center line velocity development in rectangular duct of 3 : I aspect ratio.
8.0-
5.0
-
3.0
2.0
-
-
0
1
2
3 D
X
4 x I0
6
Present B
Re 897
Expt.
t
:ig
Sparrow
o A 0 m
1233 1416 1876 1986
6
work
data n+ .-1. . . ,_ _.
7
2
Fig. 11. Pressure drop in rectangular duct of 2 : 1 aspect ratio.
(15)
M.S. Greywall. Stream wise computation of parabolic flows
87
E 2’
5
2t 20
Fig. 12. Variation of logarithmic with constant wall temperature.
0.4-
Fig. 13. Motion
I
I
50
mean
Nusselt
number
**-e-e_,--+
of the streamlines
in the cross-section
-
- Wibulwas
I
100
I
150
with Graetz
--+
work
Present
--
200
number
(16)
J
220
for flow through
a rectangular
duct
--.+
plane-figure
is for the lower-left
quarter
of the duct
5. Concluding remarks A new approach to compute three-dimensional parabolic flows has been presented. The dependent variables in the present approach are the streamline velocity and the streamline coordinates in the cross-stream plane. Strictly speaking, the present approach is not Lagrangian, as one does not follow the individual fluid particles, but rather follows the flow through individual streamtubes. This approach, however, seems closer to the Lagrangian than the Eulerian flow description. The advantage of the present approach is that one needs to solve only one set of finite difference equations---for the streamwise velocity u-rather than three sets of equations for u,, u,, and u,. The disadvantage is that extra computations are required to calculate streamline coordinates.
h4.S. Greywall.
88
Stream wise computation
of parabolic flows
Appendix A In this appendix formulae which supplement the finite difference equations given in Section 2 are presented. The lower-left corner of the duct is taken as the origin of the y-z-plane. A. 1. Mass flow rates Let
Ri.j
=
&jUi.j.
Formulae
for stream
(2,2):
= :(2R,,, + Rs.2+ R,,)
Pi.;
. ;( Yi - YJ +;(Z, - Z,) .
‘z&f = y; . 1:(3R2.2+ Rz.x)( YJ . ;(Z.i - Z2) , ‘Pyt, = yf . +(3R,,, + R?.z) . +(Y3 - Y?)(Z,) . Formulae
for stream
(i, 2), i >2:
Formulae
for stream
(2, j),j > 2:
l~l,T = :(2Rz,j + R3.j + R,.j+l) . i( Y?-
Yz) . $(Z,+I - 4) 3
Ti,t = ~7 *i(3Rz.j + Rz,j+l)( Yz) . i(Zj+l - Zj) 3 lpp4 = ry *$(3Rz., + Rz,j-I)( YL) *$(Z, - Zjm1)3 Ti.7 = i(2Rz.j + Rz,j-1 + R,,j) . $( Y.3- J’z) . $(q - &I> Formulae
for stream IJZ-C+T = WLY+
A.2.
(i j). i and j > 2:
i(2Ri.j + Ri+l,j + Ri,j+l) * $(J’+I - X) . t(4+1-
~$(2Ri,,
+R,j+,+
Ri_I,j).$(X-
$(X - K-1) * %Zj- -T-l)
$(;?Ri,j + Ri_I,, + Ri,j-1) .
IJY~,T=
i(2Ri.j + Ri,j_1 + Ri+l,j) . t( X+I -
Surface and cross-sectional areas
$(Y3 -
Yz)
+
Y,)Ax
,
4)
3
X-I>~~(G+I-Z),
IJFL,T =
A;’ =
.
9
Yi> * k
M.S. Greywall.
Streamwise
computation
89
of parabolic flows
AsJ = $(y,,, - k;-,)Ax, A;'
= &Z, - Z, + Z&Lx,
A;'
= +(Zj+, - Z,-,)Ax,
A+
r.1
=-@‘A: hX2.
A.3. Streamline coordinates Y t, ZT After the velocities and densities have been calculated used to calculate the streamline coordinates at x + Ax. Let RT, = p:iU; and
For a given value of j. by summing i = I, one obtains for j = 2,
at x + Ax, the following
the mass flow rates through
the streamtubes
procedure
is
with i = 2 to
(A-1) and for j > 2, (Zf - ZT_,)Yf
By adding
(A.l)
16P’,; = 77(3Rl,i+ R;i_l)+
16P;,;_, rr-1(3Rl,ip1 + R&)+
I l6 z? Gi’i*
64.2)
to (A.2) for j = 3 to j = J, one gets 64.3)
z:y:=c, I
where c, represents the sum of the right-hand sides. By using the duct aspect ratio, specified at x + Ax, (A.3) can be solved for Y:. With Y: known, (A.l) and (A.2) can now be solved for Y G, Yf, . . . , Y: are calculated similarly. z;, z;. . . . ) Z:. The coordinates A.4. Mass, momentum and energy correction factors For laminar flow, the streamlines adjacent to the duct walls can be chosen near enough to the walls so that the velocity and enthalpy distributions between the walls and the adjacent streamlines can be approximated by linear functions. In this case V& and V:,i, introduced in Section 2.2, are given by the general expressions for V’s given in Section 2.1, and the
M.S. Greywall.
90
correction
factors
y,
Streamwise
computation
and 6 take on the following
a,
of parabolic
flows
values
These are the values for y, cy, and 6 used in the sample flow computations given in Section 4. The finite difference equations for the enthalpy were solved by using the wall temperature as the reference temperature; for this case the p’s become equal to the corresponding cy’s. For turbulent flow the correction factors y, (Y, 6 and p, and V$ and V:j will depend on the turbulence mode1 chosen for the near wall region in three-dimensional flows. Their calculation can be carried out by using three-dimensional extension of the procedure given in [l] for the case of two-dimensional flows. A._5. Viscous forces OHstream (i, j) The viscous forces exerted on the stream (i, j) by the neighboring streams are calculated as follows. The force exerted on the stream (i, j) by the stream (i, j + 1). call it F(top), is given by = Ap
F(top)
(a24
dz) .
(A.4)
where A is the interface area between the streams (i, j) and (i, j + l), and p is the effective viscosity at this interface. By substituting the appropriate symbols for A and ,u defined in Section 2.1 (see Fig. 3) and by approximating the velocity gradient by l&j+1
-
Ui.1
Zj+, -Zj one obtains
’
from (A.4)
F(top)
= AsJp{j+l/z ‘g,+’ I >I .
By using the definition F(top)
=
V:]+I(W,j+l
Similarly one obtains call it F(bottom), F(bottom)
of V:,,,,
(AS) kj)
is rewritten
= V:j(ni,j -
as
.
for the viscous force exerted
The net force on the stream by,
-
(A.5)
I
/+I
(A.6)
by the stream
(i. j - I),
(A.7)
u,,,_,).
(i, j) by the neighboring
(i, j) on the stream
streams
(i, j + 1) and (i. j - 1) is thus given
MS.
Greywall,
Streamwise
computation
of parabolic flows
91
V{j+l(Ui,j+lUi.j)Vi.j(Ui,j - Ui,j-l)e
F(top)-F(bottom)=
64.8)
By adding to (A.8) a similar expression for the net force exerted by the streams (i + 1. j) and (i - 1, j), one obtains the required relation for the viscous forces exerted on the stream (i, j) by the neighboring streams.
Appendix
B
In this appendix the effects of the integration step size along the duct axis, and the number of streamlines employed, are discussed. In Fig. 14 is shown the pressure drop in a squaredsection duct calculated with various integration step sizes and with various numbers of streamlines. The calculations were carried out for the lower-left corner of the duct, and the numbers given in Fig. 14 for the streamlines employed are for the quarter section. The streamlines were arranged in equal numbers of rows and columns (see, for example, Fig. 13). The spacing between the streamlines was non-uniform; the general nature of the spacing used can be seen from the velocity distribution figures given in Section 4. The integration step size, Ax, shown in Fig. 14 is normalized by D, the equivalent duct diameter. All the calculations were started with an initial integration step equal to 1/512th the final size. The initial integration step was then doubled at the end of each integration step until it reached its final value. The pressure drops for the three cases: Ax = O.lD and 100 streamlines; Ax = O.lD and 225 streamlines; Ax = 0.05D and 100 streamlines, were found to be very close (A&& = 6.478,
7_
6-
6-
Number of streamlines
AX/D 0.1 0.1 0.1 0.06
--0
1
2
3
4
6
-
Fig. 14. Effects of the integration sectioned duct.
X D Re
6
x10
0.2
100 7
8
6
2
step size and the number of streamlines
on the pressure drop in a square-
M.S. Greywall.
92
Streamwise
computation
of parabolic
flows
6.444, and 6.483 at x/DRe equal to 0.09 for the three cases, respectively). The difference was too small to be shown in the figure; thus, the three cases are shown as a single curve. The calculations were carried out using FORTRAN IV programming language on IBM 370/3031 computer operating under CMS Batch. The CPU times to carry out the computations, for duct length equal to 10 times the equivalent duct diameter, were as follows: Ax = O.lD;
25 streamlines, 100 streamlines, 225 streamlines,
100 streamlines;
Ax = 0.050, Ax = O.lD, Ax = 0.20.
CPU time = 1.28 mns, CPU time = 4.00 mns, CPU time = 10.60 mns, CPU time = 4.31 mns, CPU time = 4.00 mns, CPU time = 6.00 mns.
For a fixed value of Ax, the CPU time is scaled, as expected, by the number of streamlines used. For a fixed number of streamlines, however, the CPU time does not necessarily decrease when Ax is increased. As Ax was increased, although, the number of integration steps required to reach the end of the duct decreased, the number of iterations needed for each integration step increased.
Acknowledgment We would like to thank Dr. John M. Smith of NASA-Lewis Research Center for his valuable discussions and suggestions during the course of this work. The work was supported by NASA/DOE Grant NSG-3186.
References [l] M.S. Greywall. Streamwise computation of duct flows. Comput. Meths. Appl. Mech. Engrg. 21 (1980) 231-247. [2] L.S. Caretto, R.M. Curr and D.B. Spalding. Two numerical methods for three-dimensional boundary layers, Comput. Meths. Appl. Mech. Engrg. I (1972) 39-57. [3] S.V. Patankar and D.B. Spalding, A calculation procedure for heat, mass, and momentum transfer in three-dimensional parabolic flows, Internat. J. Heat and Mass Transfer 15 (1972) 1787-1806. [4] A.J. Chorin. A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1967) 12-26. [S] J.A. Miller, Laminar incompressible flow in the entrance region of ducts of arbitrary cross section, J. Engrg. Power Trans. ASME 93 (1971) 11>118. [6] A.F. Emery, P.K. Neighbors and F.B. Gessner, The numerical prediction of developing turbulent flow and heat transfer in a square duct, J. Heat Transfer Trans. ASME 102 (1980) 51-57. [7] D.W. Roberts and C.K. Forester. Parabolic procedure for flows in ducts with arbitrary cross sections, AIAA J. 17 (1979) 33-40. [8] V.S. Pratap and D.B. Spalding, Fluid flow and heat transfer in three-dimensional duct flows, Internat. J. Heat and Mass Transfer 19 (1976) 118>1188. [9] J.F. Nash, The calculation of three-dimensional turbulent boundary layers in incompressible flow, J. Fluid Mech. 37 (1969) 625-642. [ 101 K.C. Wang, Three-dimensional boundary layer near the plane of symmetry of a spheroid at incidence, J. Fluid Mech. 43 (1970) 187-209.
M.S. Greywall,
Streamwise
computation
of parabolic flows
93
[ 1 I] F.H. Harlow and J.E. Welsh, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Phys. Fluids 8 (1965) 2182-2189. [12] B.E. Launder and D.B. Spalding. The numerical computation of turbulent flows, Comput. Meths. Appl. Mech. Engrg. 3 (1974) 269-289. [I31 R.J. Goldstein and D.K. Kreid. Measurement of laminar flow development in a square duct using a Laser-Doppler flowmeter. J. Appl. Mech. Trans. ASME 34 (1967) 81.3818. [I41 G.S. Beavers, E.M. Sparrow. and R.A. Magnuson, Experiments on hydrodynamically developing flow in rectangular ducts of arbitrary aspect ratio. Internat. J. Heat and Mass Transfer 13 (1970) 689-702. [ 151 E.M. Sparrow, C.W. Hixon, and G. Shavit. Experiments on laminar flow development in rectangular ducts, J. Basis Engrg. Trans. ASME 89 (1967) 116-124. (161 R.K. Shah and A.L. London, Laminar Flow Forced Convection in Ducts (Academic Press. New York, 1978).