Iy‘k

-." 7A.“.

 

" ~ “-u’nb-X 1"!
" ‘9
‘ h v 1 R
was)?“ ‘
,4

‘n
1w

~*jfl2r
.242”

4" 5'4

.1
,w
n

in :"RJ‘. '
t “may
1." V. :1;le :5:
' V
.

‘ my - ‘.
“Maw v‘N'ci-‘IT' gr
mu, .x'.‘

3* 43v

Bali“; ,
.72..» .

u
m

f

..
’2'
.L

“as

.-:?

‘~ 9957‘»:
"‘"21 “4‘71;-
" an'

 

:8

..
,,
fi-

1‘.

. v , a
. ' c‘
y, ’

v ,<

gt:
*3;

V
’5

 

:fit“;
I F.‘
A ,. r

's

“3:52:

i:
' t'
1;, , b

« ‘2
.1 4‘ 4 ~ 4 V t ‘
A 1" v" ‘ ' . I’ ‘ QIcS
‘52: IU’ "(351,117 ‘
a
‘i
44:.

\‘4 ’
J

Lied

3?,

I.
K " fig!) '3
r. ' N,“ h: 5..
Juvigfi};;ggn.
8' " ”It

a v. A

J"; ' 1‘»!

r v r p- .

’i'ta’aaifixvfixv

'J'J’ ‘3' ‘51 “

’ 3408'“ 114i.-
‘ n 1

‘40:”. 4 ,.-‘
"a film)“
I 4 ‘ VI“: ‘
a 1:15; ,1?!
"A ‘4'” ”U ”3H;
~ . ,. ., , N
’ ." ; (eififejsxe sate-:2.

25,814.49,“ (19’7"

12:13:}; " ’p

