fjAH. INCT OF STAND & TECH

AlllOb

NBSIR 82-1660

MBS

PUBLICATIONS

COMPUTATION OF TWO-DIMENSIONAL

TIME-DEPENDENT NATURAL CONVECTION OF

COMPRESSIBLE FLUID IN A RECTANGULAR

ENCLOSURE

H. Yamashita

V. D. Arp

National Bureau of Standards

U.S. Department of Commerce

Boulder, Colorado 80303

March 1982

100

,U5o

32-1660

idol

c, 2

NBSIR 82-1660

COMPUTATION OF TWO-DIMENSIONAL

TIME-DEPENDENT NATURAL CONVECTION OF

COMPRESSIBLE FLUID IN A RECTANGULAR

ENCLOSURE * KQW

t.tw>a«T

1 7 1982

H. Yamashitaf

V. D. Arp

. az

iQ'Z'J

. >-

Thermophysical Properties Division

National Engineering Laboratory

National Bureau of Standards

U.S. Department of Commerce

Boulder, Colorado 80303

fGuest worker, Fukuoka University, Fukuoka, Japan

March 1982

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige. Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

NOMENCLATURE

Symbol

9

kPPo

AP

Q

R

T

AT

uv

XAX

y

Ay

specific heat, [J kg" 1 K _1

] arbitrary value, eq . (15

)

acceleration due to gravity, [ms-1

] thermal conductivity., [W nr 1 K

_1

] pressure, [Pa]

pressure at original point, = R % T 0

pressure difference, = P - P 0

total heat is added into the enclosure, [W]

gas constant, [J kg -1 K" 1

] temperature, [K]

temperature difference, = Tw

- T0

velocity component in x direction, [m s" 1

] velocity component in y direction [m s _1

] coordinate along the vertical wall, [m]

grid mesh size in the x direction, [m]

coordinate along the horizontal wall, [m]

grid mesh size in the y direction, [m]

Greek Letters

&

y

P

p

rAt

coefficient of volumetric expansion, [K _1

] specific heat ratio, Cp/Cy

absolute viscosity, [Pa s]

density, [kg m~ 3

] time, [s]

time step, [s]

Subscripts and Superscript

0

1jm

nrw

value at the wall except the heated wall

subscript denoting the i th grid point in

subscript denoting the j th grid point in

average value

superscript denoting the time at Tn restricted value

value at the heating surface

the

the

i i i

di recti on

di recti on

COMPUTATION OF TWO-DIMENSIONAL TIME-DEPENDENT NATURAL CONVECTION

OF COMPRESSIBLE FLUID IN A RECTANGULAR ENCLOSURE

H. Yamashita and V. D. Arp

Thermophysical Properties Division

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

Studies of natural convection processes generally assume

an incompressible fluid wherein the density is a function of

temperature only (the Boussinesq approximation). However,

local pressure gradients caused by rapid variations in the

heated wall temperature cannot be described within this approximation. These time-varying gradients cause fluid motions

which perturb the quasi-static natural convection process. In this study, we describe a numerical analysis procedure which

includes compressibility effects and allows computation of

transient fluid motions during onset of natural convection.

Details of the computational procedure and preliminary results

for one geometry are given.

Key words: compressible fluid motion; convection; finite

difference approximation; heat transfer; natural convection;

nonlinear convection; numerical integration; transient fluid

motion; transient heat transfer.

v

INTRODUCTION

For two-dimensional time-dependent laminar natural convection about a

/ p

O C Cl

heated surface, many numerical solutions are known in the literature. 5 ’’

These solutions use Boussinesq's approximation in which the thermophysical

properties are taken to be constant and the fluid is assumed incompressible,

except when considering the body force term in the equation of motion. These

assumptions lead to a numerical formulation in which the pressure terms are eliminated from the equation of motion and a stream function is defined which

satisfies the equation of continuity. Spiegel and Veroni's (4) have presented

the conditions under which the Boussinesq approximation is applicable for

thermal convection in compressible fluids; although the Boussinesq approximation

is valid for a number of natural convecting problems, it is of uncertain

