Companion Computer Programs

N

N Companion z Computer Programs

10.1 I n t r o d u c t i o n In this chapter we describe several computer programs for calculating twodimensional laminar and turbulent incompressible flows. In Section 10.2 we first describe computer programs based on integral methods discussed in Chapter 7 and present sample calculations. The computer program in Section 10.3 is based on the differential method discussed in Chapter 8 and is applicable to both laminar and turbulent flows for a given external velocity distribution and transition location. In this section we also present sample calculations for an airfoil with the external velocity distribution obtained from the panel method discussed in Section 10.4. Sections 10.5 and 10.6 present computer programs for incompressible laminar and turbulent flows with heat transfer and for infinite-swept wing flows, respectively. Section 10.7 presents the same computer program in Section 10.3 except that it has the capability to start turbulent flow calculations with a specified initial velocity profile. Section 10.8 presents a computer program employing the SA model using the numerical procedure discussed in Problems 9.1 to 9.7, and Section 10.9 a computer program for a plane jet discussed in Problem 9.8. Finally Section 10.10 presents several subroutines discussed in problems in Chapters 8 and 9. All programs, including the three programs in Chapter 5, are given in the accompanying CD-ROM and are summarized in Appendix 10A.

10.2 Integral M e t h o d s In Chapter 7 we discussed integral methods for calculating heat and momentum transfer in two-dimensional and axisymmetric laminar and turbulent flows. In this section we describe F O R T R A N programs for them and present sample calculations.

362

10. Companion Computer Programs

10.2.1 T h w a i t e s ~ M e t h o d This method is applicable to both two-dimensional and axisymmetric laminar flows (see Problem 7.11) and has the following input requirements: NXT Total number of x-stations. KASE Flow index, 0 for two-dimensional flow, 1 for two-dimensional flow that starts as stagnation-point flow, and 2 for axisymmetric flow. KDIS Index for surface distance; 1 when surface distance is input, 0 when surface distance is calculated. UREF Reference velocity, Uref, feet per second or meters per second. BIGL Reference length, L, feet or meters. Kinematic viscosity, v, square feet per second or square meters per CNU second. Dimensionless chordwise or axial distance, x/L. If KDIS = 1, then X is X the surface distance, s. UE Dimensionless velocity, Ue/Ure f. Dimensionless two-dimensional body ordinate or body of revolution raR dius; r/L. Its output includes ~*, 0, H, cy and Ro = ueO/u together with Reynolds mumber, RS, based on surface distance and external velocity.

10.2.2 S m i t h - S p a l d i n g M e t h o d This method is for laminar boundary-layer flows with variable Ue but uniform surface temperature. Its input is similar to Thwaites' method and consists of NXT, KASE, KDIS, UREF, BIGL, CNU and PR. We again read in X, R and UE. The output includes X, S, UE and ST (Stanton number). 10.2.3 H e a d ' s M e t h o d This method is applicable to only two-dimensional incompressible turbulent flows. Its input consists of the specification of the external velocity distribution, Ue/U~, UE(I), as a function of surface distance x/L, X(I), with uoo denoting the reference freestream velocity and L a reference length. The initial conditions consist of a dimensionless momentum thickness, O/L, T(1), and shape factor H, H(1), at the first station. In addition, we specify a reference Reynolds number RL = uccL/~, RL and the total number of x-stations, NXT. In the code, the derivative of external velocity due/dx, DUEDX(I), is computed by using a three-point Lagrange-interpolation formula. The output includes X(I), UE(I),

10.3 Differential Methods with CS Model: Two-Dimensional Flows

DUEDX(I), T(I), H(I), DELST(I), Eq. (7.3.2).

363

5*/L, CF(I) and DELTA(I), 5, defined by

10.2.4 A m b r o k ' s M e t h o d This method is only for two-dimensional turbulent flows. Its input and output instructions are similar to Head's method.

10.3 Differential Two-Dimensional

Methods Flows

with

CS Model"