' 9r"
('1‘ figuu.
"’~”"I"§H

\

.4 . ‘.
I If
1.;

 

. s '0‘; ‘. ,

Ignaz; ' , am
I .
If I

. l ‘
‘-‘z1~-..’.- 1.1.4.1»
F1325” HAW”
r‘;:;g>3€7f-T?‘£€;gd '
r'lfrr‘réfiVW‘L ‘

a My;

rg»:
- Q .
v

a, .

' my. ‘1‘ ,
Effigy ..

x‘wzw

‘ ' .(4
l ll‘lf
’~' \‘Pt'f'
, ' l ..
(:3Qhé
l; ‘1‘,

>,». I. i _ <
{tn-31.9; a :3
1' 7- 1}” 5‘»: l 4

a?“ ”a ”w, «
Q‘AR‘EET “t.

"‘1'. Pf.“

w
3'

«2&1?
‘ “i; J;

a war g.“ ,u 74,313
' "’1‘"; “(2)1121 1;! a f}.
£1.39?” fa}; gag-y.
’“ ' ”“4 ":2"-
,,
,

.. A
\.

4&9?
a'l
\f/

4
.
L .
. I

‘ ‘ -‘ . .u 4". an. “If"
" , r ., .‘ *g:::¢, ;; . AII» , , ha
}l.‘){or \' 2‘" ’, x ‘ N

r ,-v

Map.
A," .

a
(Ll

{ » nlJ
(a pr 92.“

.

M
‘0‘

 

 

 

 

 

 

 

Date

gnqfing
llllHllHlHHllHll|U||l||||IIHNIH”IIIIIIIIIIIIIHHHHI

3 1293 00542 0108

 

imam?
Michigan State
University

 

 

 

This is to certify that the

thesis entitled
NUMERICAL STUDY OF NATQRAL CONVECTION
BETWEEN TWO VERTICAL PLATES WITH ONE
OSCILLATING SURFACE TEMPERATURE
presented by

WEI CHA

has been accepted towards fulfillment
of the requirements for

MASTER . SCIENCE
degree in

 

 

Major professor/

9%: L/ Vigil/OK
67

JAN. 6‘ 1989

V

 

0-7 639

MS U is an Affirmative Action/Equal Opportunity Institution

 

 

 

MSU

LIBRARIES
“

\—

 

 

RETURNING MATERIALS:
Piace in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

 

 

NUMERICAL STUDY OF NATURAL CONVECTION
BETWEEN TWO VERTICAL PARALLEL PLATES
WITH ONE OSCILLATING SURFACE TEMPERATURE

by
Wei Cha

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirement

for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering

1989

it) '0

fi"
\

x-

ABSTRACT

NUMERICAL STUDY OF NATURAL CONVECTION
BETWEEN TWO VERTICAL PARALLEL PLATES
WITH ONE OSCILLATING SURFACE TEMPERATURE

by

Wei Cha

This study utilizes a finite difference numerical method (employing the program
NUDSFAE) to simulate the natural convection heat transfer between two vertical
parallel plates, one of which has a time oscillation of its surface temperature. The nu-
merical code deals with real variables rather than stream function and vorticity trans-
formation variables and it solves the coupled equations by introducing a pressure cor-
rection scheme. The active surface temperature has a periodic time variation with
non-zero average. The dimensionless time average heat transfer rate, Nusselt num-
ber, is compared with that of a constant surface temperature case, and this constant
temperature is the same as the average of the oscillation. The results of this study
show that by oscillating the surface temperature the heat transfer rate can be in-
creased significantly. The study also investigates the effects of oscillation frequency,
amplitude and spacing aspect ratio on the enhancement of heat transfer.

The results and the study method can be applied to investigate new electronic

cooling techniques with natural convection heat transfer.

ACKNOWLEDGMENTS

I would like to thank my advisor, Dr. J. R. Lloyd, who helped me to initiate the
graduate study program, and whose advice and encourages not only helped me
through my studies at Michigan State University in the last two and half years but
will also benefit me for my whole carer life. I would also like to express my respect
and appreciation to Dr. K. T. Yang for his valuable advice. I wish to acknowledge Dr.
J. V. Beck and Dr. C. W. Somerton for their constructive comments to my thesis. I
feel very graceful to Dr. J. V. Beck, Dr. J. Foss Dr. C. W. Somerton, and Dr. Rosen-
berg for knowledge I learned from them. I would like to express great appreciations
to Dr. St. Clair for his efforts that make my studying here very enjoyable. I am very
gratitude to Dr. H. Q. Yang for his advising me on numerical method. I also like to
thank National Science Foundation for the financial support at the beginning of the nu-
merical study.

Finally, I wish to thank my husband, W. J. Gesang, for his being supportive

and understanding.

TABLE OF CONTENTS

LIST OF TABLES ...................................................................................................
LIST OF FIGURES .................................................................................................
NOMENCLATURE .................................................................................................
1. INTRODUCTION ................................................................................................
2. LITERATURE REVIEW .....................................................................................
3. NUMERICAL FORMULATION ........................................................................

3.1 Governing equations ....................................................................................

3.2 Numerical solution .......................................................................................
3.2.1 Cell structure .......................................................................................
3.2.2 QUICK scheme ....................................................................................
3.2.3 Pressure correction scheme .................................................................
3.2.3 Boundary conditions ............................................................................

4. RESULTS AND DISCUSSIONS ........................................................................

4.1 Validation test ...............................................................................................
4.2 Discretization of time step .........................................................................
4.3 Mechanism of natural convection heat transfer with an

oscillating boundary tmeperature .................................................................
4.4 Effect of oscillation frequency on Nusselt number .....................................
4.5 Effect of oscillation amplitude on Nusselt number ......................................
4.6 Effect of aspect ratio on Nusselt number ....................................................

5. CONCLUSIONS ..................................................................................................
6. RECOMMENDATIONS .....................................................................................
APPENDIX (Numerical code) ................................................................................
LIST OF REFERENCES ........................................................................................

ii

Page
iii

iv

v

1

6

1 1

11
13
13
15
16
18

22

22
23

27
46
50
50

54
55
56
81

LIST OF TABLES

Page
Table 1 Validity Test .................................................................................... 23
Table 2 Stability Test .................................................................................... 25

iii

LIST OF FIGURES

Page

Figure 1 The Importance of Heat Transfer in Electronic Design .................. 2
Figure 2 Problem Description ......................................................................... 5
Figure 3 Basic and Staggered Grids .............................................................. 14
Figure 4 Flow Chart for Pressure Correction ................................................ 17
Figure 5 Calculation Domain .......................................................................... 19
Figure 6 Surface Temperature as a Function of Time ................................... 21
Figure 7 Grid Size Effect on Accuracy ........................................................... 24
Figure 3 Time Step Size Effect on Accuracy .................................................. 26
Figure 9 Temperature and Flow Fields Responses to

Oscillating Surface Temperature ............................................ 29
Figure 10 Nusselt Number Variation with Time ............................................. 47
Figure 11 Frequency Effect on Heat Transfer Enhancement ......................... 49
Figure 12 Amplitude Effect on Heat Transfer Enhancement .......................... 51
Figure 13 Amplitude Effect on Heat Transfer Enhancement .......................... 53

iv

ASP

Ra

NOMENCLATURE

space aspect ratio

specific heat
distance between two vertical plates
boundary temperature variation frequency

gravitational acceleration

Grashof number [gBH3(T-T°°)/V2]
convective heat transfer cooefficient
height of the vertical plate

thermal conductivity

measurement of heat transfer enhancement [Nuo/Nuc]

Nusselt number [hH/k]
numericaloverflow

pressure field
Prandtl number [v/a]
Rayleigh number [GrPr=gB(T-T°°)H3/v0t]

numerically stable

dimensional time

ave,t

media characteristic time [Hz/v]

temperature
vertical velocity

nonfdimensional vertical velocity [u/uO]
characteristic velocity

horizontal velocity

non-dimensional horizontal velocity [v/uo]

horizontal coordinate

non-dimensional horizontal coordinate [x/H]

vertical coordinate

non—dimensional vertical coordinate [y/H]
Greek letters

thermal diffusivity [pCP/K]
volume thermal expansion coefficient [~pp/ T]

density
non-dimensional time [tv/Hz]
non-dimensional temperature [(T-TW)/(Tavc-T°°)]
Subscript
space average

time average
constant surface temperature case

oscillating surface temperature case

ambient

vi

1. INTRODUCTION

The increasing importance of thermal aspects of electronic equipment design is
one of the many motivations for continued interest in natural convection heat transfer
studies. It is a well known fact that electronics performance is strongly affected by
working temperature. To avoid malfunctions of electronic devices and prevent them
from being destroyed (burn out of the element or bad connection of the solder joints
due the thermal stress, etc.), the geometry of the cooling passages is a primary con-
sideration in packaging of the electric devices. On the other hand, the tendency of mi-
crominiaturization of electronic components and higher performance speed, make the
electronic design more temperature dependent. Therefore, electronic cooling problems
are drawing more attention. Nakayama [1] addressed the situation of thermal man-
agement of electronic equipment based on Hitachi Ltd’s technology and research on
thermal design of their electronic devices. In Figure 1, the thermal design is seen to
be of equivalent importance along with geometrical, environmental and economical
considerations. Due to the fact that natural convection cooling requires the least addi-
tional hardware, and for most cases the surrounding media, such as air, exist already,
natural convection is a primary mechanism for electronics cooling.

From a heat transfer point of view, basic natural convection heat transfer phe-
nomena should be studied to enable it to integrate into electronic design. Usually, the
features of interest can be classified by their geometry, their thermal properties, and

their boundary conditions. More precisely, there are channel natural convection and

5.89 u_=e...uo_m E 3.25:. .8: he 8:8..35— 2:. — «Eur.—

 

_ cozoanoi _

 

_l.| llllllllllllllllllll J

952 «o: “

 

 

32555.35 .
zoaom eozaom .moEoEooG .

83858.2. .meafi .
3:388. 55.5 5.89

_

q

“ _3_Eocoom .
.

 

 

 

 

 

 

-.
_. ....... + ............ L.

_ .oeam mocmpeorod oncoaomm _

 

 

3
single surface natural convection geometries; there are some situations where tem-

perature variations are not too big to make the Bousinesq approximation while there
are other cases that require the involvement of variable properties; there are cases
where the boundary-temperature is better known and cases where surface heat flux
is better known; and there are cases where the boundary conditions will vary with
time and space.

In electronic devices, there are usually many individual chips mounted on
boards which are placed next to each other. If the boards are close enough, channel
type natural convection flows will occur. If the boards are far enough apart from each
other, a single surface assumption is reasonable. Even for a single surface, the sur-
face curvature and heating locations can have a significant effect on the total heat
transfer.

For a single surface, depending on the surface physical properties and the func-
tion it has in the device, more may be known about the surface temperature than the
heat flux, or vice versa. For instance, if an electrically heated plate has a low thermal
initia, the heat flux may be a better known function of time than the temperature, but if
the thermal conductivity is very high, the surface temperature may be known as a
function of time. Example of the cases where surface temperature or heat flux varies
with time are found when the electric circuit switch for a device is turned on or off.
The thermal boundary conditions will have an exponential-like response to that step
change in current. More kinds of time variations can be from a periodically varying
power supply.

Boundary conditions can have a spacial variation. When individual chips are
separated sufficiently, the surface will have spots of which temperatures and heat

fluxes are higher or lower than the rest part of the surface.

4
The current study develops a numerical finite difference program to simulate

natural convection driven by electronic element, of which the surface temperatures
vary with time. The physical model of the problem studied is sketched in figure 2.
The physical model includes two vertical plates, with height H. They are placed paral-
lel to each other, separated by a distance D. Each of the two surface temperatures
can have its own function of time. The medium is air. An experiment was done pre-
liminarily to develop a general feeling for temperature and velocity fields caused by a
hot plate, especially with the time varying boundary conditions. Then a finite differ-
ence numerical method was utilized to reveal the detailed quantitative descriptions of
the temperature and velocity fields with one surface temperature varying with time
periodically. By doing this, the natural convection heat transfer and corresponding flu-
id flow driven by an oscillating surface temperature is simulated numerically.

The success in simulating the natural convection driven by the time varying
wall temperatures, provides the way to expose the magnitude of heat transfer en-
hancement achieved by oscillating the surface temperature and to develop possible

approaches to thermal control in electronic thermal design.

'—

 

 

l/llllllll/l/Ill/l’l/l/l/l/Il/’l/l/l/l/l/l/l/l/I/Illli l

Tw2(t)=const.

  

Twl (t)

Figure 2 Problem description

 

.' ilfa-

6
2. LITERATURE REVIEW

The dealing with cooling problems for electronic devices can be dated back to
60 years ago. Bergles (1986) [2] gave a historical overview of electronics thermal
control. In his paper, he pointed out that the electronic cooling problem, starting from
1942 when the vacuum tube cooling was concerned to the 80’s when the microelec-
tronic revolution is advancing rapidly, has been becoming more and more desirable.
To meet the desire, a lot of research was conducted and important papers were pub-
lished.

As described by Steinberg (1980) in his book, [3], in an electronic device,
there are all kinds of electronic elements, and heat transfer could be modeled in differ-
ent ways. For example, there are chip, package, printed wiring boards (PWB) and
system levels. Different geometries, different power levels and different working en-
vironments make the electronic cooling related natural convection very comprehensive.

Elenbaas (1942) [4] proposed a simple model of heat transfer and made some

measurements of natural convective flow between two isothermal vertical plates. His

experiments were conducted for a wide range of Rayleigh number, 0.2<Ra<105, and
a small gap width. Sparrow (1983) [5] conducted an experimental study on the en-
hancement of heat transfer due to pressure drop in electronic device, and provided a
flow visualization method. Bar-Cohen and Rohsenow (1984) [6] developed an inte-
gral formulation for fully developed laminar flow for symmetric and asymmetric ther-
mal situations. Park and Bergles (1987) [7] conducted some experiments simulating
the natural convection behavio of the microelectronic chips. As a result of their study,
the dependence of the heat transfer performance on heater size was gained. Wang
and Ban (1988) [8] analytically studied the low Rayleigh number thermal convection

in a Newtonian fluid confined between two horizontal, circular cylinders. Their inter-

7
est was in the thermal effect in solar collector receivers, compressed gas insulating

high-voltage electric transmission cables, and so on. An experimental study by Tori-
koshi, Kawazoe and Kurihara (1988) [9] focused on the heat transfer characteristics
of arrays of block type elements representing electronic components and the patterns
of air flow adjacent to the upper surface of the blocks developed along one wall of a
plat rectangular duck. Their study indicates that for the fully populated arrays of
blocks of uniform height, the per-block Nusselt number decreases monotonically with
increasing stream-wise distance. Arco, Bontoux, Sani, Hardin, Exterrnet, and
Chikhaoui (1988) [10] used the finite difference technique to simulate three-dimen-
sional buoyancy driven flows in vertical cylinders. They pointed but that for some par-
ticular values of the aspect ratio the problem can admit different spacial flow configu-
ration. They also studied the effect of the oscillatory movement in the boundary wall
on the natural convection flow pattern. But the details of the heat transfer were not
intensively studied Ramanathan and Kumar (1988) [11] conducted a numerical
study on natural convection flows between two vertical parallel plates within a large
enclosure. Their interest was in the parametric study of variation of Prandtl numbers
and channel aspect ratios. Bergman and Petri (1988) [12] obtained the numerical pre-
dictions that describe the potential advantage of using xenon-helium mixtures in the
natural convection cooling of discrete heated elements in an enclosure. This was ori-
ented by the electronic cooling problem in enclosure. Their paper indicates that with
the natural convection the xenon-helium, the electronics operating temperature can be
decreased.

As mentioned before, natural convection caused by a vertical hot (or cold)
plate, has been studied quite extensively. People worked on correlation between the

Nusselt number and Rayleigh number for a long time and have developed valuable

8
methods to analyze the problems. One of the well developed methods is the similari-

ty solution. Sparrow, Quack and Boemer [13], Sparrow and Gregg [14], Yu [15] and
Kao [16] gave the results of their investigations and showed the application of the
similarity solution series. Another method is experimental approach. Due to the ad-
vantages like less interruption and fast time response, the optical techniques have
been employed quite a lot in natural convection studies. Wirte and Stutzman (1982)
[17], Hamady (1987) [18], and O’Meara, and Poulikakos (1987) [19] presented dif-
ferent experimental methods dealing with natural convection problems. As computers
increase their calculation capability, numerical methods are used more and more in
solving natural convection problems. Since 1972, Lloyd and Yang [20], Doria [21],
Lloyd, Yang and Liu [22], Yang [23], and Yang and Yang [24] developed and im-
proved a finite difference code called UNDSAFE. UNDSAFE uses a pressure correc-
tion scheme which was initially introduced by Patankar and Spalding (1972) [25].
Yang and Lloyd’s numerical code can take care of transient natural convection in en-
closures, with or without heat sources, with or without radiation. It can even include
the effect of rotation of the enclosure. The significance of this numerical treatment is
. that it uses real variables other than some transformation variables as quite often
used in numerical convective heat transfer studies. This numerical approach provides
the major- tool to investigate the time-varying thermal boundary condition driven natu-
ral convection problem.

The most up to date papers on natural convecrion not only use and improve the
conventional theories and methods, but also introduce in a lot practical natural convec-
tion problems. Bhavnani and Bergles (1988) [26] conducted an experimental study of
laminar natural convection heat transfer from wavy surfaces. One of their results is

that the heat transfer from a wavy surface, compared to a plane of equal projected ar-

9
ea, increases with increasing amplitude-to-wavelength ratio. Lee and Yovanovich

(1988) [27] present their boundary layer type approximation method to investigate a
two-dimensional natural convection heat transfer from a vertical plate with a family of
non-similar surface heat flux variations problem. Their method is proved that it can
have results agreeing with the numerical results very well. Hwalek and Iyengar
(1988) [28] focused on natural convection mass transfer from a vertical plate at high
mass transfer rates. Their work was done experimentally.

When the transient thermal boundary value problem is of interest, two papers
should be paid special attention. One was written by Kelleher and Yang (1968)
[29]. In their paper, a linearization theory was utilized to study the heat transfer re-
sponse of a laminar free-convection boundary layer along a vertical hot plate with sur-
face-temperature oscillation. The oscillation had an average surface temperature
variation which was a power function of the distance from the leading edge. However,
high oscillation frequency or large amplitude may not be tolerated by the laminar
boundary constraint. Another paper was written by Shaw, Chen and Cleaver (1988)
[30]. In that paper, the effects of thermal sources on natural convection in an enclo-
sure were studied numerically. The vorticity and stream function were used to elimi-
nate the coupling between the momentum equation and the energy equation. This
falls into the transformation variable method category.

Based on the review above, it is clear that the natural convection associated
with a time-variation boundary temperature problem, which is likely to be encoun-
tered in electronic cooling problems, has not been solved so far. The time-variation
boundary temperature driven natural convection is much more complicated than that
with a stationary boundary thermal condition, and it is also different from that with a

physically moving boundary. The time-variation will definitely change both the heat

10
transfer and flow patterns. The significance in the time varying boundary temperature

is that the heat transfer and flow patterns can be manipulated by the variation in
boundary temperature without a physical movement in the boundary, and by doing

this the total heat transfer may be increased dramatically.

l 1
3. MATHEMATICAL FORMULATION

As shown in Figure 2, the configuration considered here is a pair of vertical
parallel plates, each of which has its own time-dependent surface temperature,

(Twl(t) and Tw2(t)). The two plates are dimensioned by height H and are separated

from each other by a variable gap distance D. The two plates are assumed} to be suffi-
ciently deep normal to the plane of Figure 2, so that the two-dimensional-flow as-
sumption is applicable. It is also assumed that the range of temperature variation is
moderate so that the Bousinesq assumption is valid. The governing equations in-
clude the conservation laws of mass, momentums (both in x and y directions) and en-
ergy, and the phase equation.

3.1 Governing equations

Based on the assumptions mentioned above, the governing equations are as

the following:
{ii-die (1)
%:.+u_g§.+vg;l=——g§.p Wit-+33- <2)
aa-Ir-S-Z+vi;§:v<33+‘3%’
.. 3::
p=pRT (5)

They can be expressed in non-dimensional form by using the following non-dimen-
sional groups:

_ 3

Ra=GrPr (7)

 

=3. (8)
H

12

 

 

 

=1 9
Y H ()
$12. (10)
H2
of“ (11)
V
v. ”H (12)
V
T—T
°‘ 13
9=T,,,-T.. ( )
where
:_1_ 2?.
B— 13 ar (14)
The resulting equations are:
aU av
_ V 15
ax ' BY=O . ( )
au an w 3211 azu
_ _ _= — 16
81+U8X+VBY W axZW) ( )
av av av 32v 82v
__= —— - 17
Wurrvay ”3x23? ‘ ’
as as as 329 329
_ __ _= 18
as” ex” av W ‘ ’
The boundary conditions are:
x=o 0<Y<ASP° e=o 9H=o 3—"'=0 (19)
’ ' ’ 8X ’ ax
_. .32 .89. 9x
x..1, O<Y<ASP. ax =0, ax =0, ax=0 (20)
Y=0, 0<X<1: e=ew1(x,z), U=0, v=o (21)
Y=ASP, 0<X<1: 9=0w2(X,1:), U=O, V=0 (22)

where the aspect ratio, ASP, is defined as the following:

ASP=DIH (23)

13
3.2 Numerical solution

Many of the existing numerical methods for solving the equations in the natural
convection problems, use the stream function and vorticity variables instead of the
physical variables of velocities and pressure. While the elimination of the pressure
coupling between the momentum and energy balance equations is accomplished by
using the stream-function-vorticity approach, the time varying boundary conditions
are hard to treat. The unsteady term in the mass balance equation may be significant,
and therefore the introduction of the vorticity transport equations makes the approach
even more complicated than that using physical variables.

The UNDSAFE numerical code using physical variables was developed about
1974 [16]. This code treats the full elliptical equations and handles a wide variety of
boundary conditions. It has been used successfully to solve quite a few buoyancy-
driven convection problems such as room fire [17], hot-walled enclosure with or with-
out radiation [18] and [17].

3.2.1 Cell structure
In order to use the UNDSAFE finite difference method, the governing equa-

tions have to be written into integral form involving a fixed control volume V with sur-

face S.
IpujanS=0 (24)
(Lpuldvl'P-£ulujnde-£(P-PJnldS-L(P—Pelgldv+£°ljnjds (25)
( [pcprdV).,=- 1pcprujnjds-[kr.jnjds+L(uPlu-puj pdv (26)

To solve the governing equations numerically, discretization of those equations
is necessary. The concerned domain is covered by a grid system of rectangular cells.

The grid system is arranged so that the cell boundaries coincide with physical bound-

14

A

 

 

V

 

 

 

 

Axi-l-l .

A
‘

r

 

 

 

 

 

 

 

 

 

 

 

’

 

 

 

 

 

 

 

ny

 

Figure 3 Basic and staggered grids

 

 

15
ary cells. This grid system is different from the "standar " grid system in the way of

discretizing the original differential equations. The "standard" grid system uses Tay-
lor series expansion to approximate the differential equations by finite difference
equations, and all properties are based on the grid points, with step change between
the points. The control volume approach uses the integral form shown in equations
(24, 25, and 26), and all properties are at the center of the cell and represent the over-
all average value. This guarantees that conservation is always satisfied over any
group of cells and in the whole calculation domain. Since the mass flux terms are
evaluated on the boundaries of each basic cell, the velocity x and y components are
needed on the same boundaries. This requires staggered grids. When the x momen-
tum is of concern, the cell centered at the west face of the basic cell is used. When
the y momentum is of concern, the cell centered at the north face is used. The stag-

gered grids for both u and v are shown in Figure 3.

3.2.2 QUICK scheme

When a set of finite difference equations which are developed using the control
volume approach are treated, the main concern is to estimate accurate values of the
dependent variables at the surfaces of the control volume with stable properties.
QUICK (Quadratic Upstream Interpolation Convective Kinematics), first developed
by Leonard (1979), combines the relatively high accuracy of central difference scheme
with the stability of the upstream scheme. This combination is achieved by using a
parabolic polynomial interpolation to fit the control volume surface value at consecu-
tive nodal positions, two nodes located on either side of the surface and the third one
on the next node in the upstream direction. In two dimensional problem, the general

quadratic function of T(X,Y) can be finally expressed as:

l6
T(X,Y)=C1+C2X+C3X2+C4Y+C5Y2+C6XY (27)
As can be seen in equation (27), a six-node grid is needed. Once the six-node

grid information is plugged into the equation, the coefficients Ci (i=1,2,3,4,5, and 6)

can be solved for that specific surface.
The whole domain can be treated in turn.
3.2.3 Pressure correction scheme
One of the characters of the natural convection problem is that the momentum

equation is coupled with energy equation by the pressure. In other words, among the

five variables, u, v, p, P and T, P depends on T through the state equation. According
to the analysis of [22], the dependence is very weak. If the state equation is used to
solve for the pressure and mass conservation equation for density, the pressure cor-
rection procedure will fail. To avoid this problem, Doria (1974) [20] presented anoth-
er procedure, which was originated by Patankar and Spalding (1972) [20]. A brief de-
scription of the idea is as following. At each time level a pressure field p* is guessed
and the velocity components U* and V* are calculated based on the guessed pressure
field. The mass conservation usually cannot be satisfied because of the incorrect
guessed pressure field. A correction of the guessed pressure by means of

P’=P-P* (28)
is therefore needed to correct the U* and V* fields as given in the following:

U’=U-U* (29)

V’=V-V* (30)
where the superscript prime is for the correction while the quantities without the su-

perscript represent the corrected values who satisfy the conservation of mass. The

corrections U’ and V’ are obtained from the residual mass. An assumption that the

 

(3..

ll

 

 

 

1:1

 

 

1

 

 

Specify
Initial Fields

 

 

17

 

A

V
t=t+dt

l

 

 

 

 

 

 

Energy Equation
(call CALT)

 

i

 

 

State Equation

 

 

l

 

 

Update Transport
Properties (if needed)

 

 

 

V

I Iter=0

A

 

 

 

 

 

‘

v
[ Iter=Iter+1 I

l

 

 

Momentum Equations
(call CALU, CALV)

 

I

 

Pressure Correction
(U’. V'. P’)

 

 

l

 

P=P’+P*
U=U’+U*
V=V’+V*

 

 

No

 
     
 

 

 

 

 

 

5Nu/N u
<0.0000 1

  

 

Residual Mass

Iter>Itrnax ?

 

 

 

 

Figure 4 Pressure correction flowchart

 

18
corrected mass flux through a surface of the basic cell is proportional to the gradient of

the pressure correction across the surface is employed. An iteration loop is designed
so that at the end of the iteration, the mass is conserved in each cell. The next time
level may start from the just corrected pressure field to calculate the velocity and tem-
perature fields. Figure 4 illustrates the iteration idea.
3.2.4 Boundary condition

The numerical representation of the boundary conditions described by equa-
tions (19, 20, 21, and 22) affects the total compatibility between the finite difference
equations and the differential equations. The two vertical walls have prescribed tem-
peratures and non-slip velocities. These are a lot easier to treat than the free bound-
aries. The top opening of the geometry has zero derivatives of u, v, T and P. The bot-
tom opening has the leading edge effect. If the entrance boundary is not well treated,
the effect of the false boundary conditions may be propagated into the downstream
calculation. In this study, the real entrance boundary conditions, both velocity and
temperature, are simulated by introducing a set of imaginary walls as sketched in Fig-
ure 5. As can be seen, the surrounding medium will flow into the calculation domain
more naturally. In other words, the entrance inaccuracy of the finite difference tech-
nique is softened by the added reservoir and the flow and temperature at the real en-
trance can be more realistically determined. The boundary conditions at the two verti-
cal walls will remain the same. The boundary conditions at the imaginary walls will
have zero-gradient characters. The extra width and height are determined by the cri-
teria given by [19 ], i.e., D’=4/3D, H’=2/3H.

To simulate the natural convection heat transfer with an oscillatory boundary
temperature, the two vertical surfaces have the following prescribed temperatures.

One of the surfaces has a uniform temperature which is the same as the environment.

19

 

Free boundary \

 

Calculation domain —\
lg

\_ R
Real Boundary, Twl _/ cal 30‘1“di . Tw2

x1

\ I I

 

 

 

 

.<v

 

 

 

ant/I

// /' l l \\\

 

 

 

a—V-o
Imaginary BoundariesT , an-

 

Figure 5 Calculation domain

20
The other surface has an oscillatory temperature as described by equations (28) and

 

(29):
9w=6avc+A98in(t) (28)
where
e _ TW - Too
w" (29)
ave Too

The average of the oscillation, 9 is above the environment temperature. The mag-

ave’

nitude of the variation of the surface temperature, A6, is less than the average value.
This will ensure that the surface temperature is always higher than the environment,
which is normally true in the reality. The oscillation frequency is lumped into 1: by

non-dimensionization addressed in section 4.1.

21

ll

535:3 2:5. 22anth coat-5 e 8:»...—

p 68:. Ecemmcofifiéoz

 

o>a

 

 

APvcmm®<+0>a®H3®

 

 

M9 ‘armerodulol soaring [euogsuomrp-uoN

22

4. RESULTS AND DISCUSSIONS
4.1 Validation test

The case of natural convection from an isothermal vertical plate in an infinite
medium has been studied very thoroughly. In the current study, this case was simu-
lated by the numerical program and compared against the book values. By doing this,
we can test the validity of the grid system, the pressure correction scheme, the
boundary model and the finite difference code.

Usually, when calculations are carried out, accuracy and efficiency are compro-
mised. So before a series of calculations were carried out, the necessary grid density
to ensure the accuracy must be found out. This was as a preliminary work to prevent
misleading by the lack of accuracy when comparisons were made. According to [16],
uniformly spaced 20x20 cell grid is too crude, 80x80 cell grid slows down the calcula-
tion speed and the 40x40 cell grid gives the accurate results not significantly different
from those of 80x80 cell grid. So all the calculations were carried out with 40x40 cell
grid.

The test case is a simplified problem of the real model shown in Figure 2. The
semi-infinite single plate is numerically realized by placing the two plates far away
from each other. One plate has a fixed surface temperature higher than the surround-
ing temperature, and the other has a temperature the same as that of the surround-
ing. Rayleigh number was varied in the range from 1,000 to 747,400. This range of
Rayleigh number started from 1,000 is because cases with lower Rayleigh numbers
will become conduction. And the Rayleigh numbers are no higher than 1000,000 so
that the calculations are done for laminar flows (the number 747,400 resulted form the
characteristic length, it is close to 1000,000). The heat transfer coefficient, Nusselt

number, results from the calculations is compared with the corresponding theoretical

23
values for laminar flow in Table 1. The Nusselt numbers of the calculations are within

4% difference from the theoretical values, except case 04, of which Rayleigh number is
747,400. This may be due to the transition from laminar to turbulent flow. These er-
rors were sufficiently small that the numerical code including the uniform grid was

thought to be suitable for treating the posted problem in section 3.

TABLE 1 Validation test

 

 

Test No. Ra Nutheoretical Num Error (%)
01 (20X20) 7474 5. 3 5.92 1 1.6%
02(20X20) 74740 7.94 7 .10 10.5%
03(40X40) 1000 3.9 3.98 2.1%
04(40X40) 7474 5.3 5.51 3.9%

05 (40x40) 74740 7.94 8.15 2.6%
06(40X40) 747400 15 .49 16.8 8.4%
07(80X80) 7474 5.3 5.35 0.94%

Figure 7 plots the error vursus non-dimensional grid size.

4.2 Discretization of time step

The oscillation of the temperature of the boundary surface is of primary signifi-
cance in the current study. For transient problems, the accuracy and stability are all
very sensitive to the size of the time step. How to choose the calculation marching
time step which takes care of accuracy and stability is another preliminary problem
which needed to be solved.

One non-dimensional parameter describing the stability is mesh Reynolds num-

ber, Re, as defined in the following:

Re=[(AX)2+(AY)2]/At (30)

 

 

24

Amuo< .Nzoofinm Jun—m3
.3958." :e .809 on? Saw 5 ensur—

 

Nag
0n00.0 ¢N0.0.0 0 0.0.0 N £00 00m0.0 0000.0
000 — flex 0 < . 0
. {Knox o
. ovhthuom b N ..
. 00vsvufla¢ n
s 4 a
.. in
a N s W
l e 10 —
. r
T q I
fl s
p p p n e P p p b

 

 

 

m—

25
The stability criteria is that the quantity, Re, keeps in a range when AX, AY or A1:

changes. As the grid gets denser and denser, the time step should correspondingly

decrease. Table 2 gives the results of a series test calculations. In these test cases,
the grid sizes, AX and AY, are first fixed while the time step size, At, varies within a

range. Then the A1: is fixed while AX and AY vary a little bit. The quantity
[(AX)2+(AY)2]/At was calculated for every test case. As can be seen in Table 2,

the numerical calculation is stable only when the [(AX)2+(AY)2]/At is within a cer-
tain range. This verifies the stability condition mentioned above.

TABLE 2 Stability test

 

 

 

2 2
Test No. AX AY At ( AX + AY ) 2811:1103;

A1:
01* 0.025 0.025 0.006975 0.179 S
02 0.025 0.025 0.06975 0.0179 0
03 0.025 0.025 0.3488 0.358 S
04 0.025 0.025 0.03488 0.0089 S
05 0.025 0.025 0.0006975 1.79 S
06 0.025 0.025 0.0003488 3.58 O
07 0.025 0.0167 0.006975 0.129 S
08 0.0167 0.0167 0.06975 0.0796 S
09 0.025 0.0125 0.006975 0.1 12 S
10 0.0125 0.0125 0.006975 0.0895 S

* CASE 01 is taken as a base case and every other case has some factors different

from those of the base.
Since the grid is 40x40 uniformly spaced, the AX and AY are fixed for the spe-

cific problem. So the range of At is limited too. On the other hand, the goal of the cur-
rent study is to calculate the oscillating surface temperature—driven temperature and
velocity fields. This implies that the time step should also be small enough, so that
the oscillation character is accurately described. This can be better illustrated by Fig-

ure 8. In Figure 8, there are three curves. Each of them represents the overall aver-

 

 

26

GH0< £20070 Jun—mac
55:58 .5 “cute on? 3% 2:2. w 2:»:—

9 6:5 Ecoficoéwéoz

 

ON 0 F 0

. _ om
.. a to...
. , .. .
r . low

I
I I .I
1 10m
.. a ..
,I

T , woo_
. H
I acoEmomlm o. .e IONF
LeoEmomlop $10 ..
acoEoomlow I

. _

 

 

 

(1) nN ‘roqulnu llosan

27
age heat transfer history of the same driving force, but each has a different time step

size. The smoothest curve resulted from the calculation with the time step 1/20 of the
oscillation characteristic time, the period. The curve that lost its smoothness at the
peeks and the valleys resulted from the calculation with time step 1/5 of the character-
istic time. The calculation with time step equal to 1/10 of the characteristic time has
improved significantly in terms of the carrying the complete information of the driving
force. With this time step the amount of calculation is reduced by half, but the infor-
mation lost is very little compared with the result of calculation with 1/20 time step
size. For the current study, each oscillation period is divided into ten step to carry

out time marching procedure.

4.3 Mechanism of natural convection heat transfer with an oscillating boundary
temperature

The character of natural convection is that the temperature difference initiates a
fluid flow and the flow can have all different patterns as a response to different ther-
mal boundary conditions. In return, the flow affects temperature field distribution
The thermal boundary condition not only arouses the fluid flow, but also can manipu-
late the fluid flow patterns. As a very expressive example, the natural convection
driven by a time varying boundary temperature shows how the thermal boundary con-
dition variation leads the fluid flow movement

The numerical simulation was canied out with the boundary conditions de—
scribed by equations (28) and (29), and Figure 6 in section 4.2.4. It can be imagined
that as the boundary temperature oscillates the temperature field adjacent to the hot
surface will responde to it in a periodic manner. But what is of the interest here is

that the fluid flow pattern changes periodically too.

28
The group of pictures in Figure 9 show the temperature field and velocity field

time history during time of two surface temperature oscillation periods. This group of
pictures are resulted from calculations for Ra=7,474, ASP=1, 1:0.001, AX=l/40,

AY=1/4O and A'c=1/10. The number on each figure means the time step in sequence,
and every ten steps make an oscillation period. It is found out that for the specific
case, starting from a quiet initial condition, it takes about 500 steps to reach a fully
developed time varying response. For the purpose of relating the temperature and ve-
locity fields to the corresponding surface temperature, a detailed set of figures in the
time range 560~600 is attached to each individual instant field pictures (temperature
and velocity) in Figure 9.

First of all, it is observed that the isotherms are not always open curves
starting from the bottom entrance and ending at the top opening, as is seen for the
constant temperature vertical plate. Instead, some loops and some heavily bent iso-
therms occur. The thermal boundary layer analysis is no longer suitable for this type
of natural convection. It is also noted that the existence and locations of the new pat-
tern of isotherms vary in time. Actually, the isotherm pattern variation has a periodic
behavior.

Secondly, the velocity field is seen to be more exited than that of a single plate
with a constant surface temperature. The fluid can flow upwards along the hot surface
or form a circulation, depending on the corresponding surface temperature. Again, like
the temperature field, the flow field changes its pattern (either parallel upwards flow
or the circulation flow) periodically, too. The transition from parallel up flow to circula-
tion actually represents a different mode during the circulation core developing histo-
ry. As can be seen from pictures in Figure 9, the circulation is always there. As time

passes- step by step, the core will move the bottom to the top, then out of the inner

29

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Elli“ FRI! 11mm ID

'0
2
\\K
\\ ~
\\
\
\\\ \

III'lllrttlllltllllltllttllTFFtTTTr
l a
E m
L1. 1 .1
. 0"
I -3,
' -2:
I _§
I n
! -d
' ‘II
l 4:.
p;
d:
C.
t§§35 ._
)—. 9:), .1?
1 fig
5 -
l -6
l
l _&
|
l dd
I 45
b
: ‘5
i -§
; du-
' 5
| .‘U
0
.‘f“ .13
l .l
l l .4
l 1' _
| -
l l -
_
d
1

 

 

 

 

. / /
’ .
41/ l l l l l l l l l 1 L1 L I J l l l J l L l l l l l l 1 l l l l [J l
‘ . ‘ . ' - :-‘°::-::o:oo?o.-..-.:..no.. ’..'..--—'o.- ::.:’.:o:: :0 . sci ' " o "
’ t: 3". ttttt ‘6 0 oooooooooooooooooo “.”::..o oooooooooooooooo 00. r-.." '.-I
It) i : l '. :
l .' i '. '. '1 .n
‘ ' O 3 I
-> i ‘ ' T
: . g
u 1
I a 9
: . s . . -'
L. I. I t.
1 0 -
I ' n
"‘ r ’I‘ . : é
_ ,z t
I- a E
l I
III I 'I
.- ' ’o’
I a” :
.-“...O..l'.-"‘ '
~ I,
I
I- I
I-t
‘
h ...........................
)-

.0

1' 1m CI!” [Um ' mm

 

 

 

 

0
4.1“

v

v

CHI-ll m

IzgnglgsllLlllll

p

 

 

 

..._. __ —— m
__ _-., ~:‘—-“
fir . ,D( fi..._. ~_ -‘

Figure 9 (1) Temperature field (a) and flow field (b) responses to
oscillating surface temperature at 1:560

(a)

 

(b)

 

30

 

 

 

 

 

 

 

 

 

 

 

 

 

I.‘£OO

 

 

 

 

 

 

 

 

(OHM FRI! timed) I'D

ll] 4/ /

‘- IILLIIIAIIIJOIllllJOllllLLllllllLll

 

' l I l r I r r r I r r . r r l l l l r r r I r I r I I r r r T T I r
Hr
.. 3
r j
)1» $
.. .1 g
a. .1 a
hot -4 O
F _ II
to -‘ i
~> '1 if
i§| §| § § 3 5
.. . .
m - a
.. _.
i) d d
I - lg
n. d 1
up | I — E
-- l j l - 2
I i i : g (a)
rt» ; '1 5
r g‘ 8 8 | l [r "
-t> 6' NI I I ‘1
. I l l ..
—(
_.
d
3
_(
..
_.
.1
'1
.4
j
l

 

 

 

 

.. ~ .-
. -o- ' —'

----- . ----. .e
._ .n...-.. 0 - _ ..---'

 

r- Infifi I"l--3--{.I.':.':J..;"{.j;,;v_11fi::{-.;..—J~-J""I J 3"T"F In? I l'Iv‘l l l l

v
—1
“fl

1 1 1 1 1 r I
0.23iltflt

71”.}! 3

l

 

p.
- ,a
‘\.--a-

 

 

 

(1))

1.1100

-"
-a'
.«4'
,.
.......
..................

-l.5

 

-l-J-iole-L—LJ l- 1.1. l.J_.l_l_l_fll_fil_l_l__L_l_L.1_l_1_l__Ld.fil_J_L_l_l_Ll
can“ warm or 0.3mm

. HLIEIZILLIJLI IILIC-I:T‘I"T_I—I—‘I—I I I I

 

Lmlflfl FBI!

.1 ._l. .1. A

     

fl
___:_“__:_:::-—=I"m=—-~._=--_

Figure 9 (2) Temperature field (a) and flow field (b) responses to
oscillating surface temperature at 1:563

 

 

 

lTTllTITirTTITITTIT.T‘T. .TTYYTT I

ILSI ”1(4)?

 

 

 

 

 

 

. 100
200

0.30““! 4” PH 1. )I :

 

 

(a)

 

 

 

 

 

 

 

 

 

 

 

‘mlflfl "NEW“ 0‘

 

100
—. 200
l .9500

 

 

 

 

 

 

 

 

000(le “'13 lo

I 641063 1P?”

1.1 l_LJ_l._LL.l_J._1__l_.l_L_L_l 1‘1 1 l l l l l l l l l l l L_1_L.1_L_L_L__L.l_l

 

' i ' ‘ l I ' I l l I ' ' t ' l l | 1 l l 3 ‘

 

.....
.....
............ “g“-....”_----oIooc-‘-‘--‘.- .O"

1T1. T‘L.) I l 1"F‘1'-¥°-F--F-+-(--¥--§"P‘T ‘l I I L L--I'"Y T T }—'T i T T [E

.-
“'..'.o... .0--.‘--‘--
.5 v.
a- ~
.o ‘-
~~

”l3.” 3 -0.l||.7l-OI

can” “HEM ‘ 1m

 

(b)

[0 1.10m

u ,
- o
..........
..........

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

CON" man om

 

 

 

 

“L

. __ 0
_-_-_;:T:.‘—*‘-'- ms——__

 

 

Figure 9 (3) Temperature field (a) and flow field (b) responses to
oscillating surface temperature at 1:565

32

a

 

 

\II
)
(am I b
I (
23.2.53 2 ”:5; .2, 3:3... 3 (it... 525... 8B... 3 3.5.5:... :3: :26. 5.833...- "3.2... and “I 123.... Giza”. RE; 3
ID.D
_ 4-44.4I4I; 4I4JIJiqla|qifiJIqI«JIJJI714IJIaIaIJl144I4IH .4 4.4 #4-. T. a d a u 4N 4 4 a a #4 _ 11
I I o.
I. T
I. T
[I Y;

 

 

 

 

 

 

 

 

 

 

1+
1 i
|
vI 1
I I I I l
I I oomtIIIIIIII [III Ioo.I.II I I I I ,
II .IwIIIH-IIIIHIIIIIMIwIIH- I- I , I, ,. / I
I.-I - .IIIIII.:--;I.IIIIIIII-IIIIIIIIIIMIIIII III I S» ,
II I I I ISTIII IIIlIIIIIIII1o2..- -HwIIIIIIH.I.II-....m. .../ . / I
III- IIIIIIIIIIIIIIIIEIIIlIzIIIIHHWXIH.U/, f
I . I-- .I.;ImHI:IIIII:IIIIInuuuummuNma,w /
95 8QIIIIIIIIIIII... , avg»
IIIIIIIHNSQ //9!
E... I .13, I! I/
mLfimWNHMWWWWMWEMVL,

 

IIIIIIIIIIIIIIII

,.-a.oo------~~s...,‘
.0

 

 

. _ _ d _ q _ . . quofijququ4J--jlql

  

3.0. g 528.4

“lllllll

 

 

 

 

 

 

HI
I H “I I .l I I I N I
W F F I I i

  

field (a) and flow field (b) responses to
face temperature at 1::567

'llating sur

OSCI

Figure 9 (4) Temperature

33

" TTI IIYV*VTIfirY r', vyyy, rvfir

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. 400
200

 

 

 

 

 

 

 

 

 

 

 

 

.- l l .1_J_J_l__L_1_L,l_L.l_J..1 .l_.l .l. J._.|._L_Ll l l l l l l l l l J_.1_l__L..L_1__J.| .I

II

II //./////

III I/ ,I/ =

. I / [x /

~ /§//

w ////

74. JllllJl'llLlllllllllllI'lIl'lL

 

_:- a“ ———_._...——- __.—'——' ———;—-" -" " r

’ 2","- ——
it _' OT

 

 

 

U.“

o...
................
.....
s.“
‘y

x 1 1 x

.n-.
.....

l l 1 r I I t J
I I I 1 L 1 1 x
.U

I
l

......

I I I I I I
I I I I I l I

l

 

 

11111lllllénlljlllllLllJnll

 

 

- . '
.go.-....-...o-..-Q.OOOOOOO 0‘ I ‘

..:.--~':7';«‘:.-:~ "~l--L._l I I' I’ll I I l l L l L L ['4 I I l

 

:
o
a
l
I
I
I

.‘.6|1‘,9[ U9

0.9M”)! -UI P I l i. ll 3

lm'lilfl "HEW”! OF

I , 'BHU

Immm mun I).ll0(fl(H)) IO

PHSJI = -o.sozq5£-0|

(III-R III'EIWI. ‘ 0.21111)

3.1000

V.

CON“ M -0.6flm

(a)

 

(b)

Figure 9(5) Temperature field (a) and flow field (b) responses to

oscillating surface temperature at 1:570

34

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I‘Ollalfl FRDH 0.110((DE ND {0 [11500

I [I IT T 7 l j ' I

I. I I “I

”P ' I "I

~- I a I 4
~- = : I -I 8
-- II ; , ~' 8:
I... : 5 f 5 ' _I g
”’ I ' ' :‘
In}. ! i I } ll
"' I I , 2
"" I | I _I E

P I I I
:_ % §I § I 8 I I I ‘1 E
I. I I I I 3 g
.. I I If I a
HI- I I 4 .1 3
~ I I i I I g
I ‘ I w
l I I -I E
I I -1' a:
I a (a)

I q .-
I . E

I L -+

I I ~

I (I _

I -1

q

J.

.1‘

I

 

"JllLILILLI'IIIJI"Ililllllll

 

 

 

"' “ IIIIIIIIIIITIIIIIIII

l

I I

JAINH-‘fl

 

 

 

cm Wt I tmi-Ol "ll-3|:

 

(b)

 

1’. La.

 

 

 

 

 

llllLlll‘nllllllllJlllllllllllLll

cum "I! '0. sun

 

II; I‘
I! ‘a‘. . Q‘s».

Lin--. . I J -. {'Isi’ 1‘1.l-l l‘ I I l.l~l-l~ 1.1.3-1-‘I‘IZI:‘I‘~'1:1‘I‘I'Ii
Figure 9 (6) Temperature field (a) and flow field (b) responses to
- oscillating surface temperature at 1:573

 

35

(a)
ll

9?. 9.3.0:: u .m .«..L 2. gm c 59 $232. c.5352 8m...— 2 pl: fl-Uwccé 20K... «HS—’5.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1411-74444 _ _ _ q a A q q q a _ 14I_IqJI4|dI-14Iq.J-4J4.41. I. 4.. .._-

Y I I I :3 . III. I- c?- .. IIw--.-;-.. I

w! I II I III I I III III; -II. - I..- I...-.-..-.

l IIIIIIII IIIIIII II III-HUM . - -

- I- II 8,.IIIIwI I-@.--IHWI-II . .

I HI. [III/[HI

I :8 I Q;-
llf’l‘lll I

- a3. ISn-IIHIIIIII/h.

ILLIFLLIPFIPLILIPLLIEI r1 _-_.. ._ -_.I_ LICIQ. u

;..-.3- u : A:L Shun-C .5 .3653— ‘53.:

(b)

:an... a =an :- EE EQI‘:

 

.fiqfid—quqdqdu_-—qquqq

;.za' ' ' '
/
I
// /
/
/ /

__—_ _ F- ‘ ,
_ _ - ’ fl»

1%. .43ajJJJIfi-q ......_

 

 

 

 

 

 

Figure 9 (7) Temperature field (a) and flow field (b) responses to

oscillating surface temperature at 1:575

36

(a)

a? u: 3%.: u 2 .27; 2v ETA. & EB... :29“.

 

 

 

 

jal— q H — q q — u fi‘ « q dId q q d d fi q q J_.-u_ _AI_ _ _ _ _t_ _ _ _I4I.qI_ _u I
T I
T I
J
L
J
I J
I 1
T 1
fl I.
T! II III II IIIIIIIII III | I. II. I
.. I I III|8~,.I l-rI {-83 I I: y I
a... II I III I III I I ll ., I II I / ,, I
r II .I III! I-I:l I I- 11; / 8v, ,. /1

 

 

2E. .

D. -an—HDAA. A. .31. i935.

 

 

 

 

 

 

 

 

 

 

 

(b)

ZYuahud- ”3.3—; 2.6 t 538. :_iu 836

2

3.9. ‘E I‘d—3.4

 

‘-

 

i~~=<=s- :
C)
6‘

’ ’——-

.’

 

 

\2’33’

u—u-_‘d«—__A.J*__q.4_—q«~_q_~_‘q_.qfide

 

 

 

 

 

'Illll‘lllll
Ihullllllulllll|IllIll

 

 

I'1--9.4"1“1. 1 1

 

 

 

‘r.1 1’1"}:1-4-v1 I t r L 1 x t"r~+:T.:

 

 

 

0
t
S
e
s
n
m.
7
$7
[a
WT
(m
we
er
am
mm
flD.
dm
an
e
)
C
@m
dr-
Mm
foe
e.m
ft
.m
fl
C
S
o

(8) Temperatu

Figure 9

37

 

IIITTIYTTTTIIllllllllYlflllTrlllT
, I

l
I

.4
_.
__

0.12278t—U8

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I
.J
-I
‘4 2
g g g :55 I I g
I I I -1
I I I I - II
; I I I .4 9
I I I I I .
I I I I I g
II? I I d
I? I; I I I j 3 (a)
I: ,I I : I ' ‘ 5
P 8 ' 8 I ' ' '1
I8 g T I I N I I :I 8
II I I I I I "I”.
I ' I I I I .JI -
I I I I I I I I :I 2
I I
I- I/III’I /: 1| g
' ;- 5/ ‘6’” ”I I g
I r/ 5-
32...- #1 I I I I ' I I l I I I L LL} I I l I 1 I I 1 1 L4 1 I I I L I II II
III / u “ “—— "———'"‘ --" 7' V 45:
II I “
III §
III a
III ,, 3
III -
III' " 3.
III
IIQ ‘
II
II a

COM” IMEM U

I
I

(b)

3.5m

4 l I I
l I I I l I

1
I

in
II:

I

Figure 9 (9) Temperature field (a) and flow field (b) responses to
oscillating surface temperature at 1:580 .

 

 

 

 

 

CON” "I“ 0.0011110) I'D

lllLllLlLLlllJLllllJlllllLllllll11111

 

_ _ -—“
-"‘:II-l-S-‘ ‘ ._——--._

38

r-

'fTTTT

)
a
(

17:: :7ugn—wé m9 I733:— Kan—Eu

WI 744.4!fi j-I1IfiI—II_-I4|4IJJIJI«JIJI4leqI«IJIqI 4II_IA--.1-4IJI«I_ .IqIfi q q .I

P. 91 flaucoé x9: K629.

V

I

‘LllLLl‘

l

'1111111LL" I

2.0 3 Eu;— uni—Eu

(b)

36

 

 

 

_..,----...

 

0-0

2:1 ‘: :: ::

O

 

 

14q._n.fi .~_q_. . qq _ qqq q _ ‘4 114-4IJJJIJIJJI141

    

~_—_

 

  

__.._

           
          

F

__‘ ”MM

‘71—:

  

. ._

0.. —1 ——.___——_~ ._—_‘

oscillating surface temperature at 1:583

 

 

 

 

 

 

‘I r J. :‘y --L-I:J::! -sB-J-.J-.J..l.-K

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

8a. 8c. //
H fl H HI H 4 - — H H H N II 1‘ fl H 1! rfi fl H H H H INJq — I f I

 

 

 

Figure 9 (10) Temperature field (a) and flow field (b) responses to

39

(a)

 

V

m? gwmmu A. u .m 4 . E 2... gas. 8 dam—z. Caz—.9“. 8mm..— 2 Qt H:~o:.a 155 I625.

 

 

W a «latqzfiifi.fitfilfildlfila§414 _ q a _ q _ _ — 4:4 4 q q d Al414i414|4l4ifilqld!filj
TI :1
F I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

33113-14 .1
5’-.- 'n" .'
1

2'. "95:33;

 

on- ‘

..
Osoq-o.-¢--'

‘~-...

~
‘.
.‘~‘
0" 1
"“

armed Kai-pm

pa:

0“...

---.-s....o..-."
oo-
"“'--~-oooooo---ao“'

‘.
‘~.-.-.

_ a
.--o—-oo-“‘.

3.: 8 $2.: :30...

mu

1

 

0.0

(1))

3a

 

u.—

34-

5E :30»

<q—fidc—1d«<—.du—V————~4qq_qqdqddq_qdqfilq_

 

 

 

 

 

=585

re at“!

field (a) and flow field (b) responses to
face temperatu

'llating sur

OSCI

Figure 9 (11) Temperature

4O

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

                     

   
     

 

) g )
a V b
( “ (
a: 33.3.: n .m 4.... 2. Wavinué 59 gr... ”30:6“. 91%... 2 322-36: 2 l9... Ina—E. two-.0 "5.3:; 36 t guns: i=0”. 36 P— g.” g g—iu
14 d u q d a q _ _ _ ficaeuqth 1~.l_14|414|41fil—14|AI41 HIqIAt 314: _ 1.. J ‘4 _. V_. _ 1.014411141qu 3114113141‘
VI IL 1.
T. IL F \.
VI .4 '1.- \\ _
T L ”do... 0 x ‘
v L on.— .
T l 1... , h
a 1 .00. M
IL 4” N
. .. r
H . .n =
.. .. a .. =
a .. .. ,. . .. a.
1 .. .. . fl .. .-
.. n . u u 2
L .. m m .. u 0. ;..
- 1 n n u a u c :_
T 1 .. . u . u .. f
u . u . ;_
. L .N. .. . u " _,__
l o... N “ “o .
u :m. ... .. _
..." \
fl 1 on.“ o?
l .w w
- 1 .“m
Y L nu."
\il l w 411 if r L n." /
... /
. .-
1 lil ll» ill [’I|!|a II I I no.
[I‘ll I 11 I. .z
I I." 0 /
l I .1111 ll! 11 J: ll n..tx-|: I!" . in ,
In) §~
I - 8m - I «loc~f|.!’- ,. 1! ... \\
T! lllllul I .l lllll!|.l[l|1yhvllr .ll.:.// / // . ,1 l s...- .x.
..,. w} I 1 lili- - I I M Huh _____ Jeff;..-/./,.,.,s.,. . .. \
. . I]!!! j ,/ .r / .u
r. 8.. will lion; lifrlltf’H/Izy/ / L \ \\
I l Illi'ililljipf/IHUI/ // // / O.O\\ 1 1--
.I I!!! I/ // lllllilll '1
o8. Semi/(Woo. / // \\ \xll:
.. , N am.
. .. ..... .|.| Ill! .1.” -Jhufi v; I
innuflflflflflflfl..InnuuulflflufluEEWmmmmuumw lulu... n1. urn ,

 

  

 

 

‘——__——_—brbpww_~_---_r~HF~_~fiH4HA~

 

     

Figure 9 (12) Temperature field (a) and flow field (b) responses to

oscillating surface temperature at 1:587

41

 

 

 

llllljlllJIIIILllllllllllllllllllllllll

 

 

 

 

 

 

 

 

 

 

 

 

. 400

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T.‘IlllfITlTllTTTlTTIITT'

YVIY‘FllYIY

 

__L11L1LlillIlllllllllllllllll[1111111111

 

l lllLllllLllll'lllllillllllll[ll
-.- ‘.'I--L‘.’.L-_4 ,,,,, 1. "1’31”. ,P'T‘I r

 

.1 Al
u.cr

UV

‘V

 

 

    

LlllllllllllzulllllllllllAnlllllllllLl

 

I).3201l[4)8

(OHM "HEM N naimm PHIJF

L9500

(Dual WON 0.0KLIIND I'D

0. ”927

0.8!!!) Pll3.3l=

[0 3mm can” IIIEM ‘

-s.1(m

CON“ Fill!

(a)

 

(b)

 
 

 

 

.4
.1

 

 

 

.m
.62;
.‘IOO

—200

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_ 600
.10)
WE 200

+444¢iildllTTlelllITlIlTlllIIllIT
:\¥_

    
 

 

 

 

. 1 , .

l J.l. l l l L l. l 1_l.J_l 1 L_l .l .1_1__L AL_.J.-_L.1_1__L_J._1_.1_J_-.L_J.-.L.1_LVL_l_l 1.1

 

'_‘J.-1"[,.l"l ‘]"T l I

 

 

 

   

.- . - .. .-~..Z:_;;;-_:-_:.::.::. :.'.':.:::..:::.::..----:::::::.I3”" ''''''
I
r- n
I : l' ‘
.- r . . \
F- a .' ' . .“
I : I. \
" ' I : - 9.
. , ‘
I . : N. (I , .
P : 1 : J «J (V. .
' l
P 1 " 2 c '
l ‘ ‘ |
b o . ‘ |
1 ' ‘
I- , " 2
‘ I
— c '. t x
. . , _
p—o : l l
.- '- '. .
I u t
I— 'l K ‘2. 00 ....
’- II ‘-
I h
r— '. ‘.
D— 1 ‘~
l— \ .‘ g] ......
P II \ ...........
h .I
n- “‘
b A. ...............
‘\ ......................
H o. a
t—l ? '
H \ ‘\ ‘
p \ .........................
\ .....................
A \
H \ \
r-q \
H k
k k 0 . G .
“I U ’ \
O") \
H . \
H \
H \\
p-
H
-___ _____ ;“;§5?EEEE===——_

     

11111]!

 

1111410!8_|||||||ll.1_Lel01_J_i_L_1_1_L_1.1—LLGI£.

 
  

 

U. l7lSnf-dl8

“NW!“ IMHM W 0.9)(IIE 4)! Plll.3|=

'J'SOO

'9"le {Fan I) Ilnrlufun la

e.m-0|

Pll3.3l 3

COM” mum U 0.3111]

3.5000

[I

can” FR“ 4.01113

(3)

 

(b)

Figure 9 (14) Temperature field (a) and flow field (b) responses to
oscillating surface temperature at 1:593

43

)
a
(

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(b)

ncuwcnmvdu ”3.3—; 3.9 * gm...— ‘l—‘u SEA 2 95rd- it I...‘...

 

 

 

 

 

 

 

 

 

 

 

a: 3789“ .3 n 3.2.; :7 9&6 Ia gr:— C‘_§u 8‘5— 9— 91539.4. l9”: 39:5.
.774. .JJ 4 _ . _ _ q _ .4 a q . A 1:14-7fi444; ._._ .414 144.44-
- L
l L
l L
- L
W. I ‘l I II- a I: II} I I / m
Illll: 1-: [SNIIIlnlllil IISN f I]: - -/
Will-.| III: [Ill (I . /
Will! _I...--- .- AH--- i-..|--Itl.lJJ¢Hw///.é.,,. ,. /
W!III.-lliul 8... il 8.»le l/iln/x. 0, ., ,. ,,
llllr [It]?! . .v I. /1/
I l//,
o8. 8o.|’ll!lfvlt///ce
.961 949;}
LL _ _ p b _ r _ _ LLLLLLL-F|_LL FLLIrFIFFFFFIWWT»! h

 

 

 

“0

ANN _ q q q _ . m? d a d _ . a .IqJIficfiJJ 7717.11.14. 0.9 A} a

l"

 

“III-'5'Til-l‘4 1.3—7 I,.L-'r I I I

u .. .,.
u .. v
.n . u
m m n
.‘cy. H mm —l \.. .
vs». H Id .-
z .z .
tum m
x ..
T."

u I] 1. k '--r:~1~1:::::.‘::::.{-

 

 

 

_ _ _ _ _ _ _ _ _ _ _ WHLI .. - IHIIHHLIHIJHHLIPL.

  

595

and flow field (b) responses to

llating surface temperature at t

(15) Temperature field (a)
oscn

Figure 9

 

 

 

 

      

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

)
a 7
( ‘
$35.3 "3.2—mag... 3.332.: 52.3 88.. 28.888... 3.: 52.3 $2886-25: a... tux-E... «2.3 884 3 Raw... .2: I269
qdfidqqadfid..4...1qd____.____qfi«q4fiqqq. .q._q_.d __q_.q_w.w_q_d_1q_..d¢fi_.q.....d
r I r
l L \\ - l/ I
r i a ...................
T l ‘\ in I.
l J .\ \
Fl I. .s ...........
r 1 t c ,9. .
.I ll \, J .
r 1 I .. . ...
r J .l ..\ , ... .n
I m t .. u.
I 4 . - ... w .n
..l I. ... N M ... .. m u
b
I l .. .. . c. . u u n .
I 1 .. ... m .. 5 T n u u “t
1 . ... .. o .. 4 .. .. .. ..
l .I .. . . .. u . ‘ ..
1| 1 .o u u o p. N N -- ~
I . . u n . c . .. t.
Y. 1 n u u n . u, . N
.l . . . . . e .. .
I l u u . a; 0 u ‘s
'[l’ T c o a. a ~ .~ ‘
rl (it! I1 .. 0 a a. 1 ..... ~ ~. .~
/ I u a. ’4 . u
//, 1 I ... ... ,. ,
8N.
TI
8v.
80.
”AH:1
Khalil
Lang.
rJlJIHJJITuI a m A . . _ H .

 

 

 

oscillating surface temperature at 1:597

Figure 9 (16) Temperature field (a) and flow field (b) responses to

45

YTlIllTTTllllllllllllllTllllll
\

i
I

 

.4
.4
.4

:1

0.19I‘IZE 438

 

 

 

 

 

(a)

 

 

 

 

 

can” "Him W 0.91015 0| 91”.“:

 

 

 

 

 

 

 

 

 

LEBW

 

 

400
‘“‘. 200 '

 

 

 

IWWIIIITIIIIIlabillllTllll

   
 

HOMER FRBH 0.0GIIIND I'I

 

llllllllllllJllJJlllllllllllillllllllll

 

lJlLLlllllllllllllllJ'l'illllllJ'll ‘

 

 

 

p.00. ......
......... .O

M
M
#
H
-
-
F
l.
H o '
H
y-
-1 O
" 8
~
_ E
‘3."
” _._
L d,
_ ........................... E
.—
~ I
F{ U
H» a ’
H .4
_1
—d
—4
.4
,4 —4
j -
d -
n —1
H —I
H d
u -1
J q
~ —
‘-

.1—4'1uo'1 Lu: 1 T
--"-- 'I"

l v i
x
O
:.
55
as
.‘—:
“l:
N
x:
:z
z:
:

 

 

(1' l l l 1 1 1 I l I
Pll3.3|= 0.5571E-0!

U

a
‘\ 3’
s a
‘0‘ '4'
.--....-o--

 

(b)

3.110!)

a
.v
,-
,r
-"
‘0'
.....
,-
......
.........................

I.

 

 

 

-Z.1!II)

v v

CHI“ rm

  

 

 

 

 

v—nggv v v 1 V

Figure 9 (17) Temperature field (a) and flow field (b) responses to
oscillating surface temperature at 1:600

46
Put the temperature and flow fields together, and the oscillations for both of the

fields seem to be generated by a physically oscillating hot surface. According to the
temperature and flow fields appearance, this pseudo oscillatory movement of the sur-
face should be vertical. This is explained by the fact that the natural convection adja-
cent to a vertical plate due to gravitation is primarily parallel to the plate surface. The
term "thermal turbulence" is introduced here to emphasize the phenomenon of a time
varying boundary temperature generated natural convection flow.

Just like turbulence can increase the convective heat transfer coefficient, the
thermally driven turbulence will dissipate more heat from the high temperature sur-
face. The numerical simulations proved that the enhancement not only exists but is
also very significant It appears that the frequency and the amplitude of the oscilla-
tion are the major factors controlling this phenomenon.

Although the analysis is done for a specific case, the thermal turbulence idea
pertains for the situations with a time oscillating thermal boundary condition.

4.4 Effect of oscillation frequency on the Nusselt number

It is easy to imagine that different surface temperature oscillating frequencies

will cause different characters to the oscillations of both temperature and velocity

fields. Therefore the heat transfer will change accordingly. As expected, a series of
numerical experiments, with Ra=7474, ASP=1, AX=1/40, AY=1/40 and Atr-l/IO,

eavezl’ and A0=0.5 unchanged, and frequency varies from 0 Hz to 1000 Hz, show

that the heat transfer is enhanced by the oscillation of the surface temperature signifi-
cantly. This is visualized in Figure 10. As defined in equation 27, the driving surface
temperature is oscillating around an average temperature. In Figure 10, two things
should be noticed:

1. The space-average over the plate is oscillating about an average value. At

47

Surface Temperature, 9w (1:)

me

o
(O
I
I

0 L0
rt) <1"
I I
I I

LO
I
I

O
I

 

 

 

.552:ch 9:352—

m. :

_ _ r! .

-.J
.J
J
3

fin? 539$ ..umwfi
Eu. oust-a 2: 3 3:38.. .35.:— ..omm=z S o
p .25. EcommcoEEéoz
o. m w n
- .

p F . _ _ _

Sur—

 

 

 

 

 

 

 

 

'3
.1
__-|
.J
.J
J
J

‘

——-——---—-——---—-——-
——-——————-————————-—
-——————————————-———
_——-.
_-.—-—
————-
_-——“
_ —————-—-———-————-——--
_ ———_-———---———-———--
.—
‘—
’-c—‘—
...-—
-"’—
. ——-.—
.._ -—————--———-———-—————
.---—-——————--—¢.-——-
‘-
——’-——
I‘-
-—‘-
-‘——
—‘—‘
————------———————————
-d
-—
——---
——“
‘—-
——‘—
-—
’
_, -————--——-——————-—--——-—
-
‘c-fi
—--‘
——-——
“-—‘
‘—’

_—v
—-———-——————-——-————-—--—
——-----——----———-—-—--——

‘
‘——
-——-——
-—-
-—‘-
--—
’—-
...‘C-

L

‘

_ _. :
... __. : :
_._.::
......I
......I
.. __ __ __
_. __ __ _.
L _“ L.
_
I .. ,__
_ _ ._
_ _
:0): -I
:2 II

‘

 

-----n-fl-——----——--—-——_q—-
...—q—
——-—
——‘
---—‘—
-u."
'—

r

_‘

 

¢
—-—-—----—---—————---—_—

     

      

 

--—--—--—————-‘--—------q

’—
-‘
u—c—
——‘—-
—’-—---
——
“—
-—’
P’
---—----—-—-—

 

(1) “N uaqumu 11%an

48
the time when the fully developed oscillatory pattern has settled down, the frequency

of this oscillation is the same as that of the surface temperature, though there usually
is a phase shift.

2. The average of the oscillation is done over one surface temperature period.
The average will approach a constant value as the fully developed oscillation is
reached. If the average temperature of the driving force is the same as that of a con-
stant-surface temperature case, the average Nusselt number is higher than that of
the constant surface temperature case. For example, for Ra=7474, the average is

Nu01=11.3, while for the constant surface temperature case, the Nusselt number is
NUC’1=5.3.

For the purpose of investigating the frequency effect on the Nusselt number,
some of the calculation results are summarized in Figure 11. A new non-dimensional
parameter N is defined as the following:

Nuo,l

= (32)
Nuc’l

 

The Nusselt number represents the non-dimensional heat transfer coefficient. The

01 and II

subscripts " denote surface-average Nusselt numbers for oscillating and

c,l

constant surface temperature cases, respectively. The "-" over Nuol denotes the av-

eraged value over an oscillation period and the measurement of the heat transfer en-
hancement is in the sense of the time-averaged Nusselt number. As is obviously
seen in Figure 11, the number N is always higher than one as long as the frequency is
non-zero. This shows that an oscillation in surface temperature results in increased

overall heat transfer effect. The heat transfer enhancement increases as the frequen-

49

an? ,.un.m<.
3259:25 .825... .3: 5 55:3...— ue Beam 3 flaw:—

ANE u .xocoscofl coca—=80 95.8380... Sufism

comp 00.0— 0mm

L I—

0.0

 

 

ITO

de

Fm.—

IQ.—

fioN

 

TN

 

N ‘roqumN luamaouequa ragsuml JeaH

50
cy increases but finally reaches a plateau after it reaches a value of 2.35. This indi-

cates that the effect of enhancement of heat transfer by varying the surface tempera-
ture is limited for a given average surface temperature. After the maximum is
reached, increasing the surface temperature frequency will not contribute any more to
the enhancement of heat transfer.
4.5 Effect of oscillation amplitude on the Nusselt number

As the frequency increases the heat transfer coefficient, the amplitude of the
surface temperature oscillation will also affect the enhancement of heat transfer.
When the amplitude is low, the heat transfer is close to that of the constant-surface
temperature situation. The heat transfer then increases as the amplitude increases,
at least in a certain range of amplitude. To test out this theory, several cases, with

Ra=7474, ASP=1, AX=1/40, AY=1/4O and At=1/10, 9 =1, and Fre.=10 Hz un-

ave
changed, were calculated, and the results are shown in Figure 12. In Figure 12, N
versus A6 is plotted. N is as defined in equation (32). Two observations are made
from Figure 12.

1. The heat transfer number, N, is always greater than one. This again verifies
that oscillation in surface temperature will increase heat transfer coefficient.

2. The effect of enhancement of heat transfer increases very rapidly at the small
amplitudes and reaches an upper limit around a value about 2.35 . After the plateau is
reached, the increasing in amplitude will not result in greater effect of heat transfer en-
hancement.

4.6 Effect of aspect ratio on the Nusselt number

The effects of both frequency and amplitude on the enhancement discussed pre-

viously are for the case when the aspect ratio is large enough so that the two surfac-

es are basically independent of each other. What will happen in terms of heat trans-

51

o—

m
_

F

o
_

Anmoenm Jun—93
EoEoogscom 3,—2.8. .8: .5 «gigs he .35— N— «....wE

c<

k

e
_

b

N
_

 

 

 

 

52
fer, if the surfaces are quite close to each other? A series of calculations were carried

out, with Ra=7,474, AX=1/40, AY=1/4O and Ar=1/10, 9 =1, A9=O.5, and Fre.=10

ave
Hz unchanged, but the aspect ratio decreases from two to close to zero. Figure 13
gives the result. When the aspect ratio decreases from two to one the heat transfer
coefficient almost does not change. This implies that for the aspect ratio greater than
one, the simulation is for single plate. When the aspect ratio decreases from one to
about 0.8, the heat transfer starts to increase, this is quite like the chimney effect for
the two constant surface temperature situation. The chimney effect will reach a peak
and then get weaker. As the aspect ratio get too small, say 0.2, conduction starts to

become the dominant heat transfer mode.

53

Guano £52":
3259558 cam—.2. .8: .5 83.. Bonus ac «ovum mu 8:»:—

a?
m; N4 mwo

b b P P — n b

Yo
_ _

b

 

 

 

 

54
5. CONCLUSIONS

A numerical study of natural convection heat transfer caused by a time-oscil-
lating surface temperature has been carried out using the UNDSAFE program. The
following conclusions are supported by the above study.

1. An oscillating surface temperature causes a natural convection flow which
is very different from the conventional pattern adjacent to a steady surface tempera-
ture surface. This was termed "thermal turbulence", and it appears in both tempera-
ture distribution and flow pattern;

2. The heat transfer enhancement has a strong relation to the oscillation fre-
quency. The increase of the heat transfer was seen to becomecd smaller as the fre-
quency was increased. There appears to be a limit for the heat transfer enhancement
for a specific average surface temperature and oscillation amplitude;

3. The oscillation amplitude has a significant effect on the heat transfer. At
high amplitudes, the heat transfer plateaus;

4. In terms of increasing heat transfer, it is unnecessary to have very high fre-
quency and high amplitude.

The success in simulating natural convection heat transfer driven by oscillating
surface temperature exposes a new approach in electronic cooling technique. Since a
time varying surface temperature is easy to achieve and the oscillation in surface tem-
perature will increase the heat transfer significantly without additional hardware, it is
very suitable for electronic cooling. Also, as pointed out by the study, the oscillation
does not have to have too high frequency or amplitude, which is desired when the fa-

ti gue is concerned for materials inside electronic devices.

55
6. RECOMMENDATIONS

To continue the current study, the following questions can lead to some inter-
esting new directions:

1. At small aspect ratios, what is the interaction of the two buoyant flows
caused by two oscillating wall temperatures? It is obvious that when two plates,
each has a time oscillation in surface temperature, are close to each other, the two
natural convection flows caused by the two oscillating surface temperatures will inter-
act with each other. This will make the heat transfer even more complicated. If the
interaction will increase or decrease the heat transfer enhancement is yet unknown.

2. When the surface temperature has a time varying spacial variation, what is
the temperature field, flow pattern, and heat transfer rate? In the electronic devices,
the hot (or cold) spot locations may vary with time. As can be imagined, this varia-
tion will cause a new pattern of heat transfer and fluid flow, and the heat transfer rate
will be different from that of the constant surface temperature situation. The details
need to be studied.

3. When the incoming boundary conditions are different from quiet ambient con-
ditions, how will the temperature and flow fields as well as heat transfer rate be af-
fected? There are all kinds of entrance boundary conditions other than the discussed
ones in the current study. Different boundary conditions may request different model-
ing technique for treatment. Different boundary conditions will result in various types

of temperature and velocity fields. Heat transfer rate will be different, too.

APPENDIX

(Numerical Code)

C

C
C
C

56

DATE: JAN. 1, 1989

NAME: WEI CHA

THIS PROGRAM IS DEVELOPED BASED ON UNDSAFE. IT IS USED TO SIMULATE NATURACL CONVECTION CAUSED

BY AN OSCILLATIN GSURFACE TEMPERATURE. I

LOGICAL‘I LBAND
COMMON/RU LBAND(9)
COMMON/R4/ X(450:90).Y(-60:90).DXX(-60:90),DYY(-60:90). DXXS(-60:90),DYYS(-60:90)
COMMON/BLlIDXDY,VOLDTMXOY.YOX,VOLDT,TI-IOT,PI
COMMON/BL7/NI,NIP1,NIP2,NIM1 MNPI ,NJP2.NJM1
COMMON/BLBMMLPINIBNIBPLNJTNTFMLNHPLNTMMNBPZ
COMMON/BLIZ/ NWRIFENTAPENTMAXONTREALTIMESORSUMII'ER
COMMON/31.16] CONST 1.CONS'I2,CONST3 ,RA.PR..NT.U0.H.UGRT,BUOY,
A: PSY,CPO.CONDO.VISO.RHOO,HR,TR.TA,DTEMP
COMMON/3L3]! TOD(-60:90.-60:90),ROD(-60:90,—60:90). UOD(-60:90.-60:90),VOD(-60:90,-60:90)
COMMON/3132! T(-60:90,-60:90),R(-60:90.-60:90),U(-60:90.-60:90),V(-60:90.-60:90),P(-60:90,-60:90)

& ,PSI(-60:90.-60:90)

COMMON/BIS?” TPD(-60:90.-60:90),RPD(-60:90,—60:90), UPD(-60:90,-60:90),VPD(-60:90.-60:90),

& PPD(-60:90.-60:90),PEQ(—60:90,~60:90)

COMMON/31.34] HEIGHT(-60:90.-60:90),SMP(—60:90,-60:90),SMPP(~60:90.60:90). DU(-60:90,-60:90),DV(-60:90,-60:90),
& PP(-60:90.-60:90)

COMMON/BU6/AP(-60:90,-60:90),AE(-60:90,-60:90),AW(-60:90.-60:90).AN(-60:90,-60:90),

& AS(-60:90.-60:90),SP(-60:90.-60:90).

& SU(-60:90.-60:90),REQ(-60:90,-60:90) ”

COMMON/BL37/ VIS(-60:90,-60:90),COND(—60:90,-60:90). CPM(-60:90.60:90),RESORM(20).TERM(20)

COMMON/BUG] ASXASY

DIMENSION RNU(1l),ANUC(3000),ANUl-I(3000),ANUM(3000),XT(3000) ,AFRE(3000),TEM(3000)
DATA XTIME.SORMAX,GRAV,GC,RAIR,1‘IMAXIO.,25.0 .2‘32.17,53.34.10l

C #1. READ IN DATA T0 INDICATE EITHER KRUN=O OR 1

READ(5.") KRUN.PR,TAO,N'IT
READ(5,‘)NM

CflREADINDATASETl-‘iDATA

READ(5,") NMAXNWRII'ENI‘APE
READ(5.‘) RA

NMN'I'I‘=NM‘N'I'I‘
W1.-3.‘((FLOAT(NI)—l.)/4.)
XNJB=1 .{FLOANND-l .)/4.

XNJT: (FIDATMI.)/4.+FLOAT(NIP1)
NIL=INT(XNIL)

NILPI=NIL+1

NILPZ=NIL+2

NJB=INTCXNIB)

NJBP1=NIB+1

NJBP2=NJB+2

NJT=INT(XNJT)

NJTM l=NJT-1

WI=N3T+1

NJ'I'P2=NJT+2

WRITE(6.‘) NIL,NIB,NJT

C ‘” INTRODUCE GIVEN PARAMETERS

DX=l .IFLOAT(NIM1 )
DY=l .IFLOAT(NIM1)
XI-I=0.0254
ALFA=2.25E-5
UO=ALFAIXH
T0=XHIUO

C ”’ GENERATION OF GRIDS

CALL GRID

C ”‘ INITIALIZE VARIABLE FIELD

=-l .5‘DYY(2)
DO 220 J=NJB,NJT
YY=YY+DYYS(I)
DO 220 I=NIL,NIPl
ROD(I,I)=1.
R(I,J')=l.
RPDGJ)=1-

57

PPD(I.I)=0.
DU(I,J)=0.
DV(I.I)=0.
SU(IJ)=O.
5mm.
mum.
AP(IJ)=O.
AWGJ)=0.
AB(I.I)=0.
AN(I,J)=0.
AS(I,J)=0.
swam.
swram.
VIS(IJ)=PR
COND(I,J)=1.0
CPM(IJ)=1.0
TOD(I,J)=-0.5+YY
“roomy-0.0
T(U)=1‘0D(U)
TPD(IJ)=TOD(IJ)

220 CONTINUE

TAO=TAOIIO

WTAOIFLOATM

TOU=TAO‘T0

DTIME=DTWDATW

TC=298.

DTHOT=250

THO= DTI'IOT+TAVE

TAVE=308.

1W0

G==9.8

MUE=1 589E-5
RA=G‘(XH”3)I(ALFA‘MUE)‘(2.‘DTIIOD/(THO-t-TC)

WRH'E(6,‘X3.XH.ALFA,MUE,DTHOT,'IHO,TC

D0 2m M=1,NIWNTI‘

2121 FORMAT(1X,'XH=’,F9.6,'TO=’,P‘9.4,’TAO=’,F9.6,'THO=’,F'7.3./.1X,

00000000 0

&'TC-—-’.F6.2,'N'IT=’,I4,'D'I'I'IME=’,F8.6.'D'I'IME=’ ,F9.6.‘NM= ’ ,IS

&J.’UO=’,F10.6,1X.‘ RA=',E10.4.‘ ASY/ASX=’,EI 0.4)

1H

PRINT‘,’RA=' ,RA

'I'I'IME=FLOAT(M)"DTI'IME
THOT=(1 .0-EXP(-'I'I'IME/I'AO))

THOT = l.0+SlN(2.‘3. 1415926‘1'I'IME/I' AO)‘DTHOT/(T AVE-TC)
THOT=l.
WRI'I'E(11,“)'I'I'IME,'IHOT
'I'I-IHOT=(298.15+5"I'HOT)
UE=14.58E-6"THHOT”1.5/(l10.4+THHOT)
CP=0.23 83 -0.79 l E-S‘THHOT+O.4834E-7‘THHOT‘ ‘2
XK=2.6482E-6‘THHOT”.5/(l .+245.4"'(1 0.“(- 12.miHOD)miIi01')
AI.FA=XK/CP/10.
RA=G‘(XH“3)/(ALFA"'UE)‘(2."DTHOT)/(THO+TC)

TCOOL=0.0

DO 564 I=1,NIP1

TUD(I,NJP1)=THOT

TOD(I,1)=TCOOL

564 CONTINUE

DO 9080 I=NIL1

58

TOD(I.NJBP1)=0.0
TUDWlHD
9080 CONTINUE
DO 9090 I=NJB,NIT
9090 CONTINUE
N1MAXO=0
C “" FOR CONTINUING RUN, READ DATA FROM TAPE
IF(KRUN ..EQ 1) GO TO 9997
GO TO 19
9997 READ(8.ENIk9998)
& TIME,N'IMAXO,'IOD,ROD,UOD,VOD,POD,CPM,COND,VIS,ITERT,
&M.NI,NIPI,NIPI,NIM1.NJM1,
& XXXXJBAND
GO TO 9997
9998 CONTINUE
REWIND 8
19 CONTINUE
BUOY=GRAV‘HI(UO* U0)
IX) 229 I=1,NIP1
DO 229 I=1,NIP1
REQGJ)=1.0
IF(KRUN .NE. 0) GO TO 229
RPDGJFREQM
RODGJFRPDGJ)
RWRPDUJ)
229 CONTINUE
228 CONTINUE
C ”" INITIALIZE U,V,T,R,P FIELD
DO 210 J=1,N.IP1
DO 210 I=1,NIPI
T(UFTODGJ)
R(LD=R0D(U)
UWUODGD
VGJFVODGJ)
PWPODCU)
210 CONTINUE
DO 211 I=NIB.NIT
DO 211 I=NIL,0
T(UFTODG-D
RWRODGJ)
UaJ)=UOD(IJ)
V(U)=V0D(U)
HUFPODGJ)
211 CONTINUE
NT =0
NTIM=0
300 CONTINUE
NT =NT+1
C IF(NT.GT.15)DTIME=0.(XX)005
C "” NTMAXO HAS THE MEANING AS "NTREAL" WHEN IT IS READ FROM
C DISK OR TAPE.
IF(NT.GT.NMAX) GO TO 303
NTREAL=NT+NTMAXO
TIME=TINIE+DTWE
551 D0 5551=1,NIP1
555 T(I,NIP1)=THOT
C FOR THE FICTIONAL WALL TEMPERATURE(5311,5312)
DO 5311 J=NJB,NJT
IF(U(NILJ).GT.O) T(NI1.,I)=TCOOL
531 1 CONTINUE
DO 5312 I=NIL.0
IF(V(I.NJE ).GT.0)T(I,NIB )=TCOOL
IF(V(I,NJT).LT.0) T(I,N1T)=TCOOL
5312 CONTINUE

 

59

C FOR THE CORNER TEMPERATURE
IE (U(NIL,NTB).GT.0 .AND. V(NIL.NIB).GT.0) T(NIL,NIB)=TCOOL
IF(U(NII.,NJT).GT. 0 .AND. V(NIL,NJT).LT.0) T(NIL.NIT)=TCOOL
C ”‘ START CALCULATION
ITER=0
JTERM=0
LITERM=0
C "” DEFINE THE UPDATED TPD(IJ). RGJ), UPDaJ) AND VPD(I,J)
C NR THE USE OF CALVIS AND SU(IJ)
DO 48 J=I,NIP1
DO 48 I=1.NIP1
1?me
RPMJFROJ)
UPDGJFUGJ)
‘0’me
48 CONTINUE
DO 499 I=NIB,NIT
DO 499 I=NIL,0
1?me
UPDGJFUGJ)
‘0’me
RPDGJFRGJ)
499 CONTINUE
29 CONTINUE
JTERM=ITERM+ l
C mu emu-tumult
CALL CALT
C to. accentuate.
DO 2“” I=NIB,NIT
DO 2000 I=NIL.NIP1
C IF(T(IJ).LE.TCOOL) T(IJ)=TCOOL
C IF(T(IJ).GE.'I'HUI') T(IJ)=THOT
2W CONTINUE

C ““ START PRESSURE CORRECTION ITERATIVE LOOP, IT IS THE MAJOR

C PART OF THE ERROR CONTROL ROUTINE
301 CONTINUE
ITER=ITER+ 1

C .0. 000......Ctt0.

CALL CALU

C it. ......‘tttt...

CALL CALV

C ‘3. ......COOOOQOO

CALL CALP
C ... ttttttttttttt.

mmfisomam) .LB. SORMAX) GO TO 49
IFaTER .EQ. 1) GO TO 302
lTERMl=ITER-l
IF(RESORM(ITER) .ua RESORMGTERM1))GO TO 302
GO TO 304

302 IFOTERM .GE. 2) GO TO 37
SOURCE:RESORM(ITER)
GO TO 39

37 IF(RESORM(ITER) .LE. SOURCE) GO TO 38
GO TO 304

38 SOURCE=RESORM(ITER)

39 CONTINUE
D0 23 1:1,Nm
D0 23 1:1 .NIPl
WDGJ)=T(IJ)
RPDGJ)=R(IJ)
UPDOJ)=U(U)
VPDCIJ)=V(IJ)
PPDGJ)=P(IJ)

23 CONTINUE
DO 233 J=NIB.NJT

60

DO 233 I=NIL,0
TPDGJPTGJ)
UPDGJFUGJ)
VPDGJFVGJ)
RPDGJFROJ)
PWPGJ)
233 CONTINUE
HTERM=0
IE(ITER .EQ. WAX) GO TO 49
IFOTERM .EQ. 2) GO TO 35
IFCITER .EQ. 4) GO To 29
35 CONTINUE
IPOTERM .EQ. 3) GO To 58
IF(ITER .EQ. 7) GO To 29
58 CONTINUE
HTERM=O
GO TO 301
304 CONTINUE
HTERM=HTERM+1
IEUTI'ERM.EQ.1)WRI'I'E(6.95) ITERRESORMGTER),SORSUM
IFOTERM .EQ.1)GO TO 41
IEOTERMHEQ 2.AND.IJTERM..EQ1.AND.ITER.NE. 5)GOTO41
GO TO 82
41 CONTINUE
DO 40 I=1.NIP1
DO 40 I=1,NIP1
R(I.J)=RPD(U)
UCUFUPDGJ)
VUJFVPDGJ)
MUFPPDGJ)
40 CONTINUE
DO 34 I=NIB.NIT
DO 34 I=NIL.0
RWRPDGJ)
UGJFUW
V(U)=VPD(IJ)
HUFPPDa-J)
34 CONTINUE
IFaTER ..EQ I'IMAX) GO TO 49
GO TO 29
82 CONTINUE
DO 43 J=1.NIP1
DO 43 I=1,NIP1
T(UFTPDGJ)
RGJFRPDOJ)
”(UFUPDGJ)
VGJFVPDGJ)
PWPPDGJ)
43 CONTINUE
DO 32 J=NIB,NJT
DO 32 I=NIL,0
T(IJ)=TPD(IJ)
U(I,J)=UPD(I.I)
VGJFVPDGJ)
RWRPDGJ)
POJFPPDGJ)
32 CONTINUE
IF(ITER ..EQ ITMAX) GO TO 49
IF((JTERM .EQ. 3 .AND. ITER .NE. 8) .OR. JITERM ..EQ 2) GO TO 49
GO TO 301
49 CONTINUE
ITERT=ITERT+ITER
C no :t:::;::::::tt:::::t::t:t:

