enhancement of flare method for predicting buoyancy

11 downloads 0 Views 598KB Size Report
Apr 25, 2007 - This article was downloaded by: [National Cheng Kung University] ..... Equation (11) becomes where. V;-l,i AY - 4/Re Pr. A i N = ZAy 2. U l' ..... C. H. Cheng, W. H. Huang, and H. S. Kou, Laminar Free Convection of the Mixing.
This article was downloaded by: [National Cheng Kung University] On: 23 May 2013, At: 01:53 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/unhb20

ENHANCEMENT OF FLARE METHOD FOR PREDICTING BUOYANCY-INDUCED FLOW REVERSAL IN VERTICAL DUCTS VIA PARABOLIC MODEL a

a

Chin-Hsiang Cheng , Shuen-Yi Huang & Win Aung

b

a

Department of Mechanical Engineering, Tatung Institute of Technology, 40 Chungshan North Road, Sec. 3, Taipei, Taiwan, 10451, Republic of China b

College of Engineering, Southern Illinois University, Carbondale, Illinois, USA Published online: 25 Apr 2007.

To cite this article: Chin-Hsiang Cheng , Shuen-Yi Huang & Win Aung (1997): ENHANCEMENT OF FLARE METHOD FOR PREDICTING BUOYANCY-INDUCED FLOW REVERSAL IN VERTICAL DUCTS VIA PARABOLIC MODEL, Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology, 31:3, 327-345 To link to this article: http://dx.doi.org/10.1080/10407799708915113

PLEASE SCROLL DOWN FOR ARTICLE Full terms and conditions of use: http://www.tandfonline.com/page/terms-and-conditions This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. The publisher does not give any warranty express or implied or make any representation that the contents will be complete or accurate or up to date. The accuracy of any instructions, formulae, and drug doses should be independently verified with primary sources. The publisher shall not be liable for any loss, actions, claims, proceedings, demand, or costs or damages whatsoever or howsoever caused arising directly or indirectly in connection with or arising out of the use of this material.

ENHANCEMENT OF FLARE METHOD FOR PREDICTING BUOYANCY-INDUCED FLOW REVERSAL IN VERTICAL DUCTS VIA PARABOLIC MODEL

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

Chin-Hsiang Cheng Shuen-Yi Huang Department of Mechanical Engineering, Tatung Institute of Technology, 40 Chungshan North Road, Sec. 3, Taipei; Taiwan 10451, Republic of China

Win Aung* College of Engineering, Southern Illinois University, Carbondale, Illinois, USA Enlumeemenl of 1M FLARE method ftwtunUJy tmp/oyed in a porobolit:~n modtl tuItl1ysis for buoyoney-illlhu:edflow rtvtrsal wilkin vtT1U:ol cluuultls is investigtlUd. ~ prtStnl study slwws Stveral modiji i; calculated along the negative y direction, where i, denotes the grids at the centerline. The computation starts from the inlet at X = 0, and then is advanced continuously in the axial direction. By following this axial marching procedure, the solutions at each cross section are obtained by using the known values at the previous cross section. The computation is terminated when the fully developed flow region is reached. The flow is regarded as fully developed when the relative changes of the axial velocity profiles between two consecutive cross sections are less than to- 4 at all grid points.

FLARE Method (B.L.

+ FLARE)

The marching procedure based on the above simple parabolic boundary-layer model generally fails when dealing with separated flows. The reasons for the divergence of the solutions may be physical or numerical. From the standpoint of numerical analysis [16J, a useful strategy to avoid divergence of iteration in solving a set of algebraic equations is to rearrange the equations so as to keep the coefficient of the largest magnitude on the diagonal at each step. This is called pivoting. However, it is observed that a negative U;-l.i may decrease the relative magnitude of the diagonal coefficients Ci P and A j p in Eqs. (14) and (15), respectively, or even make these values negative. These diagonal coefficients no longer dominate the matrix system of the algebraic equations, permitting roundoff errors to accumulate in the iterative back-substitution process used to solve the equations and resulting in numerical instability. As was described in [6, to, 11], the calculation can be advanced into the reversed flow region by introducing the FLARE approximation, which suggests the replacement of the convective term uau/ax in Eq. (9) with clulau/ax. This leads to a slightly different finite-difference form for Eq. (9): (18)

FLARE METIIOD FOR FLOW REVERSAL IN DUcrS

335

where

CiS

=

(Vi-l. j dY

+ 4/Re)

2dy 2

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

1

D.=I dX E i -

Clu;-1 jlu;-l j .. dX

+ Pi - 1 + -Gr -8·· 2Re 2 ,.}