In this section we present the computer program discussed in Chapter 8 for two-dimensional incompressible laminar and turbulent flows. Its extension to flows with heat transfer is discussed in Section 10.5 and to infinite-swept wing flows in Section 10.6, to turbulent flows with specified initial velocity profile in Section 10.7, to turbulent flows employing the SA model in Section 10.8 and to a plane jet in Section 10.9. This computer program, called BLP2, and also described in [1], consists of a MAIN routine, which contains the logic of the computations, and seven subroutines: INPUT, IVPL, GROWTH, COEF3, SOLV3, EDDY and OUTPUT. The following subsections describe the function of each subroutine. 10.3.1 M A I N BLP2 solves the linearized form of the equations. Thus an iteration procedure in which the solution of Eq. (8.2.24) is obtained for successive estimates of the velocity profiles is needed with a subsequent need to check the convergence of the solutions. A convergence criterion based on v0 which corresponds to f~ is usually used and the iterations, which are generally quadratic for laminar flows, are stopped when 15v0(= DELV(1))I < al

(10.3.1)

with sl taken as 10 -5. For turbulent flows, due to the approximate linearization procedure used for the turbulent diffusion term, the rate of convergence is not quadratic and solutions are usually acceptable when the ratio of 15vo/vol is less than 0.02. With proper linearization, quadratic convergence of the solutions can be obtained as described in [1]. After the convergence of the solutions, the O U T P U T subroutine is called and the profiles F, U, V and B, which represent the variables fj, uj, vj and bj are shifted.

364

10. Companion Computer Programs

10.3.2 S u b r o u t i n e I N P U T This subroutine prepares data for boundary layer calculations. The data includes grid in ~- and ~/-directions, dimensionless pressure gradient m(~), mass transfer fw(~) parameter and calculation of the 7tr-term in the CS model which is given by Eqs. (5.3.18) and (5.3.19). The streamwise grid is generated by reading the values at ~i. In general, the ~-grid distribution depends on the variation of Ue with ~ so that rapid variations in external velocity distribution and the approach to separation require small A~-steps (kn). For laminar flows, it is often sufficient to use a uniform grid in the 7/-direction. A choice of transformed boundary-layer thicknesses r/e equal to 8 often ensures that the dimensionless slope of the velocity profile at the edge, f"(~Te), is sufficiently small (< 10 -3) and that approximately 61 grid points are adequate for most flows. For turbulent flows, however, a uniform grid is not satisfactory because the boundary-layer thickness r/e and the dimensionless wall shear parameter Vw ( - - f ~ ) are much larger in turbulent flows than laminar flows. Due to the rapid variation of the velocity profile close to the wall, it is necessary to make much smaller steps in r/close to the wall. The program uses a r/-grid which has the property that the ratio of lengths of any two adjacent intervals is a constant, that is, hj = K h j _ l and the distance to the j - t h line is given by

K J-1 7/j=hl K-1

j=

1,2,...,J

K>

1

(10.3.2)