CALL NU(RNU .RNUCRN UH)

... AAAAAAAAAAA ALAAA AAAAAAAA A 4
r v v v rv y—vvv—v—rv-v-Vvv—r—v rvvvrrw-

 

 

61

WRITE(13,") TIMERNU H
WRITE(14.‘) TIME,THOT
C ”" PRINT OUTPUT
C ”‘ SORSUM IS THE SUM OF ERROR SOURCE "8MP” FROM ALL OF THE CELLS
C IN THE ENCLOSURE.
IF (MOD(NT, 10).NE.0) GOT O 2433
2433 CONTINUE
IF(MOD(NT.9).NE.0) GO TO 2322
IP=1P+1
ANUC(JP)=RNUC
ANUHOPFRNUH
ANUM(IP)=RNU(6)
AFREOPFIO’THOT
XTOP)=NTREAL
IFOP.EQ.1) GO TO 2435
IFOP.GT.101) GO TO 2430
AAS=AAS+ABS(ANUC(IP)-ANUC(IP-l))
GO TO 2435
2430 AAS=AAS+ABS(ANUC(IP)-ANUC(IP~1))-ABS(ANUC(IP-100)-ANUC(IP-101 ))
IF(AAS/ANUC(IP).LE..(X)1)GO TO 276
AASS=AASIANUC(IP)
GOTO 2435
2435 CONTINUE
2322 CONTINUE
500 FORMATOX, ‘NTREAL=’.I9.1X,
& ’I'I'ER.=’,12J .'SOURCE=’,
& P9.6,1X,'SORSUM='B.6,1X,'NUC=’,F8.4,1X.’NUH=',F8.4,1X,
&/,11(1X.F6.3),' THOT:',E10.4J)
C un- fitttttttttttttttttttttttttfitmitt.
c CALL 'I'IEFI'GI')
123 FORMATC ITLEFI‘ = ’,110)
IID=IT
IF(IT.LT.ITLEFT) GO TO 277
C no ::::‘ :: :;::‘ttt*‘:t::t::‘ ::::::::
IF(NTREAL .NE. N'IREAIJNW'RI'I'E‘NWRITE) GO TO 505
GO TO 277
276 WRITE(6,11 12) NTREAL
1112 FORMATC THE STEADY STATE HAS BEEN REACHED AT NT = ',I9)
277 D0 502 II=1,NIP1f2
I=II‘2
513 D0 503 II=1NPIR
I=2"II
XTEMP=T(IJ)
=R(U)
XU=U(I,J)
=V(IJ)
XP=(P(U))
XVIS=VIS(IJ)
503 CONTINUE
502 CONTINUE
505 CONTINUE
C IF(IP.GT.101.AND.AAS/ANUC(IP).LE..001)GO TO 166
C “" RESET THE OLD TIME VALUES TOD, ROD, UOD, VOD AND POD.
DO 305 J=1,NJP1
DO 305 I=1,NIP1
T‘ODGJFTGJ)
RODOJFRGJ)
UODGJFUGJ)
VOD(I,J)=V(I,J)
PODUJFPOJ)
305 CONTINUE
DO 306 J=NIB,NIT
D0 306 I=NIL,0
mDGJFTGJ)
RODGJ)=R(I,J)

 