accuracy for studies of supercritical fluid motion where the thermophysical

properties of these fluids are strongly dependent on temperature and pressure.

However, pressure changes generated by pulsed thermal input, which are significant in determining the fluid motion in early time periods, cannot be

described within such a formulation.

From this point of view, the authors have attempted to get computational

results using the program PDETWQ (7) for two-dimensional time-dependent natural

convection of a compressible fluid in a rectangular enclosure. However, stability problems arose in computational work, perhaps associated with implicit time

integration. Also, the cross differential term d^/{dxdy) in momentum

equations is neglected in the PEDETWO program, and it was not known whether this

could lead to a significant error.

In this paper, computational analysis is described on the laminar natural

convection heat transfer from an isothermal wall to compressible fluid within a rectangular enclosure, taking into account the variation of thermophysical properties and cross differential terms in momentum equations. Calculated results

for both pressure and buoyancy effects in early time periods for air are show".

1
MATHEMATICAL FORMULATION AND PHYSICAL DESCRIPTION

Consider the motion of a viscous fluid within a rectangular enclosure as shown In Fig. 1. The fluid Is Initially motionless with a uniform temperature

T0

. The enclosure walls are also at this temperature. The temperature

difference Tw-T0 is initiated at time 0 to Induce the flow within the

enclosure. The variations In all relevant physical properties are taken Into

account. However kinetic energy. Internal heat sources and Irreversible viscous

dissipation In the energy equation are not considered.

The governing equation for an compressible fluid with variable properties

will be as follows. (1)

continuity equation:

( 1 )

momentum equations:

( 2 )

where viscous forces on element are

(5)

( 6 )

( 7 )
energy equation:

4 *>

+ $. / 9P

T

+ ^ + v

,3P)

ay f

( 8 )

where

thermal expansion of the fluid.

The initial conditions are:

is the volumetric coefficient of

3P

3x

- - P9

(9)

( 10 )

( 11 )

and boundary conditions are:

t > 0 u

s v = o at all wal Is (12)

T * T0 at all walls except (13)

heated surface

I = Twat heated surface (14)

NUMERICAL SOLUTION OF THE EQUATION

Finite difference formulation

The left hand sides of eqs. (1), (2) and (3) can be manipulated by noting

that for arbitrary variable F as follows.

- & (PF) ^PUF) £<pvF) (15)

3
Using eq. (15), continuity eq. (1) can be written

9(p) + 9{(p) u)_ + = 0

3x 9 x

momentum eq. (2) is

3(pu) 9( (pu)u) 3{(pu)v)

9t 3x 3y

,4 / 3jj 3u 9

2 u\ 2 /3jj 3v\

1 ^3x 3x y 3x7 /” 3 ^3x 3yy

+ 3n 9u +

9

2 u +

3±j_3v 1 3

2 v 3y 3y P~9? 9y 9x I

y 3X3y

(17)