There are two parameters: hi, the length of the first Ar/-step, and K, the variable grid parameter. The total number of points, J, is calculated from j = In[1 + ( / ( - 1)(r/e/hl)] + 1 (10.a.3) lnK In practice, it is common to choose hi [DETA(1)] and K (VGP) so that, for an assumed maximum value of r/e (ETAE), the number of j-points do not exceed the total number (NPT) specified in the code. For example, for r/e = 50 and hi = 0.01, the number of j-points depends on K. Figure 10.1 can help in the selection of K and shows that, for example for J = 61, the value of the variable grid parameter must be less than about 1.1. With velocities known at each C-station, called NX-stations, the pressure gradient parameters m (P2) and m l (P1) are computed from their definitions for all C-stations, except the first one at which m is specified, after the derivative of due/dr (DUDS) needed in the calculation of m is obtained by using threepoint Lagrange interpolation formulas given by (n < N): n-1 n {due~ _-Ue (~n+l -- ~n)-t- -~---(~n+l ue -- 2 ~ n -I- ~ n - 1 ) \ d~ un+lA1 ~2 (10.3.4) +

A3 ( ~ n - ~ n - 1 )

10.3 Differential Methods with CS Model: Two-DimensionM Flows

365

1.24 1.20 J = 30

1.16

40

K 50

60

1.12 1.08 1.04

1.00

'

'

~

,'..'.~,..-'L-

0.2

. . . . . . . .

1

,

. . . . . . . .

,

. . . . . . . .

100

10

,

1000

10000

vlJhl x 10.2 F i g . 1 0 . 1 . V a r i a t i o n of K w i t h hi for different r/e-values.

Here N refers to the last ~n station and A: = (~n - ~ n - 1 ) ( ~ n + l

- ~n-1)

A 2 -" ( : n - ~ n - 1 ) ( ~ n + l

- ~n)

(10.3.5)

A3 = (~n+l ~n)(~n+l ~n--l) -

-

-

-

The derivative of d u e / d ~ at the end point n = N is given by -1

d~ IN

A1

(~N - ~ N - 1) +

-4- ~ ( 2 ~ N ~3

-- ~ N - 2

A2

..... (~N - ~ N - 2 ) (~o.3.6)

-- ~ N - 1 )

where now A1 = (~g-1 -- ~ Y - 2 ) ( ~ Y

-- ~Y-2)

A2 = (~N-1 -- ~ g - 2 ) ( ~ g -- ~ g - 1 ) A 3 = ( ~ Y -- ~ Y - 1 ) ( ~ g

(10.3.7)

-- ~N-2)

The 7tr-term, which accounts for the transition region between a laminar and turbulent flow, is calculated as a function of ~ once the onset of transition is specified. In this subroutine we specify Ue at ~ - 0 and the reference Reynolds number R L (RL). In addition, the following d a t a are also read in and the total number of j-points J ( N P ) is computed from Eq. (10.3.2). NXT NTR

Total number of x-stations NX-station for transition location Xtr

366

10. Companion Computer Programs

NPT Total number of r/-grid points. DETA(1) Ar/-initial step size of the variable grid system. Use A t / = 0.01 for turbulent flows. If desired, it may be changed. ETAE Transformed boundary-layer thickness, r/e VGP K is the variable-grid parameter. Use K = 1.0 for laminar flow and K = 1.12 for turbulent flow. For a flow consisting of both laminar and turbulent flows, use K = 1.12. Reynolds number, u~L RL v x Surface distance, feet or meters, or dimensionless. Velocity, feet per second or meter per second, or dimensionless. Ue 10.3.3 S u b r o u t i n e I V P L At x = 0 with b = 1, Eq. (8.2.6) reduces to the Falkner-Skan similarity equation which can be solved subject to the boundary conditions of Eq. (8.2.7). Since the equations are solved in linearized form, initial estimates of f j , uj and vj are needed in order to obtain the solutions of the nonlinear Falkner-Skan equation. Various expressions can be used for this purpose. Since Newton's method is used, however, it is useful to provide as good an estimate as is possible and an expression of the form. 3rlj 1 (rlj) 3 -(10.3.8) U j - - 2 tie 2 tie -

-

usually satisfies this requirement. The above equation is obtained by assuming a third-order polynomial of the form f ' = a + br/+ Cr]3 and by determining constants a, b, c from Eq. (8.2.7) for the zero-mass transfer case momentum equation which requires that f " The other profiles fj, vj follow from Eq. =

--

vj = 2 r/e

the boundary conditions given by and from one of the properties of = 0 at r / = r/e. (10.3.4) and can be written as

3-

(10.3.9)

~ee

(10.3.10)

10.3.4 S u b r o u t i n e G R O W T H For most laminar-boundary-layer flows the transformed boundary-layer thickness r/e(x) is almost constant. A value of r/e = 8 is sufficient. However, for turbulent boundary-layers, r/e(X) generally increases with increasing x. An estimate of r/e(x) is determined by the following procedure.

10.3 Differential Methods with CS Model: Two-Dimensional Flows

367

We always require that ~Te(Xn) >_ rle(xn-1), and in fact the calculations start with ~ e ( 0 ) = r/e(Xi). When the computations on x = x n (for any n > 1) have been completed, we test to see if Iv~l _~ ev at ~e(X n) where, say ev = 5 • 10 -4. This test is done in MAIN. If this test is satisfied, we set ~%(x n+~) = ~%(xn). Otherwise, we call G R O W T H and set Jnew = Jold + t, where t is a number of points, say t = 1. In this case we also specify values of (f~, u2, v2, b2) for the new nj points. We take the values of uj = 1, vr~ = O, f~ = (rlj - rle)urJ + fr~, and bjn = b~. This is also done for the values of f ~ - l , vr~-I and b~ -1. 10.3.5 S u b r o u t i n e C O E F 3 This is one of the most important subroutines of BLP2. It defines the coefficients of the linearized momentum equation given by Eqs. (8.2.20) and (8.2.23). 10.3.6 S u b r o u t i n e E D D Y This subroutine contains the CS algebraic eddy viscosity model in Section 5.8. For simplicity we do not include the low Reynolds number effect, roughness effect, mass transfer effect and strong pressure gradient effect (variable c~) in this subroutine. These capabilities, if desired, can easily incorporated into the formula as defined in this subroutine. The formulas for the inner and outer eddy-viscosity expressions are given by Eqs. (5.2.11) with each side of equation multiplied by ")'tr given by Eqs. (5.3.18) and (5.3.19). In terms of transformed variables (c+)i and (S+)o given by Eq. (5.2.11) can be written as (e+)i = 0.16R1/211 - exp(-y/A)]2~2VVtr

(10.3.11a)

(s+)0 = 0.0168R1/20%- fe)~tr')'