62

UODGJFUGJ)
V0D(LD=V(LD
PODGJFPGJ)
306CON'1'INUE
C fit. ::,::::CC::C:C3C‘::::::::::::::C:Ciiiifitiffiii :1: . 11:1;13111
IF(NTREAL .NE. N'I‘REAUNTAPE‘NTAPE .AND. (XTIME+DTIME*HIUO)
3 .LE. TMAX) GO To 522
WRTI'E(9)
a IMMREALTRU.V,P,CPM.COND.VISIFERT.
a H.TA,UO,CONDO,VISORHOO.NI.NJ.NIPINIPI,NIM1.NIM1,TIMES,
a WH20,wcoz.x.Y,Dxx.DYY,D)Ots,DYYs.LBAND
REWIND 9
c CALI. TLEETm')
IF(1'1'.LT.I'I'LEF1') GO To 166
C it. ..........O...‘OO‘.¢¢.‘.O......fittitttttttttittttttitltitttt
522 CONTINUE
C 1“ AC::3::1::iiiiitiiiiiiiiiiiiiii1:
c CALL 11mm)
IP(I'1‘.LT.1'I‘LEF1') GO TO 166
C it. iii: 1:1: 1 . Zitfi‘ *i‘ :*‘ “‘tttittttt
GO To 300
303 CONTINUE
WRITE(6,1111)
C WRIFE(11.1111) 1001 FORMAT(1X,' M=’.13)
1111 FORMATax.'—”m THE MAXIMUM TIME HAS BEEN REACHED maria)
GO TO 172

 

 

 

 

166 IF(N'IREAL .NE. NTREAIJNTAPE‘NTAPE) WRITE(9)
& TIMENI'REAI..T,R,U,V,P,CPM ,COND,VIS,1TERT.
& H,TA.UO.CONDO.VISO,RHOO,NI,NI,NIP1.NIP1,NIM1,NIM1,TIME3.
A: WH20,WCO2,X,Y.DXX.DYY.DXXS,DYYS,LBAND
REWIND 9
C it. :1:t?:::::i::tt‘t:::iittttttttt:3::::::C:tf::ti:::i :t*::::ii
172 CONTINUE
DO 152 I=1JP
BNUH=BNUH+ANUH(I)IIP
BNUM=BNUM+ANUM(I)/JP
BNUC=B NUC+ANUC(I)/JP
152 CONTINUE
D0 155 I=2,N'I
XX=X(I)+DXX(I)/2
156 FORMAT (2XJ4,3X.F8.4,10X,F10.4)
155 CONTINUE
CALL CALPL
IF(M. EQ.550 )GO'IO 319
IF(M. EQ.553 )GOTO 319
IF(M.EQ.555 )GOTO 319
IF(M. EQ.557 )GOTO 319
IF(M. EQ.54O )GOTO 319
IF(M. EQ.543 )GOTO 319
IF(M.EQ.545 )GOTO 319
IF(M. EQ.547 )GOTO 319
GOTO 2222
333 continue
2222 CONTINUE
STOP
END

C 1‘. 0.0....‘flittttfifittfitfittt..‘tfittttttt

SUBROUTINE CALT

C fit. .QOQQO..‘fitfittttttttfittfifitfifitttttit.

 

COMMON/RM X(-60:90).Y(-60:90).DXX(-60:90).DYY(-60:90), DXXS(—60:90),DYYS(-60:90)

COMMON/BL1IDX,DY.VOLDTIMEXOY.YOX,VOLDT,THOT.PI
COMMON/BL7/NINIP1 .NIPZNIMI ,NI,NJP1 ,NJP2,NIM1
COMMON/BL8/NILNILPI ,NIBNIBPlNJT,NJ'I'M1,NI'I'Pl,NI'I'P2,NILP2,NIB P2

63

COMMON/BLIZI NWRI'I'ENI’APENTMAXON’I'REAL.’IM,SORSUMII'ER

COMMON/ELM] CONS'I'1.CONST2,CONST3.RA.PR .NT.U0.H.UGRT,BUOY,

& PSY,CPO.CONDO,VISO,RHOO,HR,’I‘R,TA.DTEMP

COMMON/81.31] TOD(-60:90,-60:90),ROD(-60:90.-60:90). UOD(-60:90.-60:90),VOD(-60:90.-60:90)
COMMON/BL32/ T(-60:90,—60:90),R(-60:90.-60:90), U(-60:90,-60:90)N(-60:90,-60:90),P(—60:90,-60:90).

& PSI(-60:90,-60:90)

COMMON/BUN TPD(-60:90,-60:90),RPD(-60:90,-60:90). UPD(-60:90.-60:90),VPD(-60:90.-60:90).
& PPD(-60:90.-60:90).PEQ(-60:90,-60:90)

COMMON/BUM HEIGHT(-60:90,-60:90),SMP(-60:90,-60:90),SMPP(-60:90,60:90).

& DU(-60:90,-60:90),DV(-60:90,-60:90),PP(-60:90,-60:90)
COMMON/BL36/AP(-60:90,-60:90),AE(—60:90.-60:90),AW(-60:90,-60:90).AN(-60:90,-60:90),

& AS(-60:90.-60:90).SP(-60:90.-60:90).SU(-60:90,-60:90),REQ(-60:90.60:90)

COMMON/8L3” VIS(-60:90,—60:90),COND(-60:90.-60:90), CPM(-60:90,-60:90),RBSORM(20),TERM(20)
C “" CALCULATE COEFFICIENTS

D0 100 J=NIBP1NJTM1

JP2=J+2

JPl =J+l

IMl=J-l

JM2=I-2

D0 100 I=NILPINI

IF((I.GE.1 .AND. J.LE.1).OR.(I.GEI .AND. J.GE.NJP1)) GO TO 100

IP2=I+2

IPl=I+l

IMl=I-l

IMZ=I-2

DXI=DXX(I)

DXMl=DXX(IMl)

DXP1=DXX(IP1)

DYI=DYY(I)

DYM1=DYY(IM1)

DYP1=DYYOP1)