3(pv) . 9((pv)u) . 3{(pv)v>

3t 9x 3y

If. + i&i $1 +1 9

2 v\ (2 3^ 3u\ "

" 3y 3\3y 3y \3 3y 3 x/

+ lE . 9f_v + 3y 3u 1 9

2 u 3x 3x y

3x 2 3x 3y T

y 3x3y (18)

and energy eq. (8) is also written as follows.

( 19 )

4
The numerical scheme shown In Fig. 2 employs a 12x12 rectangular grid

system In the x and y direction with a total of 144 grid points and Ax, Ay are the grid size in both x and y direction respectively. This numerical method is

based on the simplest explicit finite difference approximations to the governing

differential equations which will be obtained at a finite number of grid points

having coordinate x s iAx and Y * jAY except next to the wall where the mesh

size is one-half the mesh within the cavity. All grid points are evaluated at

discrete times r

n

. The values of all physical properties at each grid point

should be thought of as average value over a small volume of fluid in Fig. 2.

The finite difference approximations on the derivatives of the arbitrary

variable F at the grid point (1, j) in advancing from time rn to the new level

Tn+1 = rn +Ar may be written as follows.

At (21)

The finite difference scheme for spatial derivatives uses central differences except next to the wall. The first order difference is

IE - ( F

i+i..1

3x 2Ax

F

niiLtil 9 ( 22 )

the second order finite difference is

2 cn

- ?Fn 1_F ( P i+U i.j

Sx2 Ax2

c

n

( 23 )

the cross finite difference is

5
d (24) 2 F dxdy

( F i4i, J+i - F

i-i .j-n

- Fn + F

n i+1 »j-l i-lJ-l)

4AxAy

The nonlinear convection terms cause the main difficulties to achieve a stable

numerical method. For some difference methods, the rate of heat removal may

differ from the rate of heat addition at steady state. Torrance (5) tested

several methods for differencing the convection term. In this calculation,

Torrance's V method (5) is employed as follows.

r

3(Fu)

3x

TTfl

= <

(F

n

.tf? . - F

n

, . I

j" .)

,J 1 ,J

•

fF

n Tf1 - F

n U

11 ) ^ i-’l . (TO . and U" ,

. < 0)

Av 1 »<J 1-1 »J —

If" TJ

11 _ p11 TF 1 ) 1+1 • J i-’U i-1.j

(u*? . U

1? ,

. < 0) • ' l.J i-l. J

1

fcfl Tin rll Tifl

Ax

(25)

(26)

(27)

where

« + l »

j

U

11m,n

Ki K °’ ^-i.o > 0)

(0? ,j > °* ^i-l ,j

<«»

1 (U

n + u

n

) 2

l Vl ,n um,n

;

When these approximations are Introduced into eqs. (16), (17), (18) and (19) we obtain

n+1

5 i,j

p" - M

~ (p<l.-i

1,J 2Ax

, (py)ri! - (pvvg-i

2Ay

,n+l

(28)

6
(P<1.

.

1 * J (pU)" , + At

* t J Ax

- <K.,i ^..i - ^.i-y ~Sy

(P? - p \ i+1,j i-1,,r n 2ax p iJ 9

+ <1.j ) J

[ 2ax 2Ax J

n f(U

n

, . - 2U

P

. + U

P

, .)1 + P

it jj

t+T.J

t Ax2 J -u^..i-^-i..i) ( vUi- v ij-in

4

(_ 2Ax 2Ay J

[ 2Ay 2Ay J

. f (u" . . - 211" • + U

n

• ,) I + U

n

. <

* 1.J+1 1.J T-J' 1 1 ’•J l W J

+ |

2Ay Tax J

+ Iy" . /

(V

i+l . j+l

" V

1-1,3>1

' V i^1,M *

(29) 3 i,o |

?AxZy JJ

'

(pv)?:] = (pV)" . + At[- ^pV

^i,j ^i,j ' ^ pV

^i-l ,j 1 * J L Ax

Ay

(p

n

- p

n

i ’ i’J*1 i J-l }

2A

y

7
*

L 2Ay 2 Ay

+ y

1 ? . (V

n

. Al - 2

V

Ay2

"

0

l 2 Ay 2Ax j

[ 2Ax 2Ax J

+ y

1 ? . ( V

" . ” 2V^ . + v

n

, .)

1 »J 1 +1.J y»J 1 Ax^

L 2ax 2Ay J 4»?i

~ U 1-l,3-Q

~ UH1,j-l

J 4AxAy

•»?!. 1-1 »]

(3o;

n+l _ Tn At Ln J

^ U«H,j

~ U i-1,j^ ^ V

i,j+1

~ V i,j-~P

] ij ‘

i , j C [‘ij I 2ax +

2Ay J

/’TTn T

n

_ TT1^ T

n

) (W*1 T

n _

sjft T

n

) - U

i »j i.i U i-l,J

1 j-1,j J - iv

i j'ij V

i J-l 1 1 , J-l ] Ax Ay

+ __j_ ;k n (<W 2 VT i-u) . (T i.^i - 2T

ij * ti.m

P

n .Cp

n

. I ’-A 4y^

*

(k

"+i..i

- kU.i> t^i.,1 -

,

< k

?..w - fc ".,i-i)

- T

2Ax 2Ax 2ay 2Ay

R"- Tn - c.i - (px)j!i..i)

,

- <<]

2Ax 2Ay

Bn

•• T

n

. + 11jl3 Lil

n r

n p iJ* Cp

i,j

L V/ L

fTTn

p

n _ Tj11 pn

) pn _ Tp p

n

) Ax Ay

- >"{ <1.1 ‘Cl ti

) 1J

l 2Ax

, ^

2 Ay J J J

8
where for an ideal gas

~ n i£ = . pi i ST —L»j- n

1 »

j

T?

BP _ D

n ST

' R - p i,j

JlJP - p.Tn 3p

R T

i ,j

3

= -1-

T

n

. T.J

c

. { . *-

•

}

1 Cp i .j J

Although the finite difference of convection terms in eqs. (29), (30) and (31) are shown for only one condition of eq. (25), appropriate finite difference of

convection terms should be used In actual computation work.

NUMERICAL PROCEDURE

The calculation proceeds by explicitly advancing p, n, v and T with

difference forms of eqs. (28), (29), (30) and (31). Also pressure P is

calculated explicitly from an equation of state using P and T. Fluid within

enclosure is initially at a uniform temperature To and at rest. Here, for

preliminary studies, we consider a rectangular enclosure of height Xmax (0.1m), width Ymax (0.05m) and vertical heated wall (0.04m) which is

located in the middle part of left side wall, as shown In figure 1.

During any one time-step, all values appearing in the right side of eqs.

(28), (29), (30) and (31) are treated as constants. In the first place, the new (pU)?+

l and (pV)?

+

! at all Interior grid points may be obtained from

1 » J • » J n+1

successive momentum eqs. (29) and (30). Then new density should be

calculated from continuity eq. (28) substituting (pU)J

+

l and (pV)"

+1

. Into eq.

9
n +i

(28). A new temperature Tv

j

is obtained from energy eq. (31) using the value

of just computed. Finally new velocities and V?

+

l are calculated

• »J i »j i

mathematically from (pU)^j and (pV)!Jj using the value of

p

n +

] as follows:

U

n+1 . (pU)

n+1

hi

> n+! i.J

(32)

New pressure P. • are calculated from the equation of state using new

1 * J temperatures and densities which are already computed at Interior grid points,

though not at the wall. Pressures at the wall are obtained by quadratic

extrapolations. This process is repeated in time, provided the time-step is

sufficiently small. The time-step Ax

f

has been restricted to 10-5$ or less

[5] in this computational work. This value corresponds to less than the time

interval for a sound wave to propagate across the mesh size y as follows

(33)

DISCUSSION OF RESULTS

Numerical calculations have been carried out for air within rectangular

enclosure (0.1 x 0.05 m2) in a time periods from 0 to 40 [ms]. The fluid

conditions and the imposed temperature differences correspond to a Grashof

number of approximately 8.8 x 10^. Figure 3-(l) ~ (10) shows the velocity

vectors at Intervals of 0.03 ms from 0.03 to 0.30 ms, where the dominant motion

Is normal to the heated wall (and gravity). The absence of a vector at a given

grid point means that the magnitude of the calculated velocity at that point was less than 5% of the maximum velocity at any of the grid points at that instant

of time. The disturbance-front separating the region of non-zero fluid

velocities from the region of essentially static fluid is seen to move away from

the heated wall with the velocity of sound. At t = 0.15 ms (figure 3-(5) ), the

disturbance front reaches the right side wall. In figures 3- (6 ) to (10), the

fluid motion becomes complicated by the sum of many phases and amplitudes of

10
motion with multiple reflection from all walls. No gravitational effects can be seen.Figure 4 — (1) to (10) shows the velocity vectors at later times, from

20.03 to 20.30 ms. The gravitational contributions to the fluid motion, which

causes assymmetry around the horizontal center line, is still quite small

compared to the motions induced by the initial expansion wave away from the

heated surface. In order to distinguish the growth of natural convection, the

velocity components Ug^

2 and Ug

2 are shown in Fig. 11 for increasing time.

Ug

2 and ug 2 are vertlca ^ components of velocity at grid points (5,2) and

(8,2) respectively. These grids points are symmetrical located about the

horizontal center line of the heated surface. The data are plotted at every

hundreth time step, which causes the apparent sawtooth character. Nevertheless,

the superposition of various amplitudes and phases, mentioned above, is clearly

evident. The dotted lines Indicate the values of-Ug

j2 , i.e. symmetrical values

of Ug

2 about bhe zero velocity. In the early time periods, less than about 5 [ms], the velocity components U

g ^ and -Ugs 2 are e Q ual within 1% or better.

After that time they begin to deviate, and this is a manifestation of buoyancy

force. The difference between U

g 2

and ~ U5,2» illustrated with shadow,

represents the growth of natural convection. The values of the shadow of

deduced from Fig. 11 are shown in Fig. 12. Also the difference of horizontal

components at both points are shown in same figure with the dotted curve. The

magnitudes of velocity component at upper point of (8,2) are larger than one at

lower point of (5,2) for either horizontal or vertical component, and this is

what should be expected on physical grounds In a natural convection heat

transfer. It is of interest to note that the natural convection flow near the

heated wall induced by the buoyancy force develop continuously and smoothly as shown in Fig. 12 and stream lines, which always close for incompressible fluid

'flow, will not do so in this case.

On the other hand, the temperature field in the vicinity of the heated wall

is essentially that of pure conduction for this range, and isotherms are practically symmetrical to the heated wall.

The pressure and velocity fluctuation at a near-mid point of the enclosure

are shown in Figs. 5 and 6 for time to 1.0 millisecond. The relation between

pressure and velocity fluctuation Is not distinct in Fig. 5 and 6 for this

physical model which has the heated wall at the middle part of the left side

wall as in Fig. 2. In order to reduce the influence of reflection at upper and

11
lower wall, the calculation was repeated with the heating surface extending over

the complete left wall of the cavity. The calculated pressure and horizontal

component of velocity are shown in Fig. 7 and 8. It is quite evident that the

frequency of pressure fluctuation is two times that of the velocity. The

average pressure in the enclosure increases as heat is added, and is

proportional to the heated area, as seen in Fig. 5 and Fig. 7:

The results of pressure and velocity fluctuation in later time periods from

30.5 to 32.0 [ms] are shown in Figs. 9 and 10 respectively. From these figures

we are not directly able to make clear the correlation between pressure and

velocity fluctuation. Therefore further study, using spectrum analysis (or

something similar) should be required to fully understand this phenomena where

the fluid motion is the sum of many phases and amplitudes of motion with

multiple reflection from the walls.

Numerical procedures for solutions of heat transfer equations In the

time-dependent domain may fall Into two categories, explicit and Implicit.

These two types of difference equations have previously been studied in which

explicit difference equations are simple to solve, but which require a large

number of time steps of limited size, and implicit difference equations do not

limit the time step but which do require iteration at each time step in the solution. Therefore, explicit procedures are convenient under conditions where a sufficiently large time step, consistent with computational stability, can be

used. In order to examine the accuracy of this computational results, the

smaller time step Ar of 10-4 [ms] which is one tenth of normal time step

have been used for early time periods from 0 to 0.3 [ms]. The agreement with

the solutions of vy^g for this calculation of small time step and the

prior one Is better than 0.1%.

The numerical calculations have been performed on a large digital computer.

The execution time at each time step is approximately 33 CP millisecond.

12
ACKNOWLEDGEMENT

This work was performed during the time that H. Yamashita was a guest

worker in the Thermophyscial Properties Division of NBS. Cooperative studies

using the method developed here will continue after Prof. Yamashita returns to Fukuoka University in Japan.

13
Figure 1. Physical model and coordinates.

14
» 'i > i

Figure 2. Schematic diagram of the numerical method.

15
X

(MM)

10.0

ao.o

ss.0

oo.o

so.o

®o.§

m,®

oo.o

n.o

me

\ / * / t

\ t / S X X / t t f

+ l V \ \

\

^ r X / /

/ \ \ X ^ x x X \ |

/NX X x X \ \

X \

0.0

“T" 10.0

•r

1 i 30.0 90.0 y (mm)

t

0.0 so.o

Figure 3-(7). Velocity field at t-0.21 ms. Vmax«0.25 mm/s.

16
X

(MM)

Figure 3— ( 6 ) . Velocity field at t*0.18 ms. • 427 mo s

.

17
X

(MM)

10.0

30.0

SS.O

«

5.0

90.6

96-0

76

.®

60.0

<

30.0

100.0

4 i //////

i /

t / / / / / / / * /

/// / / s s'^

x \ v ^ ^

I \ \ \ \ \ ^ ' s '

r \ v v

1 4 N v \ V V ^ '

o

o.

0.0 T

10.0 T30.0 30.0 T

40.0 60.0 r ( mm)

Figure 3-(5). Velocity field at t«0.15 ms. Vmax-0.428 mm/s.

18
CUM) X

Figure 3-(4). Velocity field at t**0.12 ms. vmax“0*376 mm/s.

19
Figure 3-(3). Velocity field at t-0.09 ms. Vmax-0.566 mm/s.

20
X

(HI)

Figure 3— ( 2 )

Velocity field at t-0.06 ms. Vmax-0.368 mm/ s

.

21
(UU) X

Figure 3— Cl). Velocity field at t-0.03 ms. vmax"0*521 mm/s.

22
X

(MM)

o

Figure 3-(8). Velocity field at t«0.24 ms. Vmax-0.498 mm/s.

23
Figure 3-(9) . Velocity field at t-0.27 ms. V^x-0.214 mm/s.

24
X

(MM)

Figure 3-(10). Velocity field at t-0.30 ms. Vmax-0.544 mo's.

25
X

(MM)

Figure 4-(l). Velocity field at t-20.03 ms. Vmax-0.343 mm/s.

26
Figure 4-(2). Velocity field at t-20.06 ms. Vmax**0.380 s

.

27
Figure 4-(3). Velocity field at t-20.09 ms. mm/s.

28
X

(MM)

30.0

40.0

50.0

00.0

70.0

00.0

00.0

100.0

Figure 4-(4). Velocity field at t-20.12 ms. Vmax-0.567 mm's.

29
o

Figure 4-(5). Velocity field at t-20.15 ms. Vn^x-0.339 mm/s.

30
X

(Ml)

10.0

30.0

30.0

40.0

80.0

00.0

T0.0

00.0

00.0

100.0

4 4 4 *

¥ \ \ ^

0.0 —I

10.0 —

r

ao.o —i

30.0

1 40.0 so.o Y (MM)

Figure 4-(6). Velocity field at t-20.18 ms. Vmax-0.588 mo/s.

31
X

(MM)

Figure 4-(7). Velocity field at t-20.21 ms. Vmax“0 * 2 ^ 0 mm/ 8 *

32
X

(MM)

10.0

20.0

30.0

00.0

50.0

#0.0

70.0

W.O

90.0

100.0

4

< f i k k

k

k

— X \ N \ ' * ' ' '

\

'i ^ \ \ \ i r r y

^ -s ^ 1 » v 1

H T

0.0 —f

10.0 t20.0 —I 30.0 —i

40.0 50.0 X (MM)

Figure 4-( 8 ) . Velocity field at t-20.24 ms. Vmax

-0.527 mm s

33
NWS X

Figure 4-(9). Velocity field at t-20.27 ms. V^-0.325 mm/s.

34
y

(fT!)

JO.O

30.0

30.0

iL'.O

W.O

00.0

70.0

80.0

90.0

100.0

2.0 10.0 —I ao.o 3o.o X (MM)

-i

40.0 50.0

Figure 4-(10). Velocity field at t-20.30 ms. Vmax-0.531 mm.'s.

35
References

[1] Bird, R.B., W.E. Stewart and E.N. Lightfoot, Transport Phenomena, John

Wiley, New York, (1960), pp.322.

[2] Kublbeck, K., G.P. Merker and J. Straut, Advanced numerical computation of

two-dimensional time-dependent free convection in cavities, Int. J. Heat

Mass Transfer. 23, (1979) 203.

[3] Mellums, J.D. and S.W. Churchill, Transient and steady state, free and

natural convection, numerical solutions, A.I.Ch.E. Journal 8-5, (1962)

690.

[4] Spiegel, E.A. and G. Veronis, On the Boussinesq application for a compressible fluid, Astrophysical J., 131, (1960) 442.

[5] Torrance, K.E., Comparison of finite-difference computations of natural

convection, J. of Res. of the NBS, 72 3-4, (1968) 281.

[6] Vasseur, P. and L. Robillard, Transient natural convection heat transfer in

a mass of water cooled through 4°C, Int. J. Heat Mass Transfer, 23, (1980)

1195.

[7] Melgaard, D., R. Sincovec, General software for two-dimensional partial

differential equations. Rep CS76-21, Dept. Computer Science, Kansas State

Univ. (1976).

36
Line

43

35

127

127

35

130

131

135

136

107

108

78

78

97

98

78

198

198

78

180

193

181

194

73

Appendix B

List of variables [Main Program]

These are founded on computer program list of 1.

Program Symbol Definition

A =g/(RT), to use initial density distribution

ADATE Date of computation

AFPX Dummy argument, see sub. AVEVELO

AFPY

t. .•

ATIME Time of computation

AUPX Average velocity u(i) = 0.5[u(i+l,j) + u(i,j)]

DUPY Average velocity u(j) = 0.5[u(i,j+l) + u(i,j)]

AVPX Average velocity v(i) * 0.5[v(i+l,j) + v(i,j)]

AVPY Average velocity v(j) = 0.5[v(i,j+l) +v(i,j)]

CODX dk/dx

CODY dk/dy

CONE Thermal conductivity, k

CP Specific heat, Cp

CPDX dCp/dx

CPDY dCp/dy

DFDX Dummy argument of d/dx

DFDXX d

2 /dx 2

DFDXY " d

2 /dxdy

DFDY " d/dy

DFDYY " d

2/dy 2

DPDR dP/dP

DPUDX * d(Pu)/dx in energy equation

DPVDY - d(Pv)/dx

DRDT - dP/dT

DTDX - dT/dx

DTDXX - d

2 T/dx2

DTDY - dT/dy

DTDYY - d

2 T/dy 2

DTIME -= At time step

37
Line No . Program Symbol Definition

197

| DTUDX 8(Tu)/8x difference of convection term

197

| DTVDY - 8(Tv)/dx

32

| DTW s Tw-T0 , temperature difference

| DUDX c du/ dx

| DUDXX m 8 2 u/dx2

| DUDXY sr 8 2 u/dx8y

| DUDY s 8u/dy

| DUDYY ss 8 2 u/8y 2

| DURDX = 8(Pu)/dx

| DURDY ac 8(Pu)/8y

| DVDX = 8v/dx

| DVDXX s d

2 v/dx2

| DVDXY s 8 2 v/ dxdy

| DVDY ss 8v/8y

| DVDYY = 8 2 v/8y 2

| DVRDX * 8 (Pv)/dx

| DVRDX s 8( Pv)/8y

1 DX - dx

| DY - dy

1 G - g, acceleration of gravity

| I i th grid point in x direction

| ICO - 0, printout at first time step

74,225 | ICOUNT counter on the number of time-step for printout

232

1 « integer of changing of I th order for printout

59

1 IS1 The lowest grid point of heated surface

59

I IS2 The highest grid point of heated surface

225

| I WRITE Controller integer value of printout

1 J j th grid point in y direction

1 K K has 1,2 and 3. K=3 : newest one; K=1 : old one

1 M M th grid point is correspond to Xmav

| Ml - M-l

| M2 ST M-2

1 N N th grid point corresponds to ymax

N1 - N-l

| N2 m N-2

38
Line No. Program Symbol Definition

PRDX - PRDY - PRESS pressure, P

R density, P 149 RU * Pu

RUDX d(Pu)/dx

157 RV “ PV

RVDY d(pv)/dy

T temperature

32 TB static temperature or wall temp, except heated wall

269 TMAX limitation of computing time

61 TW * Tw , temperature at heated wall

73 TYME * t + At , increment time in computation

199 T1 some term in energy equation

200 T2

• 0

201 T3 •«

202 T4

•t

204 T5

••

205 T6 M

206 T7

M

U velocity component in x direction

UDV is not used, only dimension

UVMG is not used, only dimension

143 U1 some term in mementum equation

144 U2

•t

145 U3

M

147 U4

”

V velocity component in y direction

145 VIDX = dp/dx

147 VIDY = d^/dy

vise * P, viscosity

151 VI some term in momentum equation

152 V2

"

153 V3

M

155 V4

M

45 X * X

30 XMAX “ xmax

31 YMAX ymax

39
[Subroutine]

Line

***

***

***

***

***

***

No. Program Symbol Definition

PROPER Subprogram Name [calculation for properties]

COND thermal conductivity " k

CP specific heat, « C

p

DPDR - dP/dP

DPDT - dP/dT

DRDT - dP/d

M3 - M-3

N3 - N-3

PRESS pressure, P

R density, « P T temperature, * T

vise absolute viscosity, - m

FIRSTDF Subprogram Name [difference for first order]

DFDX F is dummy value, * dF/dx

DFDY “ dF/dy

DX mesh size in x direction

DY mesh size in y direction

F dummy argument value

SECONDDF

DFDXX

Subprogram Name [difference for second order]

- d^F/tfx2

DFDYY - d2 F/(9y

DXX - (Ax) 2

DYY - (AY) 2

QDXX * 1/4 (Ax) 2 QDYY - 1/4 (AY) 2

CROSSDF Subprogram Name [cross difference]

DXY * Ax * Ay

CONVDF Subprogram Name [difference for convection term]

AUPX see page 1 of this appendix

AUPY

AVPX

AVPY

DUFDX * d(uF)/dx F is dummy argument

DVFDY d(vF)/dy

AVEVELO Subprogram Name [calculation of average velocity]

AFPX - AUPX or AVPX

AFPY • AUPY or AVPY

40
NBS-114A (REV. 2 -80

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO. NBSIR 82-1660

2. Performing Organ. Report No. 3. Publication Date

March 1982

4. TITLE AND SUBTITLE

Computation of Two-Dimensional Time-Dependent Natural Convection of Compressible

Fluid in a Rectangular Enclosure

5. AUTHOR(S)

H. Yamashita and V. D. Arp

6. PERFORMING ORGANIZATION (If joint or other than NBS, see in struction s)

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State. ZIP)

10. SUPPLEMENTARY NOTES

_ Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a si gn i fi cant bibliography or literature survey, mention it here)

Studies of natural convection processes generally assume an incompressible fluid

wherein the density is a function of temperature only (the Boussinesq approximation).

However, local pressure gradients caused by rapid variations in the heated wall

temperature cannot be described within this approximation. These time-varying

gradients cause fluid motions which perturb the quasi-static natural convection process. In this study, we describe a numerical analysis procedure which includes

compressibility effects and allows computation of transient fluid motions du-ring

onset of natural convection. Details of the computational procedure and preliminary

results for one geometry are given.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words r\ se - . c

-

: compressible fluid motion; convection; finite difference approximation; heat trar. s: c:-, natural convection; non-linear convection; numerical integration; transient :.ccc. motion; transient heat transfer.

13. AVAILABILITY 14. NO. OF

[XXUnlimited

| | For Official Distribution. Do Not Release to NTIS

Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

PRINTED PAGES

45

15. Price 20402. XX Order From National Technical Information Service (NTIS), Springfield, VA. 22161 $7.50

USCOMM-DC

US. GOVERNMENT PRINTING OFFICE. 198C-5'# .'.V • REG CNN - ?