(10.3.11b)

y

N 171/4 1/2

SeX

(10.3.12)

10.3.7 S u b r o u t i n e S O L V 3 This subroutine is used to obtain the solution of Eq. (8.2.24) with the blockelimination method discussed in subsection 8.2.3 and with the recursion formulas given in subsection 8.2.4. 10.3.8 S u b r o u t i n e O U T P U T This subroutine prints out the desired profiles such as fj, uj, vj, and bj as functions of r/j. It also computes the boundary-layer parameters, c/, 6", 0 and Rx.

368

10.4 Panel

10. Companion Computer Programs

Method

The solution of the boundary-layer equations requires the specification of the external velocity distribution which may be provided either from experiments or from inviscid flow theory. In the latter case, a practical inviscid method is a panel method. Since some of the problems in this book require external velocity distribution, a computer program based on the Hess-Smith panel method (HSPM) is given in the accompanying CD-ROM. A complete description of this method can be found in [2]. The input of this program consists of airfoil coordinates (x/c, y/c), angle of attack, a, and freestream Mach number, Mc~, with compressibility corrections. The output includes lift and pitching moment coefficients together with pressure coefficient Cp defined by p- p~ cp= !2 0v2 and external velocity distribution ue/vc~. The latter is provided separately for the upper and lower surfaces of the airfoil. To illustrate the use of this code, in the accompanying CD-ROM we present results for the NACA 0012 airfoil. The calculations are performed for 184 airfoil coordinates with the x/c and y/c values read in starting on the lower surface trailing edge (TE), traversing clockwise around the nose of the airfoil to the upper surface TE. Results are given for a = 0 ~ 2~ 4 ~ with Mc~ = 0.1.

10.5 Differential Method with CS Model: Two-Dimensional Flows with Heat Transfer The program which for convenience we called BLP2H is the same computer program BLP2 which now includes the solution of the energy equation. Two subroutines, COEF2 and SOLV2, are added to BLP2 to calculate incompressible laminar and turbulent flows with heat transfer (see Problem 8.3). Sample calculations presented in the accompanying CD-ROM represent the application of this code to Problems 8.4 and 8.5.

10.6 D i f f e r e n t i a l M e t h o d Swept-Wing Flows

with CS Model:

Infinite

This program, called BLP2ISW, is also the extension of BLP2 to the calculation of infinite swept-wing equations for incompressible laminar and turbulent flows as discussed in Problem 8.6. Again two subroutines are added to BLP2 and changes are made to the eddy viscosity subroutines. Subroutine COEF2 includes

10.7 Differential Method with CS Model: Turbulent Flows

369

the coefficients of the z-momentum equation and subroutine SOLV2 is the same solution algorithm used in BLP2H and given in Section 10.10. For three-dimensional turbulent flows, the eddy viscosity formulas require changes to those for two-dimensional flows. Here they are defined according to Eqs. (5.7.4)and (5.7.5). Sample calculations for an infinite swept wing having the NACA 0012 airfoil cross section with a sweep angle of ~ = 30 ~ an angle of attack of c~ = 2 ~ chord Reynolds number Rc - 5 • 106 and transition location at x/c = 0.10 are presented in the accompanying CD-ROM. See also Problem 8.7.