DXEE=DXXS(IP2)

DXE=DXXS(IPI)

DXW=DXXS(I)

DXWW=DXXS(IM1)

DYNN:DYYS(JP2)

DYN=DYYS(JP1)

DYS=DYYS(I)

DYSS=DYYSOM1)

YOXEzDYJ/DXE

YOXW=DYJIDXW

XOYN=DXIIDYN

XOYS=DXIIDYS

VObDXI‘DYJ

VOLDT=VOIJDTIME

CN=V(I,IP1)‘DXI

CS=V(IJ)*DXI

CE=U(IP1J)‘DYJ

CW=U(IJ)‘DYJ

CONDN1=XOYN

CONDS1=XOYS

CONDE1=YOXE

CONDW1=YOXW

CEP=(ABS(CE)+CE)‘DXE/DXI/16.

CEM=(ABS(CE)-CE)‘DXEIDXPII 16.

CWP:(ABS(CW)+CW)‘ DXW/DXMI/lé.

CWM=(ABS(CW)—CW)'DXW/DXI/16.

CNP=(ABS(CN)+CN)’ DYN/DYI/ 16.

CNM=(ABS(CN)-CN)‘DYN/DYPl/16.

CSP=(ABS(CS)+CS)‘DYS/DYMl/16.

CSM=(ABS(CS)—CS)‘ DYS/DYJ/lé.

AEGJ)=o.5‘CE+CEP+CEM‘(1.+DXE/DXEE)+CWM"DXW/DXE

AW(I,J)= .5‘CW+CWP‘(1.+DXW/DXWW)+CWM+CEP‘DXEIDXW

ANGJh-J‘CNo-CNHCNM‘(1.+DYN/DYNN)+CSM‘DYS/DYN

AS(I,J)= .5'CS+CSP‘(1.+DYSIDYSS)+CSM+CNP‘ DYN/DYS

IF (LLT.NI) GOTO 801
AEE=0.

AEER=0.

GOTO 802

801 AEE=CW8DXEIDXEE
AEER=AEE8TPD(IP2J)

802 CONTINUE
IF (LGTZ) GOTO 803
AWW=0.