j = 1,2, ... .n

where C = 0 for U(Y) < 0 and C = 1 for U(Y) > O. Note that by letting C be zero or a small value for the points with negative velocities, the value of C i P at these points will not be allowed to decrease so much as to produce the numerical instability. This tantamounts to adding an artificial convection to maintain the diagonal dominance of the algebraic system in the presence of reversed flow, which has also been noted by Pletcher [17]. However, the FlARE method still has its limitations. It will be illustrated later that the FlARE method is applicable only when the Richardson number is lower than a threshold value. Above the threshold value, a spatially oscillating solution may be observed downstream and the solution may diverge as the Richardson number becomes even higher. Improvement in the FlARE method is still necessary. Consistent FLARE Approximation (B.L. + FLAREC)

In the present study, a modified version of the FlARE method is proposed to enchance the boundary-layer model. In the FlARE procedure, the convection term Uau/ax of the momentum equation (9) is replaced with clulau/ax for solving the distribution of U. The values of V are obtained by solving the untouched continuity equation (8) and using the known values of U. That is, the approximation is introduced into the momentum equation, while keeping the continuity equation unchanged. This raises the question of consistency between these two equations. The inconsistency may cause local errors that are swept downstream and lead to the spatial oscillation observed in the numerical solutions. To introduce a consistent approximation that still yields suitably simple expression for V, the continuity equation (8) is multiplied by U to yield

au av + U - =0 ax ay

U-

(19)

336

C.·H. CHENG ET AI..

Since the term uau/ax is replaced by clulau/ax in the momentum equation (9), we now have the same treatment for the same term in this above equation, giving

au av Clul-+u-=o ax aY

(20)

There then folIows the finite-difference form for V:

V;,i =

CIU·I U".J (l!;,i-I

V;,i-t -

+ ti, -

~Y

l!;-I,i-I - l!;-I,i) 2

~X

(21a)

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

'oJ

for i < i. calculated along the positive y direction, and

V;,i =

V;,i+1

CIU·I oJ

+ U'

(l!;,i+1

+

~Y

l!;.i - l!;-I,i+1 - l!;-I,i)

2~X

(21b)

'oJ

for

i > ic calculated

along the negative y direction. The changes that are necessary to incorporate the present modifications for the continuity equation into the FLARE method are minor, but the enhancements are great. It wilI be demonstrated that the oscillation phenomenon and the divergence in the numerical solutions no longer appear for all the cases considered in the present study by using this modified method. RESULTS AND DISCUSSION

TypicalIy, a grid system of 1,601 x 41 grid points (in the X x Y directions) is distributed uniformly on the solution domain, which ranges from X = 0 to X = 40. That is, the grid sizes are ~X = ~Y = 0.025. It is important to note that the same grid sizes are used for each of the elIiptic and the parabolic models. Results of flow field and temperature distributions are displayed in this section with streamlines and isotherms, respectively. The solutions of temperature enable the heat transfer rate to be further evaluated. The dimensionless local heat fluxes on the respective walls are defined as

QI.

H

aTl -0

= - T1 - To ay

ae I

(22a)

ae I

(22b)

aY Y-O

and

aY Y_I Comparisons of Flow and Temperature Fields

Figure 3 shows the numerical solutions for flow fields based on the fulI Navier-Stokes model for various Richardson numbers at Re = 100.The occurrence

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

339

FLARE METHOD FOR FLOW REVERSAL IN DUcrS

Table 1. Comparison of solution methods, based on me data of separation point location (S.P.) and maximum stream function (Sma,)

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

D.L Re

Ri

100

2.0 3.0 4.0 4.45 4.85 4.9 5.0 5.25 5.3

200

300

500

S.P.