10.7 Differential Method with with Initial Velocity Profile

CS Model:

Turbulent

Flows

See Problem 8.9.

10.8 Differential

Method

with

SA Model

This computer program called BLPSA is the extension of BLP2 with the procedure described in Problems 9.1 to 9.7. Many of the subroutines used in BLP2 remain the same except for some minor changes. Several additional subroutines are added. Subroutine COEF contains the coefficients of linearized continuity, momentum and eddy viscosity equations (Problem 9.4), subroutine MSA (subsection 10.10.3) is the solution algorithm to solve the linear system in Problem 9.4. The accompanying CD-ROM has two test cases for this program (Appendix 10A) and also includes another program for this model in which the continuity, momentum and eddy transport equations are all solved together. This program is referred to as 5 • 5 in contrast to the other one which is referred to as 2 • 2.

10.9 Differential

Method

for a Plane

Jet

See Problem 9.8.

10.10

Useful

Subroutines

In this section we present three subroutines that are useful to solve some of the problems in Chapters 8 and 9. They are briefly described below and are given in the accompanying CD-ROM.

370

10. Companion Computer Programs

10.10.1 S u b r o u t i n e I V P T This subroutine is for generating initial velocity profiles with the method discussed in Problem 8.9. It requires the initial values of Reynolds number based on momentum thickness Re (-- ueO/~') and local skin friction coefficient c/

(_--

10.10.2 S u b r o u t i n e SOLV2 This subroutine is similar to the solution algorithm, SOLV3, in subsection 8.2.4. It is designed to solve two first-order equations with the block-elimination method subject to the boundary conditions given by Eq. (P8.2.4). See also Problem 8.2. 10.10.3 S u b r o u t i n e M S A When the system of first-order equations to be solved with the Box method becomes higher, the preparation of the solution algorithm with the recursion formula becomes tedious [1]. A matrix-solver algorithm (MSA) described in [1] can be used to perform the matrix operations required in the block-elimination method. It consists of a subroutine where we read in matrices Aj, Bj, Cj appearing in the coefficient algorithm. Subroutines GAUSS, USOLV and GAUSV then perform the matrix operations in the block-elimination method. In this subroutine IROW, ICOL correspond to number of maximum rows and columns respectively. ISROW denotes the number of "wall" boundary conditions and INP the total number of j-points in the u-direction, and BB=Bj,

YY=wj,

GAMJ=Fj,

AA=Aj,

CC=Cj

The first and second numbers in the arguments of AA, BB, CC and GAMJ correspond to the number of nonzero rows and columns in Aj (or Aj), Bj, Cj and Fj matrices, respectively (see subsection 8.2.2, for example). Note that Bj and Fj have the same structure and the last row of Bj and the-first two rows of Cj are all zero. To use this subroutine, we first define the matrices Aj, Bj and Cj by reading in all their elements, and follow the steps given below. 9 Call subroutine GAUSS. 9 Call subroutine GAMSV to compute F1. 9 Compute Aj according to Eq. (8.2.27a), call GAUSS. 9 Call GAMSV to compute/"2. 9 Repeat the above two steps for j < J. 9 Compute w0 according to Eq. (8.2.29a).

References

371

9 Define the right-hand side of Eq. (8.2.29b) and compute wj according to Eq. (8.2.29b). 9 In the backward sweep, with 5j corresponding to UM(I,J), compute ~ j according to Eq. (8.2.30a) by calling USOLV at j -- J. 9 Define the right-hand side of Eq. (8.2.30b) and solve for 5j by calling USOLV for j = J - J , J - 2 , . . . , 0 .

References [1] Cebeci, T. and Cousteix, J.: Modeling and Computation of Boundary-Layer Flows. Horizons Pub., Long Beach, Calif., and Springer, Heidelberg, Germany, 1998. [2] Cebeci, T.: An Engineering Approach to the Calculation of Aerodynamic Flows. Horizons Pub., Long Beach, Calif., and Springer, Heidelberg, Germany, 1999.