AWWR=0.
G010 804

803 AWW=-CWP‘DXW/waw
AWWR=AWW‘TPD(MJ)

804 CONTINUE
IF (ILTM) G010 805
ANN=0.

ANNR=0.
GOTO 806

805 ANN=-CNM‘DYNIDYNN
ANNR=ANN‘TPD(IJP2)

806 CONTINUE
IF (1 .GT .2) GOT O 807
ASH.

ASSR=0.
GOTO 808

807 ASS=6PDYSIDYSS
ASSR=ASS’TPD(IJM2)

808 CONTINUE
APWAEWAWMANM+ASOD +AEE+AWW+ANN+ ASS)+CONDE 1 +CONDW1 +CONDN1+CONDSI
AEGJFAEGJHCONDEI
AWOJ)=AW(U)+CONDWI
AN (”FAN WONDN 1
AS(U)=AS(IJ)+CONDSI
SWVOLDT
SUM): VOLDT‘TPDGJ)
SUWUWAEER+AWWR+ANNR+ASSR

100 CONTINUE

C 8” TAKE CARE OF B.C TI'IR AE.AW.AN,AS,....

C “" ISOTHERMAL WALLS ,FLOOR AND CEILING
DO 500 1:2,NI
SPGJFSPMMSGJ)
SUMFSU(L2)+2.’AS(I,2)‘TPD(LI)
SWPWANW
SUWUWZ‘ANGND'ITIXLNTPI)
AN (I,NJ)=0.

500 CONTINUE

C 8" ADIABATIC WALLS: RIGHT WALL
DO 600 I=2,NI

C SPQJFSP (2.1)4'AW(2J)
SWPWUHAENJ)

C‘ SU(2,J)=SU(2J)+2.‘AW(2J)*(0.0)

C AW(2,J)=0.
AE(NI.J)=0.

600 CONTINUE

Cotoootmsz WALL
DO 601 I=NIBP1 ,NJTMI
SW1 J)=SP(NILPI D-AWMPI J)
SU(NILP1 MUMPI J)+2."AW(NILP1 D‘TCOOL
AW(NILP I J)=0.

601 CONTINUE

C“" FICTION AL HORIZONTAL WALLS
DO 222 I=NILP1,1
SP(I,NJ'IM1)=SP(I.NJTM1)-AW(I,NJTM1)
SUGWI )fiUONJTMIH-Z.‘AN(LNJTMI)'TCOOL

 
 

65

AN(I,NJT)=O.
C CN(I.NJ'I')=R(I.NJ'I')‘DX’V(I.NTI')
C IP(CN .GT. 0.) GO TO 233
SP(I,NJ'BP1)=SP(I,NJBPI)-AS(I,NIBP1)
SU(I,NJBP1)=SU(I,NJBP1)+2.‘AS(I,NIBP1)"'PCOOL
AS(I.N3BP1)=O.
222 CONTINUE
C‘“ I-ICI'IONAL TWO SHORT WALLS
D0 444 J=NJBP1.1
SP(0J)=SP(O,I)-AE(0J)
SU(0J)=SU(OJ)+2.‘AE(OJ)‘TPD(0J)
AE(0,J')=0.0
444 CONTINUE
DO 445 J=NIP1.NITM1
SP(0J)=SP(0J)—AE(0,I)
SU(0,J')=SU(0J)+2.‘AE(OJ)‘TPD(O,I)
AE(0,J)=:0.
445 CONTINUE
C "" ASSEMBLE COEFFICIENTS AND SOLVE DIFFERENCE EQUATIONS
DO 300 1:230
D0 300 I=2.NI
APGJFAPGMPGJ)
300 CONTINUE
DO 301 I=N3BP1,NITM1
D0 301 I=NILPl,l
APGJFAPGJI-SPGJ)
301 CONTINUE
CALL ’I'RIDAG(NILP1NIBP1NI,NI'IM1.T,2.NJ,1)
CALL TRIDA (NMINJBPININJ’IMI,T.NILPI,0.1.NTP1)
CALL 'I'RIDAX(NILP1.NTBP1,NI,NTI'M1.T.2,NJ.1)
CALL TRIDAY(NILP1NJBP1NI.NI'IM1,TNILP1,0.1,NIP1)
DO 604 I=2,NJ
T(NIPl J)=T(NI,J)
604 CONTINUE
DO 605 I=NILP1.0
T(INIB)=T(I,NIBP1)
T(I,NJT)=T(INI'I'M1)
605 CONTINUE
RETURN
END
SUBROUTINE CALU
C OOOOCOOOOOOOOOOOOOO
COMMON/BL lIDXDYNOLDTIMEXOYXOXNOLDTIHOTJ’I
COMMON/R4] X(-60:90),Y(-60:90),DXX(—60:90),DYY(-60:90), DXXS(-60:90),DYYS(-60:90)
COMMON/BL7/NU‘IIP1N1P2MM1 ,NI,NIP1,NIP2,NJM1
COMMON/BLS/NILNILPI ,NJBMBPINJTNTI’MIMTPl ,NITPzNILPZNJBP2
COMMON/BL12/ NWRI'I'ENTAPENI‘MAXONI'REALJ'IMESORSUMII’ER
COMMON/ELM] CONSTl .CONS'I‘2,CONST3 ,RA,PR ,NT.U0.HPSY.CPO,CONDO,VISO,RHOO,HR,TR,TADTEMP
COMMON/8L3 l/ TOD(-60:90.-60:90),ROD(-60:90.60:90), UOD(-60:90,-60:90),VOD(-60:90,-60:90)
COMMON/BUZ/ T(-60:90,-60:90),R(-60:90.-60:90),U(-60:90,-60:90),V(-60:90.-60:90),P(-60:90.-60:90) ,PSI(—60:90,-60:90)
COMMON/BUN TPD(-60:90.-60:90),RPD(-60:90.-60:90), UPD(-60:90,-60:90),VPD(-60:90.~60:90),
& PPD(-60:90.-60:90),PEQ(—60:90.-60:90)
COMMON/BL34/ HEIGHT(-60:90,-60:90),SMP(—60:90.-60:90),SMPP(-60:90,60:90), DU(-60:90.-60:90),DV(-60:90,-60:90),PP(-
& 6090,6090)
COMMON/BL36/AP(-60:90,-60:90),AE(-60:90,-60:90). AWN-60:90:60:90),AN(-60:90.-60:90), AS(~60:90,-60:90),SP(-60:90.-
& 60:90).SU(-60:90,—60:90),REQ(-60:90,-60:90)
COMMON/BMW VIS(-60:90,-60:90),COND(-60:90,-60:90), CPM(-60:90,-60:90), RESORM(20),TERM(20)
C ‘“ CALCULATE COEFFICIENTS
DO 100 J=NJBP1NJTM1
IP2=J+2
JP1=I+ 1
JM 1=J- l
IM2-J-2
DO 100 I=NILPZ.NI
IF((I.GE.1 .AND. J.LE. 1) .OR. (I.GE.1 .AND. J.GE.NIP1)) GOTO 100

66

IP2=I+2

IP1=I+1

IMI=I-1

IM2=I-2

DXI=DXX(I)

DXM1=DXX(1MI )
DWDXXGMZ)

DXP] =DXX(IPI )

DYJ=DYY(J)

DYM I =DYYOM 1)

DYPI=DYYOP1)

DXEE=DXXS(IP2)

DXE=DXXS(IP1)

DXW=DXXS(I)

DXWW=DXXSGM I )
DYNN=DYYS(II’2)

DYN=DYYSOPI )

DYS=D YYSU)

DYSS=DYYSOM1 )
XOYN=DXWIDYN
XOYS=DXW/DYS

YOXE=DYIIDXI

YOXW=DYJIDXM 1

VOL=DXW‘ DYI
VOLDT=VOIJDTIME

GN= V0.1? 1 )

GNW=VWI JP 1)

GS: V(I,.I)

GSW=V(N I J)

GE=U(IPI J)

GP=U(U)

GW=U(IM1 J)

CN=(GN‘DXNI I+GNW"DXI)I(DXM I +DXI)‘DXW
C3=(GS’ DXM I +GSW‘ DXI)I(DXM 1 +DXI)‘ DXW
CE=.5" (GE+GP)‘DYI

CW=.5" (GP+GW)‘ DY]

VISN 1 =XOYN‘ PR

VISSI =XOYS‘ PR

VISE1=YOXE8PR

VISWI =YOXW‘PR

CEP=(ABS(CE )+CE)’ DXI/DXWII 6.
CEM=(ABS(CE)-CE)’DXI/DXE/l 6.
CWP=(ABS(CW)+CW)‘ DXM 1/DXWW/l6.
CWM=(ABS(CW)-CW)‘ DXMl/DXW/16.
CNP=(ABS(CN)+CN)* DYN/DYI/16.
CNM=(ABS(CN)-CN)‘DYN/DYP1/16.
CSE(ABS(CS)+CS)'DYS/DYM1/16.
CSNI=(ABS(CS)--CS)‘I DYS/DYJ/I 6.
IE (LLTM) GOTO 801

AEE=0.

AEER=0.

GOTO 802

801 AEE=-CEM‘DXI/DXP1

AEER:AEE‘UPD(IP2J)

802 CONTINUE

IF (LGT.3) GOTO 803
AWW=0.

AWWR=0.

GOTO 804

803 AWW=-CWP“' DXM l/DXMZ

AWWR=AWW‘ UPD(IM2J)

804 CONTINUE

IF (J.LT.NI) GOTO 805
ANN=0.
ANNR=O.

67

GOTO 8&5

805 ANN=-CNM’DYNIDYNN
ANNR=ANN‘UPD(UP2)

806 CONTINUE
IF (I.GT.2) GOTO 80‘!

ASS=O.
ASSR=0.
GOTO 808

807 ASS=-CSP"DYSIDYSS
ASSR=ASS‘UPD(UM2)

808 CONIINUE
AE(IJ)=-.5‘CE+CEP+CEM‘(1.+DXI/DXP1)+CWM‘DXMl/DXI+VISEI
AW(I,I)= .5‘CW+CWP‘(1.+Dml/DXM2)+CWM+CEP‘DXIIDXM1+VISW1
ANaJk-j I"CN+CNP+CNM"(1.+DYN/DYNNMCSM"DYSIDYN+VISN1
AS(I.I)= .5'CS+CSP"(1.+DYSIDYSS)+CSM+CNP‘DYN/DYS+VISSI
SP(I,I)=-VOLDT
APGJ)=AE(U)+AW(IJ)+AN(U)+AS(U)+

& AEE+AWW+ANN+ASS
C “‘ SU FROM NORMALSTRESS
SUGJ)=DYJ‘ (mun-90.1»
a +UOD(LJ)‘VOLDT+AEER+AWWR+ANNR+ASSR
& +RA‘PR*(T(IMID‘DXI+T(I,J)’DXMl)/(DXI+DXM1)*VOL
100 CONTINUE
C "‘ TAKE CARE OF B.C THR AE.AW.AN,AS.....
C “‘ FLOOR AND CEILING
D0 500 1:3,NI
SP(I.2)=SP(I.2)-AS(I.2)
SWPGMANGND
AS(L2)=0.
AN(LNI)=0.
500 CONTINUE
C “‘ RIGHT WALL
DO 600 1:2,NI
SU(NIJ)=SU(NIJ)+AE(NLD‘UPD(NIPI ,1)
W-
600 CONTINUE
C“” LEI-'1‘ WALL
DO 505 J=NIBP1,NJTMI
SU(NILP2,J)-JU(NILP2J)+AW(NILP2J)‘ UPD(NILP1,J)
AW(NILP2J)=O.
505 CONTINUE
C””” SOLID VERTICAL WALL
DO 501 J=NJBP1,1
AE(0.J)=0.

501 CONTINUE
DO 502 J=NJP1 ,NJTMI
AE(0,J)=0.

502 CONTINUE

C‘”“" HORIZONTAL FICIIONAL WALLS
DO 508 I=NILP1,1
SP(LNIBP1):SP(LNJBP1)—AS(LNIBP1)
SU(I,NJBP1)=SU(I,NIBP1)+2.’AS(I.NIBP1)‘UPD(I.NIB)
SP(I,NJTM1)=SP(INTI'M1)-AN(I.NITM1)
SU(I,NJTM1)=SU(I,NITM1)+2.’AN(I.NTI‘M1)"UPD(I,NJT)
AS(I,NIBP1)=0.
AN(I,NJTM1)=0.

503 CONTINUE

C “"' ASSEMBLE COEFFICIENTS AND SOLVE DIFFERENCE EQUATIONS
DO 301 J=NJBPLNJTM1
DO 301 I=NILP2,NI
IF((I.GE.1 .AND. J.LE.1) .OR (I.GE.1 .AND. J.GE.NIP1)) GOTO 301
DYI=DYY(J)
AP(LD=AP(U)-SP(IJ)
DU(IJ)=DYJIAP(IJ)

301 CONTINUE

68

CALL TRIDAG(NILP2,NIBP1 ,NLNITMI.U,2,NI,1)
CALL TRIDA (NILP2,NJBPI,NI,NITNII ,U,NILPZ.0,1,NIP1)
CALL TRIDAX(NILP2,NIEP1,NI,NITM1,U,2,NT,1)
CALL TRIDAY(NILP2,NJBPI,NI,NITMI,U.NILPZ,0,I,NIP1)

C WRITE(6,999) ((AP(IJ),I=1,NIP1)J=1,NIPI)

999 FORMAT (22F6.3)
DO 704 I=NLBP1 ,NTIMI
U(NILP1J)=U(NILP2J)

704 CONTINUE
DO 705 J=2,NI
Um I J)=U(NU)

705 CONTINUE
DO 7% I=NILP2,0
U(I,NIB)=U(I,NJBP1)
U(I,NIT)=U(I,NJ'IM1)

7w CONTINUE
RETURN
END
SUBROUTINE CALV

‘A‘AAAA‘AAAAAALA‘L‘LAAJAAAAAAA‘AAALAAAAAAAJAJ‘J‘ A .A ALLA

‘ w—vvrvvv—vvvvv—v-v-v—v—vvvw-vvvrvr r'vf‘vv-VTV‘VT‘VV'TTV‘varr rr-rvr-r

COMMON/BLIIDX,DY,VOL,DIIME,XOY,YOX,VOLDT,THOT.PI

COMMON/R4/ X(-60:90).Y(-60:90),DXX(-60:90),DYY(-60:90), DXXS(-60:90),DYYS(-60:90)
COMMON/BL7/NINIP1NIP2NIM1NJ,NIP1,NJP2,NIM1

COMMON/BLS/NILNILP1 ,NIB ,NIBPI ,NJT.NITM1,NJTP1 ,NITP2,NILP2,NIBP2
COMMON/BL 12/ NWRITENTAPENIMAXONTREALJIMESORSUMII’ER
COMMON/BL 16/ CONST 1 ,CONSI‘2,CONST3 ,RA,PR

COMMON/3L3 I/ TOD(—60:90,-60:90),ROD(—60:90.-60:90),

& UOD(-60:90,60:90),VOD(-60:90,-60:90).

& POD(-60:90.-60:90)

COMMON/BUN T(-60:90,-60:90),R(-60:90,-60:90).U(-60:90.-60:90),
&. V(-60:90,-60:90),P(-60:90,-60:90)

& ,PSI(-60:90,-60:90)

COMMON/BL33/ TPD(-60:90.-60:90).RPD(-60:90,-60:90).

J: UPD(-60:90,-60:90).VPD(-60:90.-60:90).

& PPD(-60:90,-60:90),PEQ(-60:90,-60:90)
COMMON/ELM] HEIGHT(-60:90,-60:90),SMP(-60:90,-60:90),
& SMPP(-60:90.-60:90),

& DU(-60:90,-60:90),DV(-60:90,-60:90),PP(-60:90.-60:90)
COMMON/BL36/AP(-60:90.-60:90).AE(-60:90.-60:90).

& AW(-60:90,-60:90),AN(-60:90.-60:90).

8: AS(—60:90,-60:90),SP(-60:90,-60:90),SU(-60:90,-60:90)
& ,REQ(-60:90.-60:90)

COMMON/81.37] VIS(-60:90,-60:90).COND(-60:90.—60:90),

& CPM(-60:90,-60:90).

& RESORM(20),TERM(20)
C ”" CALCULATE COEFFICIENTS

NJTM2=NJT-2

DO 100 J=NIBP2,NITM1

JP2=J+2

IP1=J+1

1M 1:] -l

IM2=J-2

DO 100 I=NILP1,NI

IF((I.GE.1 .AND. J.LE.1) .OR. (I.GE.1 .AND. I.GE.NIP1)) GOTO 100

IP2=I+2

IP1=I+1

IM1=I-l

IMZ=I-2

DXI=DXX(I)

DXM1=DXX(IM1)

DXP1=DXX(IP1)

DYJ=DYY(I)

DYM1=DYY(IM1)

DYM2=DYY(JM2)

DYP1=DYYOP1)

69

DXEE=DXXS(IP2)
DXE=DXXS(IPI)
wa=nxxsa)
DXWW=DXXS(IMI)
DYNN=DYYS(JP2)
DYN=DYYS(JP1)
DYS=DYYS(J)
DYSS=DYYS(IM1)
XOYN=DXIIDYI
XOYS=DXIIDYM1
YOXE=DYSIDXE
YOXW=DYSIDXW
VOL=DYS"DXI
VOLDT=VOUDIIME
GN=V(I,JP1)
GP=V(I,I)
GS=V(LIM1)
GE=U(IP1,J)
GSE=U(IP1,JM1)
GW=U(IJ)
GSW=U(I,JM1)
CN=.5"(GN+GP)‘DXI
CS=.5‘(GP+GS)‘DXI
CE=(GE" DYJ+GSE'DYMl)/(DYJ+DYM1)*DYS
CW=(GW*DYJ+GSW‘DYMl)/(DYJ+DYM1)‘DYS
CEP=(ABS(CE)+CE)’DXEIDXI/l6.
=(ABS(CE)-CE)’DXE/DXPl/l6.

CWP=(ABS(CW)+CW)‘DXW/DXMl/16.
CWM=(ABS(CW)-CW)‘DXW/DXI/l6.
CNP=(ABS(CN)+CN)* DYJIDYS/l6.
CNM=(ABS(CN)—CN)*DYJ/DYN/16.
CSP:(ABS(CS)+CS)‘DYMl/DYSS/l6.
CSM=(ABS(CS)—CS)‘DYMl/DYS/l6.
VISEl=YOXE*PR
VISW1=YOXW‘PR
VISN1=XOYN*PR
VISSl=XOYS‘PR
IF (LLT.NI) GOTO 801
AEE=0.
AEER=0.
GOTO 802

801 AEE=-CEM-DXE/DXEE
AEER=AEE*VPD(IPz.J)

802 CONTINUE
IF (LGT.2) GOTO 803
AWW=0.
AWWR=0.
GOTO 804

803 AWW=-CWP*DXWIDXWW
AWWR=Aww-VPD(IM2.J)

804 CONTINUE
n= (I.LT.NI) GOTO 805
ANN=0.
ANNR=0.
GOTO 806

805 ANN=-CNM‘DYJIDYP1
ANNR=ANN‘VPD(IJP2)

806 CONTINUE
IF (J.GT.3) GOTO 807
ASS=0.
ASSR=0.
GOTO 808

807 ASS=-CSP‘DYMl/DYM2
ASSR=ASS"VPD(I,JM2)

808 CONTINUE

70

AE(IJ)=-.5*CE+CEP+CEM"(1.+DXE/DXEE)+CWM‘DXW/DXE+VISEI
AW(I,I)= .5‘CW+CWP‘(1.+DXW/DXWW)+CWM+CEP*DXE/DXW+VISW1
AN(I.J)=-.5‘CN+CNP+CNM‘(1.+DYI/DYP1)+CSM"DYMl/DYI+VISN1
AS(I,I)= .5‘CS+CSP‘(I.+DYMl/DYM2)+CSM+CNP‘DYJ/DYMI+VISSl
AP(IJ)=AE(IJ)+AW(IJ)+AN(IJ)+AS(LI)+AEE+AWW+ANN+ASS
SP(LJ)=-VOLDT

SU(I.I)= AEER+AWWR+ANNR+ASSR

c .

SUGJFSUGJFDXI'GGJMITKLD)

& +VOD(IJ)‘VOLDT

100 CONTINUE
C ”" TAKE CARE OF B.C THR AE,AW,AN,AS,....
C ‘” FLOOR AND CEILING
DO 500 I=2,NI
AS(I,3)=O.
AN(I,NI)=O
500 CONTINUE
C “‘ RIGHT WALL
D0 600 J=3.NZI
SP(NIJ)=SP(NIJ)-AE(NU)
SU(NU)=SU(NIJ)+2.‘AE(NI.D‘V(NIP1 J)
AE(NIJ)=0.
600 CONTINUE

C‘” LEFT WALL
DO 505 I=NIBP2.NITM1
SP(NILPI,J)-£P(NILP1,J)—AW(NILP1,J)
SU(N1LP1,J)=SU(NII..P1,J)+2.‘AW(NILP1J)*VPD(N1LJ)
AW(NILP1J)=0.

505 CONTINUE

C“ SOLID VERTICAL WALLS
D0 501 I=NJBP2.1
SRIMKIJI-AHU)

AE(1J)=0.

501 CONTINUE
D0 502 I=NJP2.NJT
SP(1J)=9P(1J)-AE(1J)
AE(1,I)=0.

502 CONTINUE

0"“ FICTION HORIZONTAL WALL
DO 503 I=NILP1.1

C SKLNIBP2)=SP(I.NIBP2)-AS(ININP2)‘VPD(I,NIBP1)
SUaNJBP2)=SU(I.NIBP2)+AS(LNJBP2)‘VPD(INIBPl)
AS(I.NJBP2)=0.