6.096 5.120 4.583 4.529 4.427 4.203 4.163

1.0 1.5 2.0 2.1 2.2 2.3 2.35

12.163 11.156 10.382 9.766 9.501

1.0 1.25 1.4 1.45 1.5

20.957 16.726 15.831 15.081

0.6 0.8 0.85 0.9

30.388 27.344 25.132

B.L. + FLARE Sma,

1.0 1.0 1.00671 1.01510 /;.

x x x x 1.0 1.0 1.00659 /;.

x

x

x 1.0 l.(lO306 /;.

x

x 1.0 1.00354

x x

S.P.

6.096 5.120 4.583 4.529 4.427 4.203 4.163

12.163 11.156 10.382 9.766 9.501

20.957 16.726 15.831 15.081

30.388 27.344 25.132

Sm.. 1.0 1.0 1.00671 1.01510 1.02489 1.02621 1.02887 1.03666 /;.

1.0 1.0 1.00659 1.00994 1.01384 1.01814 /;.

1.0 1.00306 1.00951 1.01172 /;.

1.0 1.00353 1.00642 /;.

B.L+FLAREC

Elliptic

Analytic ED.

Sm..

S.P.

Sm..

s,..

6.096 5.120 4.583 4.529 4.427 4.203 4.163

1.0 1.0 1.00671 1.01510 1.02489 1.02621 1.02887 1.03666 1.03823

14.164 5.568 4.784 4.335 4.288 4.201 4.008 3.973

1.0 1.00008 1.00754 1.01617 1.02612 1.02748 1.03026 1.03825 1.03984

1.0 1.00002 1.00709 1.01547 1.02535 1.02673 1.02958 1.03719 1.03880

12.163 11.156 10.382 9.766 9.501

1.0 1.0 1.00659 1.00994 1.01384 1.01814 1.02065

28.229 11.080 10.280 9.647 9.132 8.908

1.0 1.00008 1.00754 1.01096 1.01509 1.01940 1.02214

1.0 1.00002 1.00709 1.01043 1.01439 1.01892 1.02140

20.957 16.726 15.831 15.081

1.0 1.00306 1.00913 1.01169 1.01485

42.321 18.678 15.409 14.678 14.054

1.00008 1.00419 1.01096 1.01399 1.01725

1.00002 1.00383 1.01043 1.01334 1.01658

30.388 27.344 25.132

1.0 1.00350 1.00633 1.01007

57.915 27.689 25.253 23.422

1.00008 1.00752 1.01179 1.01724

1.00002 1.00709 1.01136 1.01658

S.P.

- . No flow reversal. /;.. Oscillating solution. x, Divergent solution.

Table 1 displays the solutions obtained by various solution methods, including B.L., B.L. + FLARE, B.L. + FLAREC, and the elliptic model, for the location of separation point (S.P.), as well as the value of maximum dimensionless stream function in the fully developed region (Smax ), at various Reynolds and Richardson numbers. The analytical solutions of Sma. deduced from the fully developed velocity profiles [2, 3] are also given for comparison. Relative performance of these boundary-layer models may be evaluated based on the results given in this table. The threshold values for the instability to occur are obviously dependent on the model in use. Considering the cases at Re '= 100, one finds that the flow reversal is detected as Ri ;;. 4.0. When the simple boundary-layer (B.L.) method is used, the oscillation takes place at Ri = 4.85. The application of FLARE can delay the oscillation, but the upper limit is still found at Ri = 5.3. Note that the oscillation

C.-H. CHENG ET AL.

340

and the divergence of the solution no longer appear in all the cases when the FLAREC method is adopted. Of particular interest is that, as compared with the analytical solutions; the relative errors associated with the elliptic and the B.L. + FLAREC models are nearly of equal magnitude. However, it is found in the present study that the computation time required for an elliptic-model analysis is generally about 100-200 times that for a B.L. + FLAREC analysis. In general, on a personal computer (486-33 MHz), it takes about 3-4 h to complete a run of the program based on the

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

4O ....

.....

_ ~ _ . . , . . .

....,

--.

Re-100

p

-_..... Elliptic Model ---.

30

B.L+FlARE

- - B.L+FLAREC

20

10

0

0

20

'0

30

40

xIH (a)

10

p

R..300

• •

0r--

_._._._._ ..

ElllpllcModel - - - - _. B.L+FLARE - - B.L+FLAREC

RI-1.45

-~

xIH

"(b)

Figure 6. Variation of pressure in axial direction, predicted by different solution methods: (a) Re = 100, Ri = 3 and 5; (b) Re = 300, Ri = 1.25 and 1.45.

FLARE METHOD FOR FWW REVERSAL IN

Q..

overs

341

10

, ,

&

Q..

- --- Q. - - a.

R......

R1_0.9

ElllpClcMod.

5 --0

~--

o

-. ---. --. ---. -----------

00

-. ----- ------------- ---B.L+FLAAE

5

-

0_0 000

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

o

----. ---. ------ -----. ----------. B.L.+A.AREC

~\

--

--0 000

-. -...

__ A.

o

-------- -------------- ------- --.- ----.

..

.0

80

Figure 8. Variation of local heat fluxes in axial direction, predicted by different solution methods, for Re = 500 and Ri - 0.9.

elliptic model; however, only 90-120 s are needed for a boundary-layer analysis. The executions of the programs on a workstation (llO-MHz MicroSPARC II) show a similar proportion of time reduction. Comparisons of Pressure Variation and Heat Transfer

The discrepancy in the solutions of pressure variation between the boundarylayer and the elliptic models is observed in Figure 6. In this figure, the data for

0..

5

&



Ir'----r----,-----,----....., •.•.•.•.• EUiptio Modlll

-

0 ••

~

- - - 8.L+flARE

- - B.L.FLAREC

2

0 ..

&

8

0 ..

o o

R90500 RI-o.85

10

20

so

xIH

.

Figure 7. Variation of local heat fluxes in axial direction, predicted by different solution metbods.

C.-H. CHENG ET AL.

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

342

Re = 100 at Ri = 3 and 5 and those for Re = 300 at Ri = 1.25 and 1.45 are demonstrated in Figures 60 and 6b, respectively. Note that the values of pressure on the centerline (Y = !) predicted by the elliptic model are taken to compare with the cross-sectional average pressure calculated by the boundary-layer models. The data show exactly the same trends and reasonable differences. Figure 7 shows the comparison in local heat fluxes (QI. and Q2.), which are calculated by Eq. (22). Good agreement is found among various solution methods. Figure 8 shows the oscillation of local heat fluxes produced by the B.L. + FLARE method for Re = 500 and Ri = 0.9. However, the oscillation can be avoided by using the elliptic or the B.L. + FLAREC models. To have a deeper insight into the axial variation of the local heat fluxes predicted by these models, Table 2 displays the data of Q1< and Q2. at different axial locations for this case. The B.L. + FLARE model fails to provide a stable solution for Q2. in the downstream region (x/H > 57), where the value of Q2. is oscillating between 1.60 and 0.2. However, B.L. + FLAREC can still lead to a stable solution that agrees closely with the elliptic-model results. The influence of the FLARE coefficient C on the numerical solutions has been evaluated. The value of C has been varied from - 0.35 to 0.35 to adjust the strength of the added artificial convection. It is found that the changes of C lead to a difference of less than 1% in the pressure and heat transfer data, for all the cases considered in this study. Table 3 shows the results of streamwise pressure variations predicted by B.L. + FLAREC method based on different values of C, for the case at Re = 300 and Ri = 1.45. Only a small discrepancy in the predicted results is observed when C is changed. This evaluation demonstrates an important fmding, that as long as the iteration remains stable, the effects of the convection in the

Table 2. Local heat fluxes (Q", Qz.) predicted by different solution methods, for Re = 500 and Ri = 0.9

Qz.

Q" B.L.+

B.L+

Elliptic

FlARE

FlAREC

Elliptic

0.20 1.00 4.00 8.00 15.00 20.00 40.00

13.31516 5.92061 3.51503 2.77816 2.22327 1.97582 1.41725

13.57855 6.15045 3.79715 3.03583 2.44236 2.17016 1.52674

13.57855 6.15045 3.79715. 3.03583 2.44236 2.17016 1.52377

50.00 59.00 60.00 64.00 70.00 75.00 80.00

1.27494 1.18925 1.18156 1.15388 1.12010 1.09774 0.99999

1.35484 1.23429 1.23381 1.19346 1.14152 1.03417 1.02108

1.35117 1.24547 1.23595 1.20148 1.13913 1.03083 1.00761

»rn

U=

oscillating solution

B.L+

B.L+

FlARE

FlAREC

0.0‫סס‬oo

0.0‫סס‬oo

0.0‫סס‬oo

0.0‫סס‬oo

0.0‫סס‬oo

0.0‫סס‬oo

0.00021 0.01638 0.14461 0.27165 0.67186

0.00015 0.01318 0.12669 0.24598 0.64548

0.00015 0.01318 0.12669 0.24598 0.65533

10.761271 10.265021 11.55925 1 10.20090 1 1 1.032041 10.277031 10.869101

0.77064 0.84095 0.84713 0.86948 0.89685 0.91516 0.99018

0.78360 0.85101 0.85704 0.87883 0.90539 0.92300 1.0‫סס‬oo

--

FLARE METHOD FOR FLOW REVERSAL IN

nucrs

343

Table 3. Pressure variations in streamwise direction predicted by B.L. + FLAREC using different values of coefficient C, for Re = 300 and Ri = 1.45

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

P

x/H

C = -0.35

C = -0.2

C =0.0

C =0.2

C = 0.35

0.10 0.50 5.00 10.00 16.00 18.00 25.00 35.00 45.00 55.00 80.00

-0.07101 -0.18579 -0.38543 -0.27469 0.27462 0.54945 1.79713 4.02100 6.54269 9.21026 16.13424

-0.07101 -0.18579 -0.38543 -0.27469 0.27462 0.54935 1.79551 4.01709 6.53712 9.20359 16.12640

-0.07101 -0.18579 -0.38543 -0.27469 0.27462 0.54921 1.79333 4.01184 6.52964 9.19466 16.11591

-0.07101 -0.18579 -0.38543 -0.27469 0.27462 0.54907 1.79115 4.00657 6.52214 9.18573 16.10544

-0.07101 -0.18579 -0.38543 -0.27469 0.27462 0.54897 1.78952 4.00258 6.51647 9.17899 16.09756

reversed flow areas are so small that the solutions are not even sensitive to the changes in C. It is important to mention that the concept of FLAREC is not limited to applications in two-dimensional problems. The FLAREC method can also be readily applied for a three-dimensional flow separation. For example, the parabolic model employed by Cheng et al. [6] for the analysis of three-dimensional flow reversal inside a vertical rectangular duct can be remarkably improved by incorporating the FLAREC method. However, more studies are required to apply FLAREC to other kinds of three-dimensional separated flows. CONCLUDING REMARKS

The improvement of the parabolic boundary-layer models for the analysis of the buoyancy-induced flow reversal in a vertical channel has been studied. Enhancement of the FLARE method frequently employed in a boundary-layer analysis for the flow separation is investigated. The comparison in the solutions between the boundary-layer and the full Navier-Stokes models shows that the boundary-layer model requires much less implementation and computation efforts hut still retains acceptable accuracy of the numerical solutions. In general, for the cases considered in the present study, the computation time required for performing a run of the elliptic-model analysis is about 100-200 times that needed for a boundary-layer analysis. Results also show that as the Richardson number is higher than a threshold number, the forward-marching procedure based on the B.L. + FLARE method leads to an oscillating or even a divergent solution. The present study modifies the B.L. + FLARE method by means of removing the inconsistency between the continuity and the momentum equations. The numerical instability can be effectively eliminated by using this modified method (B.L. + FLAREC). For example, considering channel flow at Re = 500 and Ri = 0.9, the B.L. + FLARE model leads to an oscillating solution in flow and temperature fields downstream; how-

344

C.-H. CHENG ET AL

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

ever, B.L. + FLAREC can still provide a stable solution that agrees closely with the elliptic-model results. It is found that the changes that are necessary to incorporate the FLAREC method into the boundary-layer model are minor, but the improvement is significant. Note that the FLAREC method may be applied for, but not limited to, analyses of two-dimensional buoyancy-assisted flow separation. Meanwhile, it is emphasized that the full elliptic Navier-Stokes model will always be one of the major tools providing accurate predictions for the separated flows. The contributions of this study are intended for those who are searching for an alternative that may prove efficient in solving the same problems.

REFERENCES 1. W. Aung and G. Worku, Developing Flow and Flow Reversal in a Vertical Channel with Asymmetric Wall Temperatures, ASME J. Heat Transfer, vol. 108, pp. 299-304,1986. 2. W. Aung and G. Worku, Theory of Fully Developed Combined Convection Including Flow Reversal, ASME J. Heat Transfer, vol. 108, pp. 455-458, 1986. 3. C. H. Cheng, H. S. Kou, and W. H. Huang, Flow Reversal and Heat Transfer of Fully Developed Mixed Convection in Vertical Channels, AL4A J. Thermophys. Heat Transfer, vol. 4, pp. 375-383, 1990. 4. C. H. Cheng and C. J. Weng, Flow Reversal of Combined Convection in a Vertical Rectangular Duct with Unequally Isothermal Walls, Int. Commun. Heat Mass Transfer, vol. 18, pp. 127-140, 1991. 5. C. Gau, K. A Yih, and W. Aung, Reversed Flow Structure and Heat Transfer Measurements for Buoyancy Assisted Convection in a Heated Vertical Duet, ASME J. Heat Transfer, vol. 114, pp. 928-935, 1992. 6. C. H. Cheng, C. J. Weng, and W. Aung, Buoyancy Effect on the Flow Reversal on Three-Dimensional Developing Flow in a Vertical Rectangular Duct-A Parabolic Model Solution, ASME J. Heat Transfer, vol. 117, pp. 238-241, 1995. 7. O. K. Kwon and R. H. Pletcher, Prediction of Incompressible Separated Boundary Layers Including Viscous-Inviscid Interaction, ASME J. Fluids Eng., vol. 101, pp. 466-472, 1979. 8. V. I. Vasilev, Computation of Separated Duct Flows Using the Boundary-Layer Equations, AL4A J., vol. 32, pp. 1191-1199, 1994. 9. M. A. I. EI-Shaarawi and A Sarhan, Free Convection Effects on the Developing Laminar Flow in Vertical Concentric Annuli, ASME J. Heat Transfer, vol. 102, pp. 617-622, 1980. 10. T. Cebeci, A. A Khattab, and R. Lamont, Combined Natural and Forced Convection in Vertical Ducts, Proc. 7th Int. Heat Transfer Conf., vol. 3, pp. 419-424, Munich, Germany, 1982. 11. D. B. Ingham, D. J. Keen, and P. J. Heggs, Flow in Vertical Channels with Asymmetric Wall Temperature and Including Situations Where Reverse Flow Occur, ASME J. Heat Transfer, vol. 110, pp. 910-917, 1988. 12. T. A. Reyhner and I. F1ugge-Lotz, The Interaction of a Shock Wave with a Laminar Boundary Layer, Int. J. Non-Linear Mech., vol. 3, pp. 173-199, 1968. 13. C. H. Cheng and J. J. Yang, Buoyancy-Induced Recirculation Bubbles and Heat Convection of Developing Flow in Vertical Channels with Fin Arrays, Int. J. Heat Fluid Flow, vol. 5, pp. 11-19, 1994.

FLARE MEfHOD FOR FLOW REVERSAL IN DUCTS

34S

Downloaded by [National Cheng Kung University] at 01:53 23 May 2013

14. H. Schlichting, Boundary Layer Theory, 7th ed., chap. 7, McGraw-Hill, New York, 1979. 15. C. H. Cheng, W. H. Huang, and H. S. Kou, Laminar Free Convection of the Mixing Flows in Vertical Channels, Numer. Heat Transfer, vol. 14, pp. 447-463, 1988. 16. C. F. Gerald, Applied Numerical Analysis, 2d ed., chap. 2, Addison-Wesley, New York, 1978. 17. R. H. Pletcher, Prediction of 1ncompressible Turbulent Separating Flow, ASME J. Fluids Eng., vol. 100, pp. 427-433, 1978.