C SP(I.NJIM1)=SP(I,NIIMl)-AN(I,NJTMI)‘VPD(I,NIT)
SUGNJ’I‘MIFSUONITMI )+AN(I,NJTM1)*VPD(I,NJT)
AN(I.NJTM1)=0.

503 CONTINUE

C ‘” ASSEMBLE COEFFICIENTS AND SOLVE DIFFERENCE EQUATIONS
DO 300 I=NJBP2,NITM1
DO 300 I=NILPLNI
IF((I.GE.1 .AND. I.LE.1) .OR. (I.GE.1 .AND. J.GE.NIP1)) GOTO 300
DXI=DXX(D
AP(LD=AP(U)-SP(IJ)

DV(I,J)=DXI/AP(U)
300 CONTINUE

C “‘ DV ON HORIZEN'I‘AL WALLS ARE ZERO
DO 304 I=2,NI
DV(I.NJP1)=0.

DV(1,2)=0.
304 CONTINUE

999 FORMAT (22F6.3)
DO 704 J=NJBP2.NJT
V(NIL,I)=V(NILP1.J)

704 CONTINUE

DO 705 J=2,NJ

71

V(NIP1J)=V(NU)
705 CONTINUE
DO 7061=NILP1.0
V(LNJBP1)=V(I,NJBP2)
V(I,NJTM2)=V(LNITM1)
706 CONTINUE

 

SUBROUTINE CALP

C Ottttttttttttttt
COMMON/BL1/DX,DY,VOL,DTIME.XOY,YOX,VOLDT,THOT,PI
COMMON/R4/ X(-60:90).Y(-60:90),D)O{(-60:90),DYY(-60:90). DXXS(-60:90),DYYS(-60:90)
COMMON/BL7/NI.NIP1.NIP2,NIM1 ,NINIPI .NJ'P2,NJM1
COMMON/BL8/NE.NILP1NIBNIBP1,NJT,NIIM1.NIIPI,NITP2.NILP2.NIBPZ
COMMON/BL 12/ NWRITENTAPENTMAXONIREALTIMESORSUMJTER
COMMON/BL16/ CONST 1,CONS'I‘2,CONST3 ,RA,PR
& ,NT,U0.H,UGRT,B UOY,
8: PSY.CPO.CONDO.VISO,RHOO.HR,TR.TA,DTEMP
COMMONIBL31/ TOD(-60:90,-60:90),ROD(-60:90.-60:90).

& UOD(-60:90.-60:90),VOD(-60:90.-60:90),
& POD(-60:90,-60:90)

COMMON/31.32] T(-60:90,-60:90),R(-60:90,-60:90),U(-60:90,-60:90),
a; V(-60:90,-60:90),P(-60:90.-60:90)
& .PSI(-60:90,-60:90)

COMMON/31.33] TPD(-60:90.-60:90),RPD(-60:90,-60:90),

& UPD(-60:90.-60:90),VPD(-60:90,-60:90).

& PPD(-60:90,-60:90),PEQ(-60:90.-60:90)
COMMON/31.34] HEIGHT(-60:90.-60:90).SMP(—60:90,-60:90),
& SMPP(-60:90,-60:90),

& DU(-60:90,-60:90),DV(-60:90,-60:90),PP(-60:90,-60:90)

COMMON/BUélAH-éo:90,-60:90),AE(-60:90,-60:90).AW(-60:90.-60:90).
& AN(-60:90,-60:90),

& AS(-60:90.-60:90),SP(-60:90,-60:90),SU(-60:90.60:90)
& REQ(-60:90.-60:90)

COMMON/8L37/ VIS(-60:90.-60:90).COND(-60:90.-60:90).

& CPM(-60:90,-60:90).

& RESORM(40),TERM(40)

C “" CALCULATE COEFFICIENTS
DO 100 J=NJBP1NJTM1
JP2=J+2
IP1=I+1
JMI=J-l
IM2=I-2
D0 100 I=NILP1.NI
IF((I.GE.1 .AND. J.LE.1) .OR. (LGEJ .AND. J.GE.NIP1))GOTO 100
IP2=I+2
IP1=I+1
IM1=I-1
IMZ=I-2
DXI=DXX(I)
DXM1=DXX(IM1)
DXP1=DXX(IPI)
DYJ=DYY(I)
DYMI=DYY(IM1)
DYP1=DYY(IP1)
DXEE:DXXS(IP2)
DmDXXSGPl)
DXW=DXXS(I)
DXWW:DXXS(IM1)
DYNN=DYYS(JP2)
DYN:DYYS(JP1)
DYS=DYYS(I)
DYSS=DYYS(IM1)
DYS=.5"(DYJ+DYM 1)

72

VOIzDYI‘DXI
VOLDT=VOUDTIME
AN(I,.I)=DXI"DV(IJP 1)
AS(IJ)=DXI"DV(I,I)
AEGJ)=DYI‘DU(IP1J)
AWGJ)=DYI‘DU(IJ)
CN=V(I,JP1)‘DXI
CS=V(IJ)"DXI
CE=U(IPI D‘DYI
CW=U(IJ)’DYJ
SMPGJ)=-CE+CW-CN+CS
SUGJFSMPGJ)
SWIM.
1m CONTINUE
C "'" TAKE CARE OF B.C. THRU AN ,AS,AE,AW,SP AND SU
C ”" FLOOR AND CEILING
DO 500 I:2,NI
AS(I.2)=0.
500 CONTINUE
C “" RIGHT WALL
DO 611 J=Z,NI
AWWAWMAF-(NU)
611 CONTINUE
C"u LEFT WALL
DO 505 I=NIBP1,NITM1
AE(NILP1,.I)=AE(NILP1J)-AW(NILP1J)
AW(NILP1J)=0.
505 CONTINUE
C“ SOLID VERTICAL WALLS
DO 501 I=NJBPI.1
AE(1.J)=0.
501 CONTINUE
DO 502 I=NIPI,NJTM1
AE(1,J)=0.
502 CONTINUE
C“ HORIZONTAL FICTION AL WALLS
DO 503 I=NILP1,1
AN(I,NIBP1)=AN(I,NIEP1)-AS(I,NIBP1)
AS(I,NJTM1)=AS(I,NJTM1)-AN(I.NITM1)
AS(I,NIBP1 H.
AN(I.NITM1)=0.
503 CONTINUE
C "‘ ASSEMBLE COEFFICIENTS AND SOLVE DIFFERENCE EQUATIONS
DO 300 I=NIBP1,NJTM1
DO 300 I=NILP1,NJ
APUJFANWASWAEUD+AWGD5903
300 CONTINUE
CALL TRIDA (NII.P1.NIBP1,NI,NITM1,PP,NILP1,0,2,NI)
CALL TRIDAX(NILP1,NIBP1,NI.NITM1.PP2,NI,0)
C CALL TRIDAY(NILPI,NIBP1NI,NJTM1.PP,NILP1 ,0,2,NI)
C “‘ CORRECT VELOCITIES AND PRESSURE
C ”‘ CORRECTION FOR VELOCITY U
DO 600 I=3.NI
IM I=I-1
DO 600 J=2,NI
UGJFUGDTDUOD’CPPGMI JI'PPOD)
600 CONTINUE
C “" CORRECTION FOR VELOCITY V
DO 603 I=3,NJ
.IM 1=J~1
DO 603 I=2.NI
V(IJ)=V(IJ>*DV(IJ)‘(PP(IJM 1)-PP(I,J))
603 CONTINUE

73

C ”‘ CORRECIION FOR PRESSURE P
D0 606 I=2,NI
DO 606 I=2,NI
PGJ)=P(U)+PP(U)
PP(I,J)=0.
606 CONTINUE
C ‘" RECALCULATE THE ERROR SOURCE AFTER CORRECTIONS OF U. V. P
SORSUM=0.
RESORMCTERFO.
DO 700 I=2,NJ
IP2=J+2
IP1=J+1
IMP-'1 -1
IM2=J-2
DO 700 I=2,NI
IP2=I+2
IP1=I+1
1M l=I-1
IM2.=I-2
DXI=DXX(I)
DXM1=DXX(IM1)
DXP1=DX'X(IPI )
DYJ=DYY(J')
DYM1=DYY(IM1)
DYP1=DYYOP1 )
DXEE=DXXS(IP2)
DXE=DXXS(IPI)
DXW=DXXS(I)
DXWW=DXXS(IM1)
DYNN=DYYSOP2)
DYN=DYYSOP1 )
DYS=DYYS(I)
DYSS=DYYSOM1 )
VOL=DYJ‘DXI
VOLDT=VOIJDTIME
CN=V(I,JPI )‘DXI
CS=V(I,I)‘ DXI
CE=U(IP1 J')’DYI
CW :1] (LD‘ DY]
SMP(IJ)=-CE+CW-CN+CS
C ‘” SORSUM IS ACIUAL MASS INCREASE OR DECREASE FROM CONTINUITY
C EQUATUON . THIS WILL COMPARE TO SOURCE
SORSUM=SORSUM+SMP(IJ)
C ”‘ RESORM IS SUM OF THE ABSOLUTE VALUE OF SMPOJ) FOR ANY I
RESORMGTER)=RESORM(ITER)+ABS(SMP(IJ))
700 CONTINUE
RETURN
END
SUBROUTINE TRIDA(IBEGIN,JST ARTJENDJSTOPPHIIFROMJT OJ 3 .14)
C ......‘0.’00.0.00“0.0.0.000...000000000000O“..“‘t¢tt$fi00“.
COMMON/BL7/NI,NIP1,NIPZ,NIMI,NI,NIP1,NIP2,NJM1
COMMON/BL8/NII.NILP1,NIB,NIBP1.NJT,NIIM1,NITP1,N
COMMON/B L36/AP(-60:90.-60:90),AE(—60:90,-60:90),

& AW(-60:90,-60:90),AN(-60:90,-60:90).
& AS(-60:90,60:90),SP(-60:90,-60:90).SU(-60:90,-60:90)
& .REQ(-6O:90.-60:90)

COMMON/B L4I/ A(-30:100),B(-30:100),C(-30: 100),D(~30: 100)
DIMENSION PHI(-60:90.-60:90)
C “'" COMMENCE S-N TRAVERSE
D0 100 J=ISTARTJSTOP
JM1=J-l
IP1=I+1
C ‘” DEFINE ISTART AND ISTOP
IF((I .LE. J3) .OR. (1 .GE. 14)) GO TO 201
GO TO 202

74

201 ISTART=IFROM
ISTOP=ITO
GO TO 203

202 ISTART=IBEGIN
ISTOP=IEND

203 CONTINUE
ISTMl=ISTART-l
A(ISTM1)=0.
C(ISTMI)=0.

C ”‘ COMMENCE W-E SWEEP
DO 101 I=ISTART,IST OP
IMI=I~1
A(I)=AE(I.I)

B(I)=AW(I,J)
C(I)=AN(IJ)‘PHI(I.JP1)+AS(IJ)*PHI(LJM1)+SU(IJ)
D(I)=AP(U)

C ”" CALCULATE COEFFICIENT OF RECURRENCE FORMULA
TERM=1 J(D(I)-B(I)"'A(IM1))
A(I)=A(I)"TERM
C(I)=(C(I)+B(I)‘C(IMI))"TERM

101 CONTINUE

C ‘“ OBTAIN NEW PHI’S
PHI(ISTOP,J)=C(ISTOP)
ISTARI=ISTART+1
DO 102 II=ISTAR1.ISTOP
I=ISTUP+ISTART~II
PHI(IJ)=A(I)’PHI(I+1,I)+C(I)

102 CONTINUE
100 CONTINUE
RETURN

END
SUBROUTINE TRIDAGOS'TARTJBEGINJSTOPJENDIHIJFROMJIO ,I3)

4444 4444444444444444444

Vvvvvvvv‘f‘v vvvv‘rv v—vw—vvvv

COMMON/BL7/NI,NIP1,NH’2,NINII,NI,NIP1,NIP2,NJMI
COMMON/BLWIMBNIBPINTMTMIMTPIMTPZNHHNEPZ
COMMON/BL36/AP(-60:W,«60:90),AE(-60:90.$90).
& AW(-60:90,-60:90).AN(-60:90,-60:90),
& AS(-a):90,-60:90),SP(-60:90,-60:90),SU(-60:90.-60:90)
& ,REQ(-60:90,-60:90)
COMMON/BIA“ A(-30:100),B(-30:100),C(-30:1m),D(-30: 100)
DMENSION PHI(-60:90,-60:90)
C ”* COMNIENCE W-E SWEEP
DO 100 I=ISTART.IST OP
IP1=I+ l
IM1=I-1
C ”‘ DEFINE ISTART AND ISTOP
IF(I .GE. I3) GO TO 201
GO TO 202
201 JSTART=IFROM
ISTOHTO
GO TO 208
202 ISTARTQE EGIN
ISTOP=IEND
203 CONTINUE
JSTM 1 fl ST ART- 1
AUSTM 1)=0.
C(ISTM1)=0.
C ”’ COMMENCE S-N TRAVERSE
DO 101 I=JSTARTJSTOP
IM 1:] -1
AOFANGJ)
B(D=AS(U)
C(J)=AE(IJ)’PHI(IP1J)+AW(IJ)‘PHI(IM1,J)+SU(I,J)

DO)=AP(IJ)
C m CALCULATE COEFFICIENT OF RECURRENCE FORMULA

TERM=1J(D(I)~B(I)‘A(IM1)) 75
A(J)=A(.T)‘TERM
CO)=(C(I)+B(I)‘C(IM1))’TERM
101 CONTINUE
C “" OBTAIN NEW PHI'S
PHI(I.J STOP)=C(IST OP)
JSTAR1=JSTART+1
DO 102 JI=ISTAR1JSTOP
I flST OP+J ST ART-II
PHI(IJ')=A(J)‘PHI(IJ+1)+C(I)
102 CONTINUE
100 CONTINUE
RETURN
END
C :::::*::::::::::::::::::::::t::::::::::::::::::::::::::::::::.
COMMON/BL7/NI,NIP1,NIP2,NIM1.NI.NIP1.NIP2,NIM1
COMMON/BL8/NIL,NILP1 ,NTBNIBPINJTNJTM 1,NITP1.NJTP2,NILP2,NIBP2
COMMON/BL36/AP(-60:90,-60:90),AE(-60:90,-60:90),
«It AW(-60:90,-60:90),AN(-60:90.-60:90),
& AS(-60:90.-60:90).SP(-60:90.-60:90).SU(-&:90.-60:90),
& REQ(-60:90,-60:90)
COMMON/3L4” A(-30:100),B(-30:100).C(-30:100).D(-30: 100)
DIMENSION PHI(-60:90,-60:90)
C ‘” COMNIENCE W-E SWEEP
DO 100 I=IST ARTJST OP
IM1=IoI
IP1=I+1
C ‘”” DEFINE JSTART AND ISTOP
IF(I .GE. 13) GOTO 201
GOTO 202
201 ISTARTrJFROM
ISTOP=ITO
GOTO 203
202 JSTART=IBEGIN
ISTOHEND
203 CONTINUE
181‘ GP] =IST OP+1
BOSTOPI)=0.
C(I STUP1)=O.
C ‘” COMMENCE S-N TRAVERSE
DO 101 H=ISTARTJSTOP
I =1 ST OP+J ST ART-II
IP1=J+1
AC)=AN(I,J)
BG)=AS(IJ)
CO)=AE(IJ)"PHI(IPI J)+AW(I,J)‘ PHI(IM1,J)+SU(I,J)
DOFAPGJ)
C ”‘ CALCULATE COEFFICIENT OF RECURRENCE FORMULA
TERM=1J(D(I)-A(J)‘B(IP1))
BO)=B(I)‘TERM
C(J)=(C(J)+A0)‘C0P1))‘TERM
101 CONTINUE
C ”" OBTAIN NEW PHI'S
PHI(IJ START)=C(J ST ART)
JST AR 1 =1 ST ART+1
D0102 IdSTARlJSTOP
PmaJ)=B(D’Pm(U-1)tC(D
102 CONTINUE
100 CONTINUE
RETURN

END
SUBROUTINE TRIDAY(IBEGIN,J ST ARTJENDJ STOP,PHI.IFROM ,IT 0,] 3 ,J 4)

C‘ft‘iitttfittiiititiitiiC:ilfititttat aautttt‘tfiltifilttttttttttt

COMMON/BL?INI,NIP1,NIP2.NIM1,NI..\'J Pl NJPZNIMI

C

76

COMMON/BWIMBNIBPINJTNJTMINJTPINITPZWPZ
COWON/BL36/AP(-60:90,-60:90),AE(-60:90,-60:90),
& AW(-60:90,-60:90),AN(-60:90.-60:90).
&. AS(-60:90,-a):90),SP(-60:90.-60:90),SU(-60:90,—60:90),
& REQ(-60:90.-a):90)
COMMON/BLAH] A(-30:100),B(-30:100),C(-30:1(X)).D(-30: 100)
DIMENSION PHI(-60:90.-60:90)
C ”* COMMENCE S.N TRAVERSE
DO 100 JI=JSTARTJSTOP
J=ISTART+ISTOPJJ
IMlaI-l
JPld+l
C "'” DEFINE ISTART AND ISTOP
IF (( I.LEJ3) .OR. (I.GE. 14)) GOTO 201
GOT O 202
201 ISTART=IFROM
ISTOfiITO
GOTO 203
202 ISTART=IBEGIN
ISTOP=IEND
203 CONTINUE
ISTM l=ISTART-1
A(ISTM1)=0.
C(ISTMI)=0.
C "'“ COMMENCE W-E SWEEP
D0 101 I=ISTARTJST OP
IM1=I-1
A(I)=AE(IJ)
B(I)=AW(IJ)
C(I):AN(I,I)"PHI(IJP1)+AS(IJ)"PHI(IJM1)+SU(IJ)
D(D=AP(U)
C '” CALCULATE COEFFICIENT OF RECURRENCE FORMULA
TERM=l./(D(I)-B(I)‘A(IM1))
A(I)=A(I)‘TERM
C(I)=(C(D+B(I)‘CCM1))‘TERM
101 CONTINUE
C ”‘ OBTAIN NEW PHI’S
PHI(ISTOPJ)==C(ISTOP)
ISTAR1=ISTART+1
D0 102 II=ISTAR1.ISTOP
I=ISTUP+ISTART0II
PHI(IJ)=A(I)’PHI(I+1J)+C(I)
102 CONTINUE
100 CONTINUE
RETURN
END

4.4.44.44_44_4A44 LA‘A‘AAAAAAAAAAALJ 444..

Qt. ‘ U“
'73- 3-3333771333373333131,77 3,---

SUBROUTINE NU(RNU,RNUC,RNUH)

444.444LA4444 44.444.- 4 L4 A+A*A‘AA_LAA_AA 4444 4 4

COMMON/BLl/DXDY VOL,DTIME ,XOY. YOX, VOLDT, THOT. PI
COMMON/BL7/NI,NIP1,NIP2,NIMI ,NJ,NIP1 ,NIP2,NIM1
COMMON/BL8/NIL.NILP1,NJBNIBPINJTNIT‘MLNITPINITPZMLPZNIBW
COMMON/BL36/AP(-60:90,-60:90),AE(-60:90,-60:90),

& AW(-60:90.-60:90),AN(—60:90,-60:90),
& AS(-60:90,-60:90),SP(-60:90,-60:90),SU(-60:90.-60:90)
& ,REQ(-60:90,-60:90)

COMMON/BLIZ/ NWRITENT APENTMAXONTREALJIMESORSUMJT ER
COMMON/ELM] CONST 1,CONST2,CONST3 ,RA.PR
& ,NT‘,U0,H,UGRT,BUOY,
& PSY,CPO,CONDO,VISO,RHOO,HR,TR,TA,DTEMP

COMMON/BUZ/ T(-60:90,-60:90),R(-60:90,-60:90),U(—60:90,-60:90),

& V(-60:90.-60:90),P(-60:90,-60:90)

& ,PSI(-60:90,-60:90)

COMMON/3L3” VIS(-60:90.-60:90),COND(-60:90,-60:90),

& CPM(-60:90,-60:90).RESORM(20),TERM(20)

77

COMMON/R4] X(-60:90),Y(-60:90).DXX(-60:90).DYY(-60:90).

& DXXS(-60:90).DYYS(—60:90)

DIMENSION RNU(11)

DO 100 K=1, ll
RNU(K)=0.

100 CONTINUE

RNUH=0.

RNUC=O.

DT=T(2.NIP1)-1.

DO 20 I=2,NI
Y1=DYYS(2)f2.

Y2=YI+DYYS(3)

DIDYC=(Y2‘Y2"(T(L2)-T(Ll ))-Y1 ‘Y1 ‘ (T (I,3)-T(I, 0))” 2

& /(YZ-Yl)/Y1
Y1=DYYS(NIP1)/Z.

Y2=Y1+DYYS(NI)
UTDYH=—(Y2‘Y2‘(T(I,NI)-T(I,NIPI ))-Yl ‘Yl"(T(I,NJMl)-

& T(I,NIP1)))/Y2/(Y2.Y1)/Yl
IF(I.EQ.20) GOTO 999
GOTO IMO

999 XNU=DTDYH"DXX(I)
1(D0 DTDYCI =2. ‘(T (12)-Ta, l ))IDYYS(2)
DTDYHI=2.‘(T(I,NJP1)-T(I,NJ))/DYYS(1~UP1)

RNU(10)=RNU(10)+DTDYC1*DXX(I)

RNU(I l)=RNU(1 l)+DTDYHl"DXX(I)
RNUC=RNUC+DTDYC’DXX(I)
RNUII:RNUH+DTDYII‘DXX(I)

DO 10 K=l.9

J=2‘K+2
IMlaI-l
1P1 dd
GS: V(IJ)‘DIO((I)

GSP=(GS+ABS(GS))‘ DYYSONDYYOMI )
GSM=(GS-ABS(GS))"DYYS(J)/DYY(I)
Q=-5‘(T(IJ)+T(UM1))"GS

& -1./16.‘((T(U)-T(IJ'M1))

& {IGJMI}T(IJ-2))‘ DYYSOIIDYYSOMI)

& )‘GSP
Q=Q-1.II6.‘((T(IJP1)-T(IJ))‘DYYS(J)

& /DYYS(IP1)—(I'(IJ)-T(IJM1)))"GSM
DTDY=('T(IJ)-T(IJM 1))IDYYSO)
DQ=(Q-DTDY'DXX(D)

RNU(K)=RNU(K)-DQ

10 CONTINUE
C DLOC(I)=DTDYII
20 CONTINUE

RETURN
END
BLOCK DATA

C 0.0.00.000

LOGICAL‘I LBAND
COMMON/BL?M,NIP1.NIP2,NIM1,NJ,NIP1NIP2,NIM1
COMMON/BL8/NILNILP1 ,NIBNJBPI .NJT.NITM1.NITP1,NJTP2.NILP2,NJBP2
COMMON/BLI6/ CONST1,CONSI2.CONST3 ,RA,PR

a ,NT.UO.H,UGRT,BUOY,

& PSY,CPO.CONDO.VISO,RHOO,HR,TR,TA,DTEMP
COMMON/RU LBAND(9)

DATA NIP2,NIP2,NIP1 ,NIPI ,NINI ,NIMI.NJM1/43 ,43 ,42,42.41 ,4 1 ,40,40/
END

C .OOOOOOCOOOOODOO

COMMON/R4]X(-60:90).Y(-60:90).DXX(~60:90),DYY(-60:90),
& DXXS(—60:90),DYYS(-60:90)
COMMON/BL32/ T(-60:90,-60:90),R(-60:90,-60:90).U(-60:90,-60:90),

78

& V(-60:%,-a):90).P(-a):90.-60:90).
& PSI(-60:90.-w:90)
COWON/BL16/ CONST 1.CONST2.CONST3 ,RA,PR
& ,NT,U0,H,UGRT,BUOY,
& PSY.CPO,CONDO,VISO,RHOO.HR,TR.TA,DTEMP
COMMON/BL7/NI,N1P1,NIP2.NIM1.NI.NIP1 W1
COWON/BIS/NWINIBNIEPINITNTMIMTPIWWIBW
COMMON/BWAP(-60:90,-6090).AE(-60:90,-$:90),AW(—60:90,-60:90),
& AN(-60:90,-60:90),
& AS(-60:90,-60:90),SP(-6):90,-60:90),SU(-60:90,-60:90)
& ,REQ(-60:%,-60:90)
DIMENSION UU(-6090.-6090),W(-6090,-6090)
DO 307 I=I.NIP1
DO 307 J=1NP1
307 CONTINUE
DO 23 I=1,NIP1
DO 23 I=1NP1/2
C T(IJ)=2.0+DTEIHP-T(NIP1+1-I,NIPI+1-I)
23 CONTINUE
DO 24 I=1.NIP1
DO 24 I=2,NIP1I2+1
C V(IJ)=-V(NIP1+1-I.NIP1+2-I)
24 CONTINUE
DO 25 I=2,NIP1
DO 25 1:1,NIP1/2
C U(IJ)=—U(NIP1+2-LNIP1+1-J)
25 CONTINUE
C‘“”CALCU LATE STREAM FUNCTION DISTRIBUTIONS
DO 701 I=2,NI
DO 701 I=2,NI
PSIWFRBQQZ)‘(U(IJ)+U (1+ 1.2))‘ DYYQ)
IFO.LT.3) GO TO 701
D0 702 K=3J
702 PSIaJ)=PSI(IJ)+<RBQ(LK)‘(U(I.K)rU(I+l.K))‘DYY(K)
& +REQ(I,K-1)‘(U(I,K-1)+U(I+1.1(-1))‘DYY(K-1))‘0.5E0
701 CONTINUE
DO 1010 J=2,NIP1
PSI“ J)=PSI(2.D
PSI(NIP1 J)=PSI(NIJ)
1010 CONTINUE
DO 723 I=2,NIP1
DO 723 I=2,NI
723 CONTINUE
VMAX=0.
DO 501 I=2,NI
VV(I,1)=0.
VV(I,NIP1 )=0.
DO 502 I=2,NJ
WGJ):((V(U)*V(U+1))IDXX(I'1)+(V(I-1.D+V(I-1J+1))‘DXX(I))/
& (Z‘CDXXG'IHDXXGW
UUGJFUGJ)
VT=‘SQRT(UU(IJ)”2+VV(IJ)“2)
IF (VT .GT.VMAX) VMAX=VT
502 CONTINUE
501 CONTINUE
WRITE (6,1)
C DEFINE LEFT AND RIGHT WALLS FOR V
DO 503 I=2,NI
VV(2J)=VV(3J)
VV(NIP1 J)=VV(NI,I)
503 CONTINUE

  
 

79

WRITE (6,1)
GR=0.(X11
REL=0.(X)1
FRA=1.
c CALL CALCOM(UU,VV,T,PSI,VMAX,REL,H,UO,NI,NI.FRA,PSY)
RETURN

 

 

COMMON/R4/ X(-60:90),Y(-60:90),DXX(-60:90),DYY(—60:90)

& ,DXXS(-60:90),DYYS(-60:90)
COMMON/BL1IDX,DY,VOL,DTIME,XOY,YOX,VOLDT,THOT,PI
COMMON/BL7/NLNIP1NIP2NIM1MMPI MP2,NIM1
COMMON/BL8/NMILP1NIBNIBP1NJTNITM1NTTPINITP2NILPZMBP2
COMMON/BL12/ NWRITENTAPENTMAXONTREALJ'IMESORSUMJTER
COMMON/BL16/ CONST 1,CONST2,CONST3 ,RA.PR

& ,NT.U0,H.UGRT,BUOY,

& PSY.CPO,CONDO,VISO,RHOO,HR,TR,TA,DTEMP
COMMON/BUN TOD(-60:90,-60:90),ROD(-60:90,-60:90),

& UOD(-60:90,-60:90),VOD(—60:90,-60:90).

& POD(-60:90,—60:90)

COMMON/BL32/ T(-60:90,-60:90),R(—60:90,—60:90).U(-60:90,-60:90),

& V(-60:90.-60:90).P(-60:90.-60:90).

& PSI(-60:90,-60:90)

COMMON/BL33/ TPD(-60:90,—60:90),RPD(-60:90.-60:90),

& UPD(-60:90,-60:90),VPD(-60:90,-60:90).

& PPD(-60:90,-60:90),PEQ(-60:90,-60:90)
COMMON/BUM HEIGHT(-60:90,-60:90),SMP(-60:90,-60:90).

& SMPP(-60:90,-60:90),

& DU(-60:90.-60:90),DV(-60:90,-60:90).PP(-60:90.-60:90)
CONIMON/BL36/AP(-60:90.-60:90),AE(-60:90,-60:90),AW(~60:90,-60:90),

& AN(-60:90,-60:90),

& AS(-60:90,-60:90),SP(-60:90,-60:90).SU(-60:90,-60:90)

& REQ(-60:90,60:90)

COMMON/3L3” VIS(-60:90,-60:90),COND(-60:90,—60:90),
& CPM(-60:90,-60:90),RESORM(20),TERM(20)
COMMON/BL40/ ASXASY
C DEFINE ASPECT RATIO
ASX=1.0
ASY=1.00
C ‘” NX1,NY1 THE POINT IN BOUNDARY LAYER
C WRITE(11,1111)ASX,ASY
C1 11 FORMAT(1X,'ASX=’,F5.2,1X,'ASY=’,F5.2)
C ”‘ NX2,NY 2 THE POINT OUT BOUNDARY PAYER
NY 1 =5
NY2=NJM1f2-NY1
NX1=5
NX2=NIMl/Z-NX1
C “‘ DELXIDELYI BOUNDARY LAYER THICHNESS ASSUMED
C ”" DELXZDELYZ OUT OF BOUNDARY LAYER

DELYI =0.0650

DELY2=(.5-DELY1 )" ASY

DELX1==0.25(X)

DELY1=DELY1‘ASY

DELX2=O.5‘ASX-DELX1 "' ASX

DELX1=DELXPASX

NII =NIP1I2+ l
N12=NII+1
X(2)=0.

DO 11 I=3,NIl
DELX=DELX l/NX 1
IF (I.GE.(NX1+3)) DELX=DELX2INX2
X(I)=X(I- l)+DELX
ll CONTINUE

80

DO 10 I=N'12,NIP 1
K=NIP1-I+2
XCI)=1.0"ASX-X(K)

10 CONTINUE
X(1)=‘X(3)
X(NTP2)=2."X(NIP1)-X(NI)

C Y COORDINATE

N11 =NIP1/Z+1
WN11+ 1
Y(2)d).
DO 21 I=3,NII
DELY=DELY 1/NY1
IF (I.GE.(NY1+3)) DELY=DEI.Y2/NY2
YO)=Y(J-1)+DELY

21 CONTINUE
DO 20 I=NI2,NIP1
K=NIP1J+2
YU)=1.0‘ASY-Y(K)

20 CONTINUE
Y(1)='Y(3)
Y(NIF2)=2."Y(NTP1)-Y(NI)
WRITE (6,103) ASX,ASY

103 FORMAT (2X,' ASX= ',F8.5 ,' ASY=',F8.5)

C uninten- UNEORM GRID

DX= 1 JFLOAT(NINI1)‘ ASX
DY=1 JFLOAT(NTM1)‘ASY
DO 5 I=NIL.NIP2
X(I)=DX"(I-2)

5 CONTINUE
D0 6 I=NIBNITP2
YO)=DY"(I-2)

6 CONTINUE

C ”m“ END OF UNIFORM GRID
DO 9 1‘:le
IP1=I+ 1
DXX(I)=X(IP1)-X(I)

9 CONTINUE

DXX(NIP2)=DXX(NIPI)

DO 17 I=NILP1.NIP2

IM1=I-1
DXXS(I)=.5‘(DXX(I)+DXX(IM 1))

17 CONTINUE
DXXS(NIL)=DXXS(NILP1)
DO 14 I=NIB,NJT
IP1=I+1
DYYCI)=Y(IP1)-Y(I)

14 CONTINUE
DYY(NITF2)=DYY(NITP1)
DO 15 I=NIBP1.NITP2
IM1=J-1
DYYSO)=.5‘(DYYO)+DYY(IM1))

15 CONTINUE
DYYS(NIB)=DYYS(NIEP1)
DO 13 I=NIL.60

13 CONTINUE
RETURN
END

LIST OF REFERENCES

 

81

LIST OF REFERENCES

l. Nakayama, W., "Thermal management of electronic equipment: a review of technol-

ogy and research topics", Applied Mechanics Review, Vol. 39, No. 12, Dec. 1986.

2. Bergles, A.E., "The evaluation of cooling technology for electrical, electronic, and

microelectronic equipment", ASME Heat Transfer in Electronic equipment, 1986.

3. Steinberg, D.S., "Cooling techniques for electronic equipment", Wiley-Interscience

Publication ISBN 0-471-04403-2. 1980.

 

4. Elenbaas, W., "Heat dissipation of parallel plates by free convection," Physica,
Vol. 9, pP.1-28, 1942.

5. Sparrow, E.M., "Enhanced and local heat transfer, pressure drop, and flow visual-
ization for arrays of block—like electronic Components," Integpap'enel leprnel pf Heat
Transfer 3.1.151 Mass Transfer, Vol.26, No.5, pp.689-669, 1983

6. Bar-Cohen and Rohsenow, "Thermally optimum spacing of vertical natural convec—
tion cooled parallel plates," Journal of Heat Mass Transfer. Vol. 106, pp.116-123,
1984.

7. Park, K.-A. and Bergles, A.E., "Natural convection heat transfer characteristics of

simulated microelectronic chips," ASME Journal of Heat Transfer Vol.109, pp.90-96,

 

Feb. 1987

82
8. Wang, Y.Z. and Bau, H.H., "Low Rayleigh number convection in horizontal, eccen-

tric annuli", Proceeding pf ASME Annpel Meefing, Natupal and Mixed convection in
electrenie equipment eeqling, edited by Wirtz, R.A., Nov. 1988.

9. Torikoshi, K., Kawazoe, M. and Kurihara, T, "Convective heat transfer characteris-
tics of arrays of rectangular blocks affixed to one wall of a channel’, Proceeding of AS-

ME Annual Meeting, Natural and Mixed convection in electronic equipment cooling,
edited by Wirtz, R.A., Nov. 1988.

10. Arco, E. C., Bontoux, P., Ssni, R.L., Hardin, G., Extermet, G. P. and Chikhaoui, A,
"Finite difference solutions for three-dimensional steady and oscillatory convection in
vertical cylinders-effect of aspect ratio", Proceeding of SME Annual Meeting, Heat

Transfer Division Nov. 1988.

 

11. Ramanathan, S. and Kumar, R., "Correlations for natural convection between heat-

ed vertical plates", in Ann 1M 'n N r1 11 Mix onvec-

tien in eleetrenie eqpm’ men; emling, edited by Wirtz, R.A., Nov. 1988. Nov. 1988.

12. Bergman, TL. and Petri, J. G., "Natural convection cooling of discrete heated ele-

ments in an enclosure using pressure helium-xenon gas mixtures", Proceeding of

SME Annual Meeting, Heat Transfer Division, Nov. 1988.

13. Sparrow, E.M., Quack, H. and Boemer, C.J., "Local non-similarity boundary-layer
solutions," AiAA Journal, Vol.8, No.11, 1936.

14. Sparrow, EM. and Gregg, J.L.,"laminar free convection from a vertical plate with

uniform surface heat flux," Journal of Heat Transfer, pp.3435-440, 1956.

15. Sparrow, EM. and Yu, H.S., "Local non-similarity thermal boundary-layer solu-
tions," ASME Journal of Heat Transfer pp.328-334, Nov. 1971.

 

 

 

 

83
16. Kao, T.T, "Locally non—similarity solution for laminar free convection adjacent to a

vertical wall," ASME Jewel ef Heet Transfg, pp.321-322, May 1976.

17. Wirtz, RA. and Stutzman, R.J., "Experiments on free convection between vertical
plates with symmetric heating," ASME Journal of Heat Transfer, Vol.104, pp.501-
507, Aug. 1982

18. Hamady, F., "Experimental study of local natural convection heat transfer in in-
clined and rotating enclosures," Dissertation for Degr_ee of Ph.D., Michigan State Uni-
versity, 1987

19. O’Meara, M. and Poulikakos, D., "Experiments on the cooling by natural convec-
tion of an array of vertical heated plates with constant heat flux," Heat and Fluid-
Flow, Vol.8, No.4, Dec. 1987

20. Lloyd, JR. and Yang, K.T., "Technical memorandum on accuracy, convergence and

stability of UNDSAFE," NSF GRAND 67-9002 May 1978

 

21. Doria, M.L., "A numerical model for the prediction of two dimensional unsteady
flows of multi-component gases with strong buoyancy effects and re-circulation,"

NSF RANN GRANT GI-37191 and ATA 73-07749A01, Nov. 1974

22. Lloyd, J .R., Yang, K.T. and Liu, V.K., " A numerical study of one dimensional sur-

face, gas and soot radiation for turbulent buoyant flows in enclosures," First National

Conference on Numerical Methods in Heat Transfer, Sept. 1979

23. Yang, H.Q., "Laminar buoyant flow transition in three dimensional rectangular en-

closures," Dissertation for the Degree of PhD. Department of Aerospace and Me-

chanical Engineering, Universig of Notre Dame, 1986

84
24. Yang, HQ. and Yang, K.T., "Natural convection suppression in horizontal annuli

by azimuthal baffles’, Internap'enal Journal ef Heal; Mass Transfer, Vol. 31, No. 10,
pp.2123-2135, 1988.

25. Patankar, S.V. and Spalding, D.B., "A calculation procedure for heat, mass and
momentum transfer in three-dimensional parabolic flow", International Journal of Heat

mass transfer Vol. 15, pp.1787, 1972

 

26. Bhavnani, SH. and Bergles, A.E., "An experimental study of laminar natural con-
vection heat transfer from wavy surfaces", ASME Proceedings of the 1988 National

Heat Transfer Conference edited by Jacobs, H.R., HTD-Vol. 96, Vol. 2, 1988

 

27. Lee, S and Yovanovich, M.M., "Laminar natural convection from a vertical plate

with variation in surface heat flux", ASME weedings ef the 1288 Naa'onal Heat
Tpaasfer genferenee, edited by Jacobs, H.R., HI'D-Vol. 96, Vol. 2, 1988

28. Hwalek, J and Iyengar, V, "Natural convection mass transfer from vertical plate at

high mass transfer rates", ASME Prgeedings of the 1988 National Heat Transfer
Conference. edited by Jacobs, H.R., HTD-Vol. 96, Vol. 2, 1988

29. Kelleher, M. D. and Yang, K. T, "Heat transfer response of laminar free-convec—
tion boundary layer along a vertical heated plate to surface temperature oscillations,"

Journal of Applied Mathematics and Physics, Vol. 19, 1968

30. Show, H. J ., Chen, C. K. and Cleaver, J. W., "The effects of thermal sources on nat-

ural convection in an enclosure," International Journal of Heat and Fluid Flow3z,Vol.
9, No. 3, Sept. 1988.

 

 

"I7'71IITTII'IITIIIW