NBSIR 83-2804

Thermal Flanking Loss Calculations

for the National Bureau of

Standards Calibrated Hot Box

R. J. Onega

Department of Mechanical Engineering

Virginia Polytechnic Institute and State University

Blacksburg, VA 24061

and

P. J. Burns

Department of Mechanical Engineering

Colorado State University

Fort Collins, CO 80523

February 1 985

I

Prepared for

QC

100

. U56

Department of Energy

Oak Ridge National Laboratory

inergy Division

)ak Ridge, TN 37830

83-2804

1985
NATIONAL INSTITUTE OF STANDARDS &

TECHNOLOGY

Research Information Center

Gaithersburg, MD 20899
NATIONAL BUREAU OF STANDARDS

LIBRARY V . A l At

NBSIR 83-2804

t Hr f / Jf / A A

THERMAL FLANKING LOSS CALCULATIONS

FOR THE NATIONAL BUREAU OF

STANDARDS CALIBRATED HOT BOX

R. J. Onega

Department of Mechanical Engineering

Virginia Polytechnic Institute and State University

Blacksburg, VA 24061

and

P. J. Burns

Department of Mechanical Engineering

Colorado State University

Fort Collins, CO 80523

February 1985

Prepared for Department of Energy

Oak Ridge National Laboratory

Energy Division

Oak Ridge, TN 37830

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

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
FORWORD

This report is one of a series documenting NBS research and analysis efforts

in support of the Department of Energy/Oak Ridge National Laboratory/National

Bureau of Standards' Building Thermal Envelope Systems and Insulation Materials

program. The work reported in this document was performed under the Calibrated

Hot Box project (Task 4) supported by DOE/NBS Interagency Agreement No. DE-AI05-

780R061 13

ii
EXECUTIVE SUMMARY

Two frame coefficients are defined FCMCA and FCHXS which describe the fraction

of the heat lost through the frame. FCMCA is related to the actual heat flux

that goes through the wall specimen and FCHXS is related to the heat flux that

goes through the specimen if it was one-dimensional.

The code FLANK (see Appendix) was written to obtain these coefficients, the temperature distribution and the heat fluxes across various surfaces in both a steady state or transient mode. The code is a two dimensional, finite difference

code. The number of nodes that can be used is variable and is limited only by

computer memory and the amount of time and money available for computer runs.

If one considers the steady-state conditions for the metering chamber, the

heat input to the metering chamber Qjfj is related to the heat out by

QlN = ^specimen Qfrarae QMC walls where Qframe -®- s *-he heat that goes out

through the frame and Qmc walls is the heat lost through the metering chamber

walls. It is easily seen that

Qfrarae = QiN (1 “ specimen

) _ wallg

Qin

= Qin FCHXS - Qmc walls

if the heat through the specimen is assumed to be ks As (Tecj-Tmc].)AYs where ks and As are the thermal conductivity and area of the specimen, AYS is the thickness of the specimen and Tg^j and T^ci are the surface temperatures at the

center of the climatic and metering chamber sides of the wall. FLANK provides

FCHXS (and FCMCA) for a given set of conditions.

Typical applications for a polystyrene wall indicate that FCHXS will be between

5 and 15 percent in steady-state applications. If the temperature of the climatic chamber is equal to the temperature of the metering chamber, it can be seen that the frame coefficients will be very large (-1). In the transient

case, the polyurethane will take about 40 or 50 hours to reach steady state.

The frame coefficients likewise will achieve their steady state values in

about 50 hrs or so. Also, the frame coefficients can be very large (=50 percent) in the transient situation due to the fact that the polyurethane frame

itself is being heated.

iii
TABLE OF CONTENTS

Page

EXECUTIVE SUMMARY ill

1. INTRODUCTION 1

1.1 INTRODUCTION 1 1.2 BACKGROUND 1 1.3 THE NBS-CHB SPECI FI CATIONS 2 1.4 SUMMARY OF TESTS TO BE RUN 2

1.5 CHB SENSORS 2 1.6 SUMMARY OF THERMOPHYSICAL DATA 6

1.7 FLANKING LOSS DEFINITION (FRAME COEFFICIENTS) 6

2. FLANKING LOSS MODEL 9

2.1 THE MATHEMATICAL MODEL 9

2.1.1

Heat Transfer Coefficients 12

2.2 STEADY- STATE ANALYSIS 14

2.3 TRANSIENT ANALYSIS 14

3. RESULTS 20

3.1 CODE VERIFICATION AND DISCUSSION 20

3.2 STEADY-STATE RESULTS AND ANALYSIS 20

3.3 TRANSIENT ANALYSIS AND RESULTS 22

4. SUMMARY 34

5. REFERENCES 35

iv
LIST OF TABLES

Page

Table 1.1 Specifications of the Walls of the CHB 3

Table 1.2 Test Specimen and Specimen Frame Specification 3

Table 1.3 Operating Conditions that may be Imposed on a Specimen

in the NBS-CHB 4

Table 1.4 Climatic Conditions to be Simulated 5

Table 1.5 Sensors in the NBS-CHB 6

Table 1.6 Nominal Thermophysical Parameters for CHB Construction

Materials and Some Specimen Materials 7

v
LIST OF FIGURES

Page

Figure 2.1 Geometry (Region A) of the CHB Frame, wall specimen,

gaskets and portions of the MC and EC walls 16

Figure 2.2 The various regions and the number of nodes in each region

for the CHB as coded in FLANK 17

Figure 2.3 Nodal structure of the frame, specimen, gaskets, and

portions of the MC and EC walls 18

Figure 2.4 Heat transfer coefficient due to forced convection that

applies to the metering and environmental chamber

specimen and frame surfaces 19

Figure 3.1 The node structure used in the calculations presented 24

Figure 3.2 The influence of node structure on the accuracy of heat

fluxes across various surfaces 25

Figure 3.3 Heat flux numbering system 26

Figure 3.4 Sensitivity of the temperature in the FRP to the time

step 27

Figure 3.5 The influence of a variation of the climatic chamber

temperatures on the frame coefficient 28

Figure 3.6 The variation of the frame coefficients with the thermal

conductivity of the wall specimen 29

Figure 3.7 The frame coefficient as a function of position from the

metering chamber face of the 0.3048 m wide frame 30

Figure 3.8 The variation of the frame coefficient with the wall

thickness 31

Figure 3.9 The transient response of the polyurethane frame and the

specimen to abrupt changes in boundary conditions 32

Figure 3.10 The transient response of the frame coefficients 33

vi
1

INTRODUCTION

1.1 INTRODUCTION

A computer code called FLANK (see Appendix) was written to calculate the flanking

loss for the NBS Calibrated Hot-Box (CHB). This code is a two-dimensional

finite difference dynamic thermal simulation of the test frame and specimen

when subjected to either steady-state or time-varying boundary conditions.

The code is specifically designed to do the calculations for the NBS-CHB but

the code is general enough to be used with any CHB.

This report presents some background information of the NBS-CHB since the code

was written so that those quantities which are measurable or controllable are either inputs or outputs of the code. Following this is a description of the

mathematical model, a description of the code and results. A users manual

for FLANK is also available.

1.2 BACKGROUND

The purpose of the NBS Calibrated Hot-Box (CHB) is primarily to measure the

heat flow through large components of a building envelope, such as a wall. A

test specimen is constructed in a frame which is then placed between two chambers. One chamber is called the metering chamber (MC) and is designed to simulate inside conditions and the remaining chamber is called the environmental

chamber (EC) which simulates outside conditions. Wider temperature variations

are possible in the EC than the MC.

The following purposes can be enumerated for the NBS-CHB:

1. Measure the heat, air and moisture transport through large homogeneous

or composite wall system under steady-state and dynamic conditions.

2. Provide the methodology in CHB measurements that can be used for

accrediting various testing laboratories.

3. Generate thermal transport properties for full scale wall systems

which may include windows, doors, etc.

4. Provide a mechanism for traceability to NBS in the measurement of

thermal properties of wall specimens through a concept such as the

Standard Reference Material (SRM) program.

5. Verify models of heat, moisture, air transport and flanking loss

through walls by comparing model results to experimental data.

A flanking loss correction must be made to provide a comparison of theory and

experiment. The flanking loss is a function of the wall thickness, material,

placement of the wall in the frame, the test frame construction and the temperature in the MC, EC and the environment. Since air temperatures, heat transfer

coefficients, etc., can be experimentally imposed, the model must be able to

accommodate these variables as inputs.

1
1.3

THE NBS-CHB SPECIFICATIONS

The CHB is to test walls with a thermal resistance ranging from an R of

0.352 m2*K/W to 8.81 m^»K/W (R-2 to R-50 hr*ft^ °F/Btu) under a variety

of temperature, moisture and air infiltration conditions. The MC Walls are made of polyurethane board about 0.38 m (15 in) thick surrounded on all outside

surfaces by a water jacket designed to minimize the temperature difference

between the outside and inside of the walls. The interior and exterior surfaces of the polyurethane MC walls are covered with about 0.003175 m ( —1 / 8 in)

fiberglass reinforced plastic (FRP) of type 45 fiberglass polyester 45 percent

glass matt base with a fire-retardant filler. Specifics of the wall construction construction are presented in table 1.1.

The EC walls consist of 0.454 m (17-3/4 in) polyurethane covered with FRP. The

FRP is finished with two coats of epoxy-base paint. The R-value for the EC

walls are about 19.02 m^»K/W (R - 108 h ft^*°F/Btu).

The wall specimen is mounted in a polyurethane frame. Details of the frame and

specimen are presented in table 1.2. The wall specimen itself is 3.048 m high

(10 ft), 4.572 m wide (15 ft) and can vary in thickness up to a maximum of

0.6096 m (2 ft). Neoprene seals are used to sandwich the test frame between

the MC and the EC.

Table 1.3 provides a list of operating conditions available in both the MC and

the EC.

1.4

SUMMARY OF TESTS TO BE RUN

Various test specimens will be put into the CHB and subjected to various

climatic conditions. These walls may be homogeneous or composites and may

contain windows and doors. Table 1.4 is taken from reference [1]. Three

different types of parameters can be varied: a) thermal quantities - air

temperatures in both chambers; b) air infiltration parameters - static pressure

difference across the wall; c) moisture parameters - relative humidity in both

chambers. These parameters determine the heat flow through a wall with given

thermal characteristics [2]. Both static and dynamic tests are to be run.

The control systems are to maintain constant dry-bulb temperatures to within

+0.11 K of their setpoints in both chambers. A static pressure differential

across the test specimen is to be maintained to within +2.5 Pa of its set

point. The air flow rates can be maintained constant and will generally be

directed in the direction that natural convection would take place.

1.5

CHB SENSORS

The results of computer models must be compared with experimental data to

determine their accuracy and range of validity. It is, therefore, important

to know what sensors are used, their accuracy and precisely what quantity is

being measured. Table 1.5 provides a list of these sensors. The sensors in

the specimen itself will depend upon the tests and the specimen.

2
Table 1.1 Specifications of the Walls of the CHB

QUANTITY METERING CHAMBER ENVIRONMENTAL CHAMBER

Outside height 4.3 m 4.3 m

Outside length 6.7 m 6.7 m

Outside width 1.50 m 1.57 m

Inside width 1.117 m 1.117 m

Inside height 3.048 m 3.048 m

Inside length 4.572 m 4.572 m

Wall thickness 0.38 m 0.4572 m

Wall construction polyurethane polyurethane

Inside volume 15.6 m^ 15.6 mJ Thermal resistance 19.0 m^ »K/w

Table 1.2 Test Specimen and Specimen Frame Specification

QUANTITY VALUE

Construction Rigid polyurethane foam covered

with a layer 0.003175 m of

fiberglass

R-Value Reinforced plastic (FRP)

0.352 to 8.81 m^'K/w

Inside height 3.048 ra

Inside width 4.572 ra

Frame width 2 frames (0.3048 m and 0.6096 m)

Frame thickness 0.4572 m

Specimen area 13.935 m

2

Specimen weight/area 700 kg/m2

Air leakage rate (max.) 255 m3 /hr

Moisture transfer rate (max.) 3.89 kg/s

Static pressure diff. (max.) 125 Pa

3
Table 1.3 Operating Conditions That May be Imposed on a Specimen in the NBS-CHB

PARAMETER METERING CHAMBER ENVIRONMENTAL CHAMBER

Dry bulb temperature range 283 to 338 K 233 to 338 K

Dew point temperature range 277.4 at 283 K D.B.

to 293 at 297 K D.B.

and above

230 at 233 K D.B

to 293 K D.B.

and above

Heat input rate (max.) 7.3 kW 12.3 kW

Heat removal rate (max.) 6.4 kW 12.3 kW

Velocity of air curtain 0.25 to 0.75 m/s 0.41 to 4.1 m/s

Maximum diurnal

temperature amplitude

17°C 56°C

Rate of temperature rise

at 27°C (min.)

3.06 x 10"3 °C/s 3.06 x 10~3 °C/s

Rate of temperature

decrease at -7°C

2.72 x 10"3 °C/s 2.22 x 10"3 °C/s

4
Table 1.4 Climatic Conditions to be Simulated

Transmission Phenomena

SEASON Thermal Air

Direction3 S.S or Dynb Direction3 S.S or D;

Winter Out S.S Nonec Winter Out Dyn None

Summer In S.S None

Summer In Dyn None

Winter Out S.S Recirce S.S

Summer In S.S Recirc S.S

Winter Out S.S Recirc S.S

Winter Out Dyn Recirc Dyn

Summer In Dyn Recirc Dyn

Winter Out S.S In S.S

Winter Out S.S Out S.S

Summer In S.S In S.S

Summer In S.S Out S.S

Winter Out Dyn In S.S

Summer In Dyn In S.S

Winter Out Dyn Out S.S

Summer In Dyn Out S.S

Winter Out S.S Out S.S

Winter Out Dyn Out S.S

a. "In" and "Out” refer to the direction of heat, air, or moisture transfer

as viewed from the Metering Chamber.

b. "S.S" stands for steady-state. "Dyn" means dynamic or variable.

c. "None" means no static pressure difference across the specimen.

d. "None" means that no vapor pressure difference will be maintained across

the specimen.

e. "Recirc" means two-way flow of air can take place through some wall

specimens due to stack effect.

5
Table 1.5 Sensors in the NBS-CHB

QUANTITY

1. Temperature

a) specimen surface

temperature

b) air temperature

i) dry bulb

ii) dew point

temperature

2. Heat flux

3.

Pressure

4. Power input

a) heat to metering

chamber

5. Air velocity

METERING

(ACCURACY)

CHAMBER

RTD

30 heat-flow meters

flat, 2-1/4", installed

between polyurethane

and FRP in interior

of chamber

(+0.5 percent of

reading)

ENVIRONMENTAL

(ACCURACY)

CHAMBER

RTD

same as M.C.

40, 24 gage Cu-Co 40, 24 gage type T

thermocouples - type T thermocouple

(+ 0.5 K) (+ 0.5 K)

1.6 SUMMARY OF THERMOPHYSICAL DATA

Steady-state and dynamic modeling require knowledge of thermophysical

parameters. Table 1.6 is a list of the relevent parameters for some of the

materials from which the CHB is constructed. In addition to therraophysical

parameters, air infiltration and moisture parameters are also required.

1.7 FLANKING LOSS DEFINITION (Frame Coefficient)

The flanking loss is loosely defined as the amount of heat that leaves the MC

which does not go through the wall and ends up in the EC. Lavine et al. [4]

define the flanking loss as "heat flow also occurs between the metering chamber and the climatic chamber through the specimen frame. This heat flow is

called the flanking loss."

According to the definition of Lavine, et al., the flanking loss is only that

heat which flows from the MC, through the frame and then into the EC.

6
Table 1 .6 Nominal Thermophysical Parameters for CHB Construction Materials

and Some Specimen Materials

This definition result in several questions which are listed below:

i) How is the heat that flows through the edge of the test wall accounted

for in this definition?

ii) How is the heat that flows from the MC, through the frame and out into

the surroundings accounted for in this definition?

iii) How does one uniquely calculate flanking loss as defined above?

Because of these difficulties, the following definitions of a quantity termed

"frame coefficient" are preferred. Both of these definitions depend upon the

deviation from one-dimensional heat flow for their basis.

The total heat flow out of the MC (for steady-state) is observed from the

difference between the amount produced from the electrical heater and the

amount removed by the MC cooling coils, i.e.,

Qin = QEH “ Qcc • (!•!)

The heat going through the center point of the specimen is one-dimensional but

that which goes through the specimen near the frame has a component in the z

7
direction as well as the y direction (see figure 2.1). The heat flux in the y

direction is

( 1 . 2 )

where the subscript "s" refers to the specimen. This is the amount of heat

that would flow through the wall if it were very large so that edge effects

(frame coefficients) could be ignored.

The heat flux is the amount of heat that crosses a unit surface area per unit

time. Using equation (1.2), if the expression is evaluated at the geometrical

center of the wall in the z direction and the front surface (y = 0), a definition of the frame coefficient

If the average heat flow through the front surface of the wall is used,

FCHXS = [Qln - Qsy (y*0, z=H/2)]/Qin (1.3)

H

FCMCA = 1 w

‘/o VTs (z >y°0) ~ TMC ]dz

Qin (1.4)

where w is the width of the wall. Other definitions can also be given for the

flanking loss or frame coefficients as well.
2. FLANKING LOSS MODEL

2.1 THE MATHEMATICAL MODEL

FLANK calculates the temperature distribution, the heat fluxes across various

boundaries and the frame coefficients. The geometry is shown in figure 2.1.

Since both transient and steady-state runs are to be made in the CHB, the code

must be able to be run in either mode.

All the therraophysical properties such as thermal conductivity, specific heats,

etc., are assumed to be independent of temperature. The heat balance equation

for the region, "A", shown in figure 2.1 is

pc il = JL (kil) + JL (kil) , (2.i)

P 9t 9y 9y 3z 9z

where T is the temperature at any point (y, z) at time, t, p is the density, Cp

is the heat capacity and k is the thermal conductivity. The boundary conditions

are of two types. Since the CHB has coils in the MC wall to maintain a constant

temperature, the constant temperature boundary condition

T(t) - T0 on S2 (2.2)

where S 2 is a line where the coils exist for heat rejection if needed.

The other type of boundary condition is a specification of the heat transfer

coefficient and ambient temperature along the line Sj

k|I + h[T <t) - 1^(01 = 0 on Sly (2.3)

S iy

and

9z

“iz

+ h[T

Siz

(t) - T^Ct)] =0 on Siz (2.4)

Here S-^y and S^ z are boundary lines along the y and z axis respectively and h is

the heat transfer coeficient.

The region designated A consisting of the frame, specimens, and a portion of

the MC and EC walls is shown in figure 2.2. This was broken into nodes as shown

on figure 2.3. A finite difference formulation of the heat conduction equation

was made in order to calculate the temperature at the node, heat fluxes across

lines and the frame coefficient. Heat balances at the various nodes were put

into a matrix equation formulation which resulted in the following set:

(P+1) (P) , (P+1) (P) (P+1) (P)

C(T - T )/At + a K T + ( 1“«) I 1

= aB +(l-a)B (2.5)

9
where £ and K are (n x n) matrices, (n) is the number of nodes, and jr( p ) is the

temperature at the various nodes at time step "P" . The vector arises

because of boundary conditions and a is the time integration parameter such that

0 < ot < 1 . If a * 1, the solution is purely implicit (backward difference formulation), if o 0, it is an explicit solution scheme and if a*l/2, the solution

method is a Crank-Nicholson scheme. If a-1 is used, no numerical instability

problems will arise due to time-step sizes. The £ matrix is the thermal capacity

matrix and K is the thermal conductivity matrix.

Equation (2.5) can be rearranged to yield

(P) (P) (P+1) (P)

C < DD + aK (DT) = a B + (l-a)B + K T (2.6) “At ” “ “

where

DT( p ) = tP+1 - T

p (2.7)

Equation (2.6) is the form that is actually used in FLANK. The £ is a diagonal

matrix where the i

t diagonal element is given by Ay^ Az^ and Ay^ and

Az^ must be evaluated for both center as well as corner nodes properly. The

K matrix consists of elements for the ij node as

Ay. Az.

Kj j = k. —1 + k, —J- (2.8) 1 Azj J Ay±

and, of course, k. Ay, and Az must be properly evaluated for each of the

different types of node.

The matrix equation (2.6) is solved using three IMSL subroutines [5], LUDATF

which decomposes JC and K using the Crout Algorithm, LUELMF which carries out

the elimination part of the solution, and LUREFF which refines the solution

to achieve more accuracy.

The Crout reduction [6] or decomposition of a matrix K works basically as follows.

The matrix K is written as the product of two matrices, i.e.,

K = L U (2.9)

where L, is lower triangular and U is unit upper triangular, so the ij^ element

can be expressed as

min (i , j

)

ki-j = Z ujm (i,j = 1,2,3.. .n) . (2.10)

Jl=l

Since ]1 is a unit upper matrix, the diagonal elements are forced to unity so u^^ =1. Therefore,

kij = x il U 11 = x il (i = 1,2,.. .n) ( 2 . 11 )

10
Also

k ij = X 11 u lj

or

u lj = k lj/Ml for J = 2,3,4,.. .n . (2.12)

The first column of ^ is determined by (2.11) and the first row of u is

determined by (2.12).

After the first m-1 columns of L and m-1 rows of £ are established, the relation

m-1

k im = Mm ^ M& u£m 1 = m » m"*‘l » • • • > n

Xr — 1

is used, or solving for X^m ,

m-1

Mm = kim — ^ M& u£m (2.13)

1=1

The relationship which permits the computation of uroj

is

1 m-1

umj =

X f kmj “ ^ ^mJt u£j (2.14)

nTM SL=1

Notice that the mtk row of IJ and the mt ^1 column of L_ are determined by eqs.

(2.14) and (2.13), respectively.

The Crout algorithm requires about n^/3 multiplications to accomplish this

decomposition so that for a 246 x 246 matrix about 5 million multiplications

are required. If an error message IER = 129 is printed by FLANK, this came

from LUDATF indicating that the matrix is "algorithmically singular” indicating

a small pivot element was encountered in the decomposition process. If the

message IER = 34 is printed, the computed solution is in error by more than can

be accounted for by the uncertainty of the data.

LUDATF, thus, triangularizes the matrix K into _L and Ut . The decomposition is

always possible provided the leading principal minor determinants are all

nonzero, i.e.

11
det kji + 0

IMl

1

k i2

l k21 k22

det K * 0

This triangular decomposition is unique if the diagonal elements of either L_ or

IJ are specified.

To obtain a first estimate of the temperature vector Tj , set

L U T i = B

and then

L X = B

X can be solved directly by forward substitution by noting

i-1

X± = (bi - Z Xii XO/Xu

J-l

The temperatures are then obtained by back substitution of the matrix equation

U T i = X

Once the decomposition has been accomplished, LUELMF calculates the temperatures,

so that a first estimate Tj is obtained as indicated. LUREFF is then called to

refine the solution. This routine basically works by calculating a residual rj

from

LI = b - K Ti (2.15)

The new solution is obtained (using the same decomposition from LUDATF) from

I XI = LI (2.16)

which is then added to Tj to obtain a better estimate T^, i.e.,

12 = !l + Z1 (2.17)

12
This computation requires about 2n^ multiplications. This procedure is continued

by now obtaining

£2 = b - K T2 , (2.18)

computing a y^?> obtaining a T/3 , etc.

Basically, the procedure continues until the accuracy specified by IDGT is

satisfied or an error is produced. Again IER = 129 means the matrix is too illconditioned to achieve the desired accuracy. For example, if IDGT is set to

three, an accuracy test is performed to be sure that when KT is calculated, the

first three digits of each element is the same as the first three digits of _b.

If IDGT = 0, this test is bypassed.

2.1.1 Heat Transfer Coefficients

The heat transfer coefficient must be read into the code. The CHB permit various

air velocities and temperature conditions to be imposed on the specimen. The

range of heat transfer coefficients is between 1.0 to 25 W/m3 »K. The heat transfer correlation used for laminar flow is [7]

where

and

h = 0.664 /Re Pr 1 / 3 (k/L), (2.19)

Re = the Reynolds number,

Pr = the Prandtl number,

k = the thermal conductivity of air (W/m» K)

L = the length of the wall in the direction that

the air is moving,

h is the average convective heat transfer coefficient (W/m3, K)

.

This expression is used if Re < 5 x 10 3 .

If the Reynolds number is greater than 5 x 10 3

, the flow becomes turbulent and

there are mixed boundary layer conditions in that for some of the wall specimen

the flow is laminar and turbulent flow dominates for the remainder of the wall.

In this case

h

L

“ i [ /

C

^laminar dx + / turbulent dx) < 2 - 20 > O Ac

and it is necessary to calculate the transition distance Xc

. The result of the

integration is

h

L

= (£) [0.037 Re4/5 - SSbjPr 1 / 3 (2.21)

Figure 2.4 is a plot of the forced convection heat transfer coefficients. These

will apply to both the MC and EC.

13
If the air velocity is very low, then natural convection becomes important and even though there is slight forced air motion, both types of convection must be

considered

Air infiltration is not taken into account in the heat diffusion equation. The

FRP is essentially impermeable to air and moisture but the wall specimen may

not be. FLANK would require extensive modifications for the specimen region in

order to take into account infiltration since convection may be an important

aspect of the heat transport mechanisms.

2.2 STEADY-STATE ANALYSIS

FLANK can be used either for steady-state calculations or dynamic calculations.

From equation (2.5) it is seen that if a is set to 1 and At is very large

(1C)30), then the static equations for heat transfer result. In FLANK, if

ITRAN = 0, these steps are automatically performed.

The boundary conditions are of the convective type or else the temperature may

be specified at various boundary points.

The equation that is solved in the steady-state case is

where JB arises because of the boundary conditions imposed.

The output from the static calculations consists of:

a) echo of the input data;

b) the temperature at each of the nodes;

c) the heat transfer across various surfaces; and

d) the frame coefficients.

2.3 TRANSIENT ANALYSIS

If the control parameter ITRAN is set to unity, a dynamic calculation will be

made. The initial temperatures and the boundary conditions imposed must be

part of the input file. In this case equation (2.6) is rewritten as

K T = B ( 2 . 22 )

[ L C + ctK](DT) p = aB( ?+1 > + (1 - a)B< p ) - K Tp

At = _ _ — — ~ — (2.23)

or

A(DT) P = _CP (2.24)

The boundary conditions may be changing in. time so that JB may be a time-dependent quantity. Equation (2.24) is solved for (DT) P and then the

actual temperatures are obtained at time (P + 1) by using the relation

T(P+1) = T( p ) + (DT) P (2.25)

14
I

I

The output of FLANK is as follows: a) echo of input data; b) the temperature

at each of the nodes at each time step; c) the heat transfer across various

surfaces at each time step; and d) the frame coefficients at each time step.

15
16

Figure

2.1

Geometry

(Region

A)

of

the

CHB

frame,

wall

specimen,

gaskets,

and

portions

of

the

MC

and

EC

walls
NZWALL

REGION 3 co

•Hbo

<u

t-. JZ

CJ

CO0)

03

<U 73Oc

4-1ou

a> JZe3

17

Figure

2.2

The

various

regions

and

them

for

the

CHB

as

coded

in

FLANK
(NTYWALL,

NTZWAL)

18

Figure

2.3

Nodal

structure

of

the

frame,

specimen,

gaskets

and

portions

of

the

MC

and

EC

walls
25.0

w •

O...O o

o «n o CN r-t —4

(X -

z m/n) ^uayoTjjaoo aajsueix

t Q

Figure

2.4

Heat

transfer

coefficient

due

to

forced

convection

that

applies

to

the

metering

and

environmental

chamber

specimen

and

frame

surfaces
3. RESULTS

3.1 CODE VERIFICATION AND DISCUSSION

FLANK was run for various situations to obtain an indication that it is working

properly. For example, if all the temperatures imposed on all boundaries are

equal, then the frame and specimen must go from the initial temperature to the

temperature specified by the boundary conditions. This test was done and the

program provided the correct answers.

Various other tests were done such as varying the node structure to determine

the "minimum acceptable" number, looking at the time steps and the accuracy of

the calculation, etc. Additionally, heat balances for given regions are performed to determine whether the heat entering a region is equal to that leaving

in a steady-state calculation.

The node structure used in calculating all the results in this report is shown

in figure 3.1. The accuracy of the calculation as a function of the node size

is shown in figure 3.2. Figure 3.3 gives the surfaces for which the heat fluxes

are calculated. It is found that 240 nodes or greater (if they are chosen

reasonably) will provide good results. Of course, the more nodes chosen the

greater the memory requirements. The time step used provides different accuracy

for different nodes. Generally, if the temperature of a node does not change

rapidly (as in the center of the polyurethane frame), longer time steps can be

used. Figure 3.4 shows the effect of the time step for node 121 which is one

of the most sensitive nodes since the FRP changes temperature rapidly. These

results are essentially the same for a = 1/2 and a = 1.0. This means that for

highly accurate results in the FRP, time steps of the order of 5 minutes or so are required. However, the accuracy in the polyurethane region is good even

for time steps of the order of an hour or so.

3.2 STEADY-STATE RESULTS AND ANALYSIS

The frame coefficient correction to the total energy balance of the CHB will

generally be less than about ten percent of the total heat entering or leaving

the MC side of the CHB. However, this is not necessarily the case since as Tj^c

approaches Tg^ with T^mb different, then essentially no heat flows through the

wall specimen and the frame coefficient can approach unity as shown in

figure 3.5.

An overall heat balance on the MC is given by the relation (assuming T^c > Tgc > tAMB

where

Qin Qout Qspecimen + Qframe (3.1)

Qin = The watts supplied by the heater and the blower minus the

watts taken by the chiller,

Qout = The heat that leaves the metering chamber.

20
It is assumed that the heat lost through the rest of the metering chamber walls

is negligible since they are maintained at a constant temperature by coils on

the outside. If this is not the case, another term QMCwalls must be added to

equation (3.1). The heat loss through the MC walls is given by the relation

QMCwalls = ^McA-eff -4= (3.2)

lMC

where k^c = the effective thermal conductivity of the MC walls (W/m-K)

Lmc = the thickness of the MC walls, m

Tamb = the temperature on the outside of the walls

A

e ££

= the effective area normal to the heat flow (m^).

It is assumed that the surface resistances are negligible compared to the

resistance of the walls in equation (3.2). Using conduction shape factors to

determine Agff [8] , the shape factors for the edge of two adjoining walls is

0.54 x (length of the edge) and for a corner, the shape factor is 0.15 x (the

thickness of the wall). Therefore,

5 8

Aeff = 2 Ain,i + 0-54 lMC * 2 e i + 4 * (0.15 )Lmc (3.3)

i=l i=l

where = length of edge i (m) and the 4 in equation (3.3) is because there

are four corners in the MC. The inside dimensions of the NBS-CHB MC are 1.111 m x 4.572 ra. Therefore,

Ae f f = 13.94 + 2(5.079) + 2(3.386) + 0.54(0.4572) • [2 x 4.572

+ 2 x 3.048 + 4 x 1.11] + 0.60(0.4572) = 36.00 m2

Assuming that there is a 1 °C difference between the inside and the outside of

the MC box,

QMCwalls

0.022 W (36 m2

) x 1 ( C)

m* °C 0.457 (m)

= 1.73 watts

Going back to equation (3.1), assuming temperatures Tjici and Tg^l on the specimen

surfaces (for a uniform temperature),

Qframe = Qin “ Qspecimen

=

Qin [l_ ^specimen

j (3.4)

Qin

= Q±n • FCHXS

,

21
if

Qspeciraen ~ (^ECl - T{^cj,)/AYg (3.5)

or

Qfratne = Qin * FCMCA , (3.6)

if

z wall

Qspeciraen = ^*572 / [T(z) - T^c] h(z) dz , (3.7)

z FRP

where the 4.572 is the width of the wall specimen in meters. T(z) is the surface

temperature of the wall as a function of height on the MC side of the specimen.

If QMCwalls I s not much smaller than Qframe* t^en eqs. (3.4) and (3.6) must be

modified to take that into account, i.e.,

Qframe Qin “ Qspecimen QMCwalls

- Qln (1 - Vclrnen

, . QMCwall3 . (3.8)

Qin

The program FLANK computes the frame coefficients FCHXS and FCMCA. Various

steady-state runs were made to obtain the temperature distribution, heat fluxes,

and frame coefficients. All runs were made using 246 nodes. Figure 3.6 illustrates the variation of the frame coefficients with the thermal conductivity of

the wall specimen. In this case, the specimen was placed 0.127 m (5 in) from

the MC side of the frame.

Figure 3.7 shows the influence of the position of the specimen on the wall

frame. The thermal conductivity of the wall specimen was taken to be 0.036

W/m ,0 K for this series of runs. The frame was 0.3048 m (12 in) wide. The

frame coefficients increase with position from the frame edge as expected. The

thickness of the specimen is 0.1016 m (4 in).

The thickness of the wall specimen also will influence the frame coefficient.

Figure 3.8 shows FCHXS and FCMCA both increase with wall thickness as expected.

The 0.3048 m wide (12 in) frame was used. It appears that FCMCA is starting to

bend toward a less than linear increase with thickness.

The frame coefficient as defined in this report is, strictly speaking, a steady-state concept even though it will be examined in time as well. If the

boundary conditions vary in time or if the system is not initially in equilibrium, the transient or a dynamic analysis must be performed.

3.3 TRANSIENT ANALYSIS AND RESULTS

A transient analysis differs from a dynamic analysis in that the latter has

driven or time-dependent boundary conditions. A transient analysis examines

the temperature response due to the fact that the system does not start out in

an equilibrium temperature distribution. The system approaches the steady-state

22
asymptotically so that "time constants" of the system can be obtained. A

series of runs with this objective were made.

Figure 3.1 shows the node numbering system used for all the transient

calculations shown. Figure 3.9 shows the temperature at nodes 68 and 226 as a function of time. The initial temperatures were all assumed to be 20°C and at

time zero, an MC temperature of 38° and a room temperature of 25°C was impressed on the system. The climatic chamber was taken to be 20°C also.

These boundary conditions were held constant in time and the equilibrium temperature approached as shown. Note that in the specimen (made of polystyrene)

the temperature reaches equilibrium in about two or three hours whereas deep

in the polyurethane, it takes about 50 hours to reach equilibrium.

If one examines the thermal dif fusivities , a, one can see the reason for these

times. The thermal diffusivity of the polyurethane is k/pc (see table 1.6)

which is (0.022 w/m» K)/(32 kg/m^) (921 J/kg» K) or 7.5 x 10”^ m^/s and

0. 036/(27. 64)(1214) or 1.08 x 10“6 m^/s is the thermal diffusivity for the

beaded polystyrene. The Fourier number is a dimensionless time and is

defined as

Fo = at/L

c

2 (3.9)

where Lq is the characteristic length of the system. The time dependence can roughly be related to e _S^°

where Bi is the Biot number (= h L^/k) . For the

frame Bi - 2.33 (0.04)/0.022 = 4.2 using the thermal conductivity of the polyurethane instead of the FRP. The time constant ("1/e" value) is thus

Fo • Bi = 1 (3.10)

or for Lc = 0.229 m, the time to reach about 60-70 percent of the steady-state

value is 5 hours. This is a very crude estimate but even without using FLANK,

it can be seen that time lags of the order of tens of hours will be required

for the polyurethane to reach equilibrium.

The approach of the frame coefficients to their equilibrium value is shown in

figure 3.10. In this case, the frame and specimen were initially at 20°C and

suddenly the metering chamber temperature was switched to 38°C and the ambient

(room) temperature was at 25°C. The time behavior of the frame coefficients

is shown to decrease to their equilibrium (steady-state) values in about 40

hours. Initially a large amount of heat is required to heat the frame material

itself and then, after steady-state temperatures are achieved, the frame

coefficient approaches the steady-state value.

23
235 I

L_J

J

240

o

O'

ft

VO VOft

O'oCM

oo m oo m cm

24

Figure

3.1

The

node

structure

used

in

the

calculations

presented
Ratio

of

Heat

Fluxes

to

the

344

Node

Case

0 100 200 300 400 500

Figure 3.2 The influence of node structure on the accuracy of heat

fluxes across various surfaces

25
HT

(6) (\

HT(5)

26

Figure

3.3

Heat

fluxes

numbering

system
ao

a>

4-1 3c

a

4-1 to

<u£

T*H

O O

-a- -j- ao

CO

vO <r <ni co c* co oCO

(0o

) (121 3Q0N) I

27

Figure

3.4

Sensitivity

of

the

temperature

in

the

FRP

to

the

time

step
o

o

o

oo o o o

vo <r oi o

(Z) ^uajOTJjaoD

28

Figure

3.5

The

influence

of

a

variation

of

the

climatic

chamber

temperatures

on

the

frame

coefficient
29

Figure

3.6

The

variation

of

the

frame

coefficients

with

the

thermal

conductivity

of

the

wall

specimen
(Z) auajojjjaoD aareii

30

Figure

3.7

The

frame

coefficient

as

a

function

of

position

from

the

metering

chamber

face

of

the

0.3048

m

wide

frame
TMC

»

38°C

(2) auaTOTjjaoD auiBJi

31

Wall

Thickness

(Meters)
o

oVO

o«n

o

oc*9

o

CN

O

09U

aoX

e

•H

O X n cs

'O «» tN © <S CM Cst

Oo

) aan^Bjadmax

32

Figure

3.9

The

transient

response

of

the

polyurethane

frame

and

the

specimen

to

abrupt

changes

in

boundary

conditions
(X) 3U3J3TJJ3O0 aarea.1

33

Time

(Hours)

Figure

3.10

The

transient

response

of

the

frame

coefficients
4 . SUMMARY

An initial analysis of the NBS-CHB has been done using the code FLANK (see

Appendix). The temperature distribution, heat fluxes across various boundaries,

and frame coefficients have been calculated for steady-state and transient

situations

.

The time constants for the frame coefficients are of the order of 40 hours.

The heat transfer coefficients are held constant during all transient analysis

but it is believed that this caused little error in the results. The frame

coefficient will "normally" be less than 20 percent unless the temperatures on both sides of the specimen are within about 5°C of each other.

The code can be used to calculate the coefficient for any specimen but the more deviation from uniformity the specimen the greater the error will be. Currently,

it is necessary to assume or calculate an effective thermal conductivity for

the specimen. One area of improvement for the calculations is to incorporate

into FLANK a more general model for the specimen.

The code is currently general enough to calculate the various thermal quantities

of interest for any homogeneous specimen subjected to constant or time varying

temperature boundary conditions. The specimen can be located anywhere on the

frame and be of any thickness less than the size of the frame itself.

34
5 REFERENCES

1. Achenbach, P. R., "Design of a Calibrated Hot-Box for Measuring the Heat,

Air, and Moisture Transfer of Composite Building Walls,” in Proceedings

of the ASHRAE/DOE - OENL Conference on Thermal Performance of the Exterior

Envelope of Buildings , December 3-5, 1979, Kissimmee, Florida, p. 308.

2. American Society for Testing and Materials, 1984 Annual Book of ASTM Standards , C976, (American Society for Testing and Materials 1916 Race Street, Philadelphia,

PA, 19103, 1984), Part 18.

3. Incropera, F. P. and DeWitt, D. P., "Fundamentals of Heat Transfer," John

Wiley and Sons, 1981, pp. 327-329.

4. Lavine, A. G., Rucker, J. L., and Wilkes, K. E., "Flanking Loss Calibration

for a Calibrated Hot Box,” Research and Development Division, Owens-Corning

Fiberglas, Technical Center, Granville, Ohio, 1982, Private Communications.

5. IMSL Subroutine Set, "International Mathematical and Statistical Library."

June, 1982, a proprietary package from IMSL, Inc., Houston, Texas.

6. Stewart, G. W., "Introduction to Matrix Computations," Academic Press,

New York, 1973, pp 135, 206.

7. Holman, J. P., "Heat Transfer," McGraw-Hill Book Company, 4th Edition,

1976 p. 171.

8. Incropera, F. P. and DeWitt, D. P., "Fundamentals of Heat Transfer,"

published by John Wiley and Sons, 1981, pp. 140-145.

35

APPENDIX

USER’S MANUAL FOR COMPUTER PROGRAM - "FLANK"

P. J. Burns

Department of Mechanical Engineering

Colorado State University

Fort Collins, Colorado 80523

R. J. Onega

Department of Mechanical Engineering

Virginia Polytechnic Institute and

State University

Blacksburg, Virginia 24061

A-i
Preface

This appendix presents a computer program FLANK which was written in ANSI Standard

FORTRAN to calculate the temperature distribution and various heat fluxes

through the specimen support frame for the National Bureau of Standards Calibrated

Hot Box (CHB). The code is a 2-D transient, finite-difference formulation of

the heat conduction equation. The mesh spacing is determined by the user. Either steady-state or dynamic calculations are possible by changing an input

parameter into the code. Constant temperature or heat-flux boundary conditions

can be imposed on various boundaries of the CHB frame and/or specimen. The

transient calculations can be done purely implicitly, explicitly or anywhere

between these two methods. Running the code in an implicit mode will eliminate

any instability problems.

A-ii
TABLE OF CONTENTS

Page

1. DESCRIPTION OF THE PROGRAM FLANK A-l

1.1 INTRODUCTION A-l

1.2 SYSTEMS OF UNITS A-l

1.3 COMPUTER USAGE A-2

1.4 FILES A-2

2. THEORETICAL FORMULATION A-3

2.1 THE HEAT DIFFUSION EQUATION A-3

2.2 TEMPERATURE BOUNDARY CONDITIONS A-4

2.3 CONVECTION BOUNDARY CONDITIONS A-4

2.4 TRANSIENT SOLUTION SCHEME A-5

2.5 STEADY-STATE SOLUTION SCHEME A-6

2.6 ACCURACY A-6

2.7 STABILITY A-7

2.8 NUMERICAL LIMITATIONS A-7

2.9 ORGANIZATION OF THE PROGRAM A-8

2.10 RESTART CAPABILITIES A-8

2.11 DATA CHECK ONLY A-9

2.12 COORDINATE SYSTEM AND GEOMETRY A-9

2.13 NODAL DATA A-9

2.14 MATERIAL REGIONS A-9

2.15 HEAT TRANSFER COEFFICIENT REGIONS A-13

2.16 HEAT FLOW LINES A-13

2.17 ENERGY BALANCE REGIONS A-13

2.18 MATERIAL SUBREGIONS A-13

3. INPUT DATA STRUCTURE A-16

3.1 PROBLEM HEADING A-16

3.2 CONTROL DATA CARD A-16

3.3 DEBUG INFORMATION A- 18

3.4 HEAT BALANCE CHECKS A-18

3.5 HORIZONTAL NODE NUMBER SPECIFICATION A-19

3.6 VERTICAL COORDINATE SPECIFICATION A-19

3.7 HORIZONTAL NODE NUMBER SPECIFICATION A-20

3.8 VERTICAL COORDINATE SYSTEM A-20

3.9 HEAT TRANSFER COEFFICIENT INPUT A-20

3.10 MATERIAL INPUT A-21

3.11 TEMPERATURES BOUNDARY CONDITION INPUT A-21

3.12 TIMING INPUT A-22

3.13 INITIAL TEMPERATURE INPUT A-22

3.14 BASE PLOT DATA A-23

3.15 CONTOUR PLOT DATA A-24

3.16 TIME HISTORY PLOT INPUT DATA A-25

3.17 TRANSIENT DATA INPUT A-26

A-iii
TABLE OF CONTENTS (Continued)

Page

4. RUNNING THE PROGRAM A-27

4.1 INTERACTIVE EXECUTION A-27

4.2 ERROR TERMINATION A-27

4.3 USEFUL UTILITY ROUTINES A-28

4.3.1 The Rename Routine A-28

4.3.2 The Delete Routine A-28

4.3.3 The Print Routine A-29

4.3.4 The Compile Routine F0RTR7 C A-29

4.3.5 The Compile/Library Load Routine FORT7 CR A-29

4.3.6 The Linker LINK33 . CSS A-29

4.3.7 The DELELIB Utility Routine A-30

4.4 CRT MESSAGES A-30

4.4.1 Input Data Incorrect A-30

4.4.2 Underflow Errors A-31

4.4.3 Overflow Error A-31

4.4.4 Program Operating Messages A-31

4.4.5 Time Step Indication A-32

5. EXAMPLE PROBLEMS A-33

5.1 A STATIC PROBLEM A-33

5.2 A DYNAMIC PROBLEM A-33

6. REFERENCES A-34

7. OODE LISTING A-35

A-iv
LIST OF FIGURES

Page

Figure 2.1 Coordinate system showing the input dimensions required

for FLANK and the various material regions A-10

Figure 2.2 The node structure for the NBS calibrated hot box A-ll

Figure 2.3 The various node designations for each of the geometrical

regions of the wall specimen and frame A-12

Figure 2.4 Heat transfer coefficients and associated temperature

designations A-14

i Figure 2.5 The heat flow lines and heat balance region designations ... A-15

A-v

1. DESCRIPTION OF THE PROGRAM FLANK

1.1 INTRODUCTION

The program FLANK is a FORTRAN computer program devloped at the National Bureau

of Standards (NBS) in Gaithersburg, Maryland, during the summer of 1982. The

program was written to calculate temperature and heat transfer quantities associated with the flanking loss for the Calibrated Hot Box (CHB) at the NBS.

Researchers at the Owens-Corning Fiberglas Corporation [1] coined the terra "flanking loss" for the heat that short circuits around the specimen in the

CHB, therefore the name of the program became FLANK.

The program is written to calculate the 2-dimensional conductive heat transfer

in the CHB geometry. Thus, the geometry is tied to that of the CHB, but the

node structure is general in each of the rectangular elements comprising the

CHB geometry. All of the material properties are taken as constant, i.e.,

independent of time or temperature. Boundary conditions include specification

of temperature, or specification of heat transfer coefficients and environmental

(ambient) temperatures. The boundary conditions may be constant or functions

of time, but not functions of temperature. Either transient or steady-state

problems may be simulated.

The program has the capability to calculate a variety of specific heat transfer

quantities of interest to the experimentalist such as the heat transferred

through the specimen, frame coefficient*, etc. Note that the above assumptions

render the problem "linear", although if the heat transfer coefficients are functions of time, it is necessary to reformulate the £ and K matrices** (to be

defined later) at each time step.

This manual is organized to present the theoretical formulation of the problem,

to detail the input quantities necessary to run the program and to present example problems. Chapter 1 is an overview of the code. Chapter 2 presents the

theoretical formulation of the problem. Chapter 3 details the specific input

necessary to run the problem, chapter 4 presents some interactive commands and,

finally, chapter 5 presents two example problems.

1.2 SYSTEMS OF UNITS

Where units are specified herein, the SI conventions are used. This is

merely for purposes of illustration as any consistent set of units may be used.

Also, the user is asked to note that the calculations are done in dimensional

form in the program (i.e., no nondimensionalization of the geometry equations

is done).

* The Frame Coefficient is defined as the fraction of the metering chamber

(MC) heat flow that does not go through the specimen.

** Nomenclature: Two bares under a denote that it is a 2-dimensional matrix;

one bar denotes a 1-dimensional vector and the remaining quantities are sealars

A-l
1.3 COMPUTER USAGE

The program was written and tested on the Perkin-Elmer 3230 computer of the

Center for Applied Mathematics (CAM) at the National Bureau of Standards. The

various subroutines were written in ANSI Standard FORTRAN, and thus, should be

transportable to other machines. For implementation on another machine, the

main program, and the various FORTRAN subroutines must be compiled and loaded

together with the IMSL Library (3 IMSL linear equation solving routines are used in FLANK as well as an IMSL plot routine) to generate an executable task

controller. Alternatively, the task code already existing on the 3230 may be

executed interactively by typing "FLANK" on the user's CRT provided that the

card image file "FLANK. INP" (the input file) resides on disk. More details

concerning these points are presented below.

1.4 FILES

The following files are necessary to run the computer code on the Perkin-Elmer

3230: (1) FLANK. CSS - the command file, (2) FLANK. TSK - the executable task

code, (3) FLANK. INP - the input file, and (4) FTEMP.INP - the restart tempratures*

.

Upon execution, the following files are generated: (1) FLANK. OUT - the output

file, (2) IMSLS .OUT - the IMSL error message file, and (3) FTEMP.OUT - the

ending temperatures (see section 3.16). The logical unit association is as follows

:

File Logical Unit No.

IMSLS. OUT 2

FLANK. INP 5

FLANK. OUT 6

FTEMP.INP 1 FTEMP.OUT 3

* The file FTEMP.INP is only necessary to run the code if the restart option is

in effect (see section 3.16).

A2
2. THEORETICAL FORMULATION

2.1 THE HEAT DIFFUSION EQUATION

The diffusion of heat in the CHB is completely specified as follows:

3^T -3T

pcrft- - k +J.

2

I]

3z 2

- k II h(T-T ) on S,

3n = 00 1

and

T = T(t) on S 2

subject to T = T0 (y,z) at t = 0

( 2 . 1 )

( 2 . 2 )

(2.3)

(2.4)

where y and z denote coordinate directions and t denotes time. Equation (2.1)

is the heat balance equation (T is the temperature at position (y, z) at time t)

for an isotropic material with constant properties. Equations (2.2) and (2.3) are the convective and temperature specification boundary conditions on surfaces Si and S 2 , respectively. The entire boundary is composed of the lines

Si and S 2

» Equations (2.1) to (2.3) are the equations that are solved numerically

to yield temperatures at the nodal points versus time, viz:

T=T(y,z,t) (2.5)

The equations (2.1) to (2.3) are assembled element by element into a finite

difference format and are solved either implicitly (0 < a <, 1) or explicitly if

a = 0. In matrix notation, the capacitance matrix £ and the conductance matrix

K may be formulated as follows:

£

= PiCpi Azi/At = Cn (2.6)

where the ii denotes that £ is a diagonal matrix, and no sum is intended on the

i's, and, denoting either Ay or Az by AXi, the conductance matrix is

, AX, AX, =K = k,

J

{ — i+_ 1}=K Ax

j

Axi

ij (2.7)

The finite difference (both in space and time) form of equation (2.1) may thus

be written for all the nodes as:

£(Tt +At " It) "I 1

= B. ( 2 . 8 )

The right-hand side is the thermal load vector LB. No mention is made here of

the time at which the T for the conductive term is to be taken, as this issue

AwilL be addressed in section 2.3. The boundary conditions roust also be

implemented in the solution scheme. These points will be addressed in the

ensuing sections of this manual.

2.2 TEMPERATURE BOUNDARY CONDITIONS

Temperature boundary conditions are handled via a penalty function formulation

for generality and ease of implementation. This approach also saves the expense

(in coding time) of reformulating the equations by taking the known temperatures

to the right-hand side (to the load vector). Thus, if N total node points are specified and M temperature boundary conditions are specified, this approach

still requires that N equations always be solved and not N-M. Thus, an increase

in formulation and solution time is accepted for ease of implementation.

A brief explanation of the penalty function follows. Suppose that node L is to

have its temperature set to T®. In the present scheme, the capacitance matrix

has the following term added to the diagonal of row L:

C(L,L) = C(L,L) + BIG (2.9)

where BIG is the penalty function (for the present program BIG = 10^0) . Since

the precision of the present machine is about 10^, BIG completely overwhelms

the original values on that diagonal such that in memory,

C(L,L) = BIG (2.10)

The same strategy is adopted for the load vector _B, i.e.:

B(L) = B(L) + BIG T^ (2.11)

B then is overwhelmed also so that the equation for node L now appears in memory

as follows:

BIG T(L) = BIG T® (2.12)

which has the solution that T(L) = T® as desired.

The user is simply required to input (1) the vector of node numbers that have

temperature boundary conditions and, (2) the vector of the values of the

temperatures at the boundary condition nodes. This specific input is described

in section 3.10.

2.3 CONVECTION BOUNDARY CONDITIONS

The convection boundary conditions contribute a terra to the diagonal of the

conductivity matrix, and a term to the load vector. The contribution to the

conductivity matrix is:

K(I,I) = K(I,I) + hx AXX , (2.13)

A-4
and the contribution to the load vector _B is

B(I) = B(I) + h x T*AXi . (2.14)

where T* is the ambient temperature and AXj is the length of the surface element

exposed to the heat transfer coefficient, hj. The user is required to input

the six values of the heat transfer coefficients and ambient temperatures as defined in section 2.15. The specific input for heat transfer coefficients is

described in section 3.9.

2.4 TRANSIENT SOLUTION SCHEME

Letting the solution matrix, C* , be defined as

C* = C + K + BIGii + (hAXi)!! (2.15)

where the (*) indicates a formuation as appropriate for the solution phase, and

equations (2.1) to (2.3) are to be resolved consistent with equations (2.6) to

(2.13). This is accomplished as follows.

Define the temperature multiplying the conductivity matrix at the a-point in

time, viz:

K T = K { a Tt+At + (1-a) Tt } (2.16)

Defining

AT = Tt+At " It (2.17)

The conduction equation becomes

CuT+aKAT=-KTt +B, (2.18)

were the terra on the right-hand side is known. Incorporating the solution

matrix which includes the convective and the temperature boundary conditions

leads to the expression

C* A T = - K T

t + £ hj AXj T? + Z BIG (ATB

) ^ (2.19)

where the summations indicate assembly into the proper row, j , of the load

vector. Note that the penalty function now is multiplied by the incremental

temeprature, AT?. It is convenient to define the load vector which includes

the convection and temperature boundary conditions, and the previous time step

conduction contribution, viz:

B* = - + Zjhj AXj Tj + E BIG(AT?) (2.20)

Thus Eq (2.19) can be written in the shorthand notation as

C* AT = B* (2.21)

A5
At the end of each time step, the temperatures are updated via equation (2.17).

There are two solution modes that the transient algorithm can adopt. The first

mode to be discussed is the one where the ambient temperatures and temperature

boundary condition temperatures are constants, i.e., they are independent of

time. In this case, JB* is constant with time except for the K terra^ and

this constant part is formulated (assembled) at the beginning of the solution

phase of the program and never reformulated. This assembly is done only once

but K_Tt is subtracted at every time step. This solution mode is invoked by

settling ITRAN = 1 on control card 1.

The second mode to be discussed is the one wherein the ambient temperatures and

temperature specifications change with time. In this instance, the solution

matrix ji* is completely reformulated every time step. If this solution mode

is desired, it can be invoked by setting ITRAN = -1 on control card 1. During

this mode of operation, new values for heat transfer coefficient boundary

temperatures and surface temperatures must be provided in the input file for

each time step. The input necessary to run a transient problem is presented

in section 3.12. No provision has been allowed for the heat transfer coefficients to change with time as this requires reformulation of the C* matrix at

each time step. This can consume considerable computer time so this solution

technique becomes prohibitively expensive.

2.5 STEADY-STATE SOLUTION SCHEME

Setting ITRAN = 0 on control card 1 will yield a steady-state solution. During

this mode of operation, the transient solution scheme is applied for one time

step, the previous time step nodal temperatures are all set equal to zero (0),

and the integration parameter, a, is set to one (1) so that the incremental

temperatures become the steady-state temperatures. The capacitance matrix _C is effectively removed from the calculation by dividing through by 10^^, and

the time step is set to one (1) so that the heat transfer output becomes a rate instead of a difference.

2.6 ACCURACY

Accuracy is a function of degree of the time discretization and the spatial

discretization of the problem. For the fixed geometry, the spatial discretization is determined solely by the number of node points in the various regions.

The user is encouraged to successively increase the number of node points

until the quantities in which he is most interested (e.g., heat fluxes or temperatures) change between successive mesh refinements by much less than what

the user considers an acceptable amount. This always involves playing off

computer time and memory versus accuracy. The user should be aware that smaller

mesh spaces result in a quadratic decrease in the stable time step if the time

integrated parameter, a < 1. In general, the increase in accuracy in a region

will vary between linearly and quadratically with the number of mesh points in

that region.

Time integration errors are dependent upon both the time step and the time

integration parameter. The time step should be reduced just as the mesh must

A-6
be refined, and the comments in the previous paragraph apply to this point also.

The time integration parameter also influences the accuracy. At a = 0

and a =1, the time integration scheme is accurate to order At. The accuracy

gradually shifts from At to At^ as a is moved from either a = 0 or from a = 1 to a = 1/2. The scheme for a = 1/2 is generally referred to as the

Crank-Nicholson scheme [2], and is recommended provided the time step is

already small enough to be stable.

2.7 STABILITY

To avoid a geometric increase in the magnitude of the temperatures (and to

conform to the second law of thermodynamics), the time step is limited to

At < fin 1-^ (AX) 2

1 k

1 ( 2 . 22 )

where the "min" denotes the minimum value overall the grid points and thermal

conductivities

.

Actually, equation (2.22) only applies at the interior mesh points. At the

external mesh points, the convective coefficients exert a debilitating influence

on stability. However, the following blanket statement applies: there are no stability problems whatsoever at a = 1. So, if stability is critical, the user

is encouraged to set a = 1 and to accept the slightly lessened time integration

accuracy.

2.8 NUMERICAL LIMITATIONS

The calculations are being performed on a 32-bit machine in the single precision

mode. The structure of a floating point number limits the accuracy to about 7 decimal places and the limit of the magnitude is about 10^5. in some instances,

7 places of accuracy will be lost during the IMSL solution of the linear set of

equations even though full pivoting is allowed. If this occurs, the IMSL routine will return a nonzero value of the IER parameter, and the program will

automatically terminate after writing a brief explanation of the error condition

into the output file FLANK. OUT. A more complete description of the problem can usually be found in the file IMSLS.OUT wherein the IMSL routine writes its own internally defined comments.

If this situation occurs, the user has two options: (1) decrease the number of

node points and rerun the problem or (2) reload the program in double precision

mode

.

During the course of decomposing the solution matrix, the IMSL routine LUDATF

computes the determinant to test for singularity (when the determinant = 0).

It is very likely that the determinant will exceed 10^5 i n value and generate

an overflow of the exponent register (this does not affect the solution as this

value is never used). When this happens, an interactive message is sent to the user at his CRT. It is likely that this is the place where the overflow occurs,

but it might not be the only place. Consequently, the user is encouraged to

check his answers very carefully whenever this happens.

A7
It is also likely that an underflow will be generated during the decomposition

and/or refinement phase of the IMSL solution. A message to this effect is

displayed upon the user's CRT and may be ignored.

2.9 ORGANIZATION OF THE PROGRAM

The program is logically organized into 3 basic elements. These elements are capable of being overlayed, but were not as very little can be gained from

this. The 3 elements are as follows:

(1) Main routine

(2) Input routine

(3) Solution routine

The main routine sets the storage for the common arrays, vectors and constants,

and then controls the program flow by calling the input and solution routines

sequentially.

The input routines consist of a variety of subroutines to read and echo all of

the input parameters. The reader may consult the code listing for more specific

information.

The solution routines perform a variety of functions, including the following

basic ones:

(1) assembly of the solution matrix, C*

(2) assembly of the load vector, B*

(3) performance of the energy balances and tallies

(4) updating and printing the solution vector, JT(t+At)

and

(5) solving the linear algebraic equations including:

(a) decomposition of the solution matrix, C*

(b) back substitution to obtain I’Ct+At)

and

(c) refineraenet to obtain a more accurate T(t+At).

For the finer details, the user is referred to the code listing.

2.10 RESTART CAPABILITIES

Occasionally, the user may wish to run the program in the transient mode for a short period of time, stop, examine the answers and start the program from

where the last calculation stopped (restart). In order to facilitate this

process, FLANK has the option to write a file of nodal temperatures upon

termination (file name = FTEMP.OUT). This file may later be renamed and read

as the initial temperatures for the subsequent run (file name = FTEMP.INP).

This capability eliminates the laborious manual editing necessary to feed these

A8
temperatures into the input file so that the user need only change the starting

and ending times (and some control parameters) in the input file. The specific

input required to exercise this option is detailed in sections 3.1 and 3.12.

2.11 DATA CHECK ONLY

One of the common errors encountered when using the computer occurs when the

data are placed incorrectly in the input file. This is particularly critical

when the input data control the memory location of subsequent input datums

(e.g., supposing that a temperature boundary condition is to be supplied to

node 5 in an 15 format, suppose further that the input is not right-justified

and that 50000 is read in lieu of 5, the value of the boundary temperature read

subsequently will then be stored in the location of the start of its array +

50000 and this will certainly exceed the length of the array with the result

that something, either a storage value or even executable code, will get overwritten). To aid the user in debugging such errors, FLANK has the capability

to read in and echo the input data, and to then stop so that the user may print

out the output file and check it for such errors. This capability is invoked

by setting IDATA = 1 on control card 1 as is described in section 3.2.

2.12 COORDINATE SYSTEM AND GEOMETRY

The coordinate system for the problem is shown in figure 2.1 with X]_ = y and

X2 = z. The origin exists as the bottom left-hand corner of the base of the

CHB. The geometry is 2-dimensional with a horizontal y-axis and a vertical z- axis. All of the dimensions shown in the figure are adjustable and must be

input by the user as described in sections 3.6 and 3.8.

2.13 NODAL DATA

The associated nodes are shown in figure 2.2. Note that the total number of

nodes must be specified in the first horizontal and the first vertical sections (NZPOLY and NY5, respectively) so that nodes are specified at the starting

boundary line and the ending boundary line in these sections. In all other sections, one line of boundary nodes already exists, from the specification of the

previous section. Stated another way, in the first horizontal section and first

vertical section:

Entering N nodes = N-l grid spaces

In all other sections

Entering N nodes = N grid spaces.

Input of this type is detailed in sections 3.5 and 3.7.

2.14 MATERIAL REGIONS

The material regions are also shown in figure 2.1. Normally, the following

materials will be specified by number:

A-9
Figure 2.1 Coordinate system showing the input dimensions required for FLANK

and the various material regions

A-10
(rrm.irrzwAL)

(l.HTZFRP) (NTY1L.NTZFRP

(i.hzpoly;

( 1 . 1 ) (HY5.1) (NTYU.l) (HTYlL.l)

( HTYWALL .HTZUAL

)

( IfTYWAL .NTZFRP

)

(HCOL.RTZFRP)

(NCOL.HZPOLY)

(MTYWAL.l) (miR.l) (rna.i) (1C0L.1)

Figure 2.2 The node structure for the NBS calibrated hot box

A-l 1
Figure 2.3 The various node designations for each of the geometrical

regions of the wall specimen and frame

A-12
Materials No IXEe

12

3

4

5

6

Polyurethane

Fiberglass Reinforced Polyester (FRP)

Fixed Gasket (Neoprene)

Inflatable Gasket

Wall

Null

It is required to input values for the thermal conductivity, the density and

the specific heat. Only the first five (5) materials may be input; the sixth

is internally set. This is detailed in section 3.10.

2.15 HEAT TRANSFER COEFFICIENT REGIONS

The heat transfer coefficients are shown in figure 2.4 by region. Note that a counterclockwise sequencing has been adopted. Normally, the heat transfer coefficents for regions 2, 5, and 8 will be set to zero to represent an insulated

(zero heat flow) condition. The heat transfer coefficient input data are described in section 3.9.

2.16 HEAT FLOW LINES

There are 13 lines along which heat flows have been defined. These lines are shown in figure 2.5 where the heat flow has been defined as positive in the

positive direction of the coordinate axis (except for convective heat transfer).

The input to invoke the energy balance routines is described in section 3.2 and

3.4.

2.17 ENERGY BALANCE REGIONS

There are 5 regions wherein steady-state energy balances may be performed.

These are also shown in figure 2.5 as enumerated areas. The user is cautioned

that specifying a temperature boundary condition on or in a region adds or subtracts energy, and that this is not explicitly accounted for in the energy

balance routine resulting in a large net flow for that region (i.e., the energies

seem not to balance for this reason)

2. 18 MATERIAL SUBREGIONS

The subregions are the result of the grid application in that everywhere an area is bounded by coordinate lines, the grid inside the coordinate lines is

regular. These regions are smaller than the material regions detailed in section 2.14, and, hence, have been termed subregions. The B, U, L, R and C

specification has been adopted for bottom, upper, left, right, and center

wherever possible. These subregions will ordinarily be transparent to the

user, but have been included here for completeness.

A- 13
*

5

. T5

Figure 2.4 Heat transfer coefficients and associated temperature

designations

A-14
Figure 2.5 The heat flow lines and heat balance region designations

A- 15
3. INPUT DATA STRUCTURE

The following sections detail the specific Input required to run the program.

The Input Is organized Into the following sequential areas:

(1) Problem Heading - 1 Card

(2) Control Information - 1 Card

(3) Constant Data - {(12 + NUMTBC) cards, steady-state; (13 + NUMTBC)

cards + Initial temperature data (variable number of cards), transient

(4) Transient input required If the heat transfer coefficients and ambient

temperatures change with time -2 (number of time steps) cards.

The control information is so termed because these data control both the number

of data items read and the execution of the program. The input follows in the

sequential fashion that is necessary for Input into the program.

3.1 PROBLEM HEADING

Conditions - None

Number of Cards to be read - 1.

Cols . Format Item(s) Note(s)

1-64 16A4 Problem heading to appear 1 on each page

Note(s): 1. The user may wish to include the time and date in the heading for

each problem.

3.2 CONTROL DATA CARD

Conditions - None

Number of cards to be read - 1.

Cols

.

Format Item(s) Note(s)

1-5 15 NUMTBC - number of temperature

boundary conditions, NUMTBC > 0.

6-10 15 NUMIT - number of initial

temperatures

NUMIT < 0 - uniform initial

temperature * TEMP0

NUMIT > 2 - NUMNP initial

temperature read from file

FTEMP.INP

1,2

A- 16
Cols Format Note(s)

11-15

16-20

21-25

26-30

31-35

36-40

41-45

46-50

15

15

15

15

15

15

15

15

Item(s)

IDATA - data check only code

IDATA = 0 - normal execution

IDATA £ 0 - data check only

IPR - print interval, temperature

output is printed every

IPR time steps

IRFINE - refinement control 3

parameter

IRFINE = 0 - no refinement is done

IRFINE ^ 0 - refinement is done

IDGT - IMSL accuracy test parameter 4

IDGT = 0 - no accuracy test is done

IDGT > 0 - accuracy test is performed

IEB - heat flow parameter

IEB = 0 heat flows or energy balances

are not computed

IEB ^ 0 energy balances and heat flows

are computed

IQPR - energy tally print interval,

energy tallies are printed every IQPR

time steps.

ITRAN - transient solution code 5

IRAN = 1 - transient problem, constant

parameters

ITRAN = 0 steady state problem

ITRAN = - 1 - tans lent problem,

variable parameters

IPLT - plotting parameter

IPLT = 0 no plots done

IPLT = 1 plots done

Note(s): 1. If ITRAN = 0, no initial temperatures are read as the initial

temperatures are automatically set equal to 0.

2. If NUMIT _< 2, the file FTEMP.OUT which includes the nodal point

temperatures is written. The name of this file can be changed to

FTEMP.INP (provided any existing FTEMP.INP is first deleted) which

can be read on a subsequent problem if NUMIT >2.

3. During refinement, the output vector obtained from back substitution

is iterated upon with the original solution matrix (the decomposed

A-l 7
solution matrix Is used to obtain the output vector by back

substitution) to achieve greater accuracy. This endeavor Is

usually worthwhile as the cost Is not great.

4. Upon decomposition, a check is made to determine if IDGT’s of

accuracy were retained in the decomposed mode. If IDGT’s of accuracy were not retained, an error termination results. This is not

a crucial test, as refinement tends to improve the answers apart

from the decomposition mode. Therefore, IDGT - 0 is recommended.

5. For ITRAN - -1, new values are required for all temperatures at the

temperature boundary condition nodes and for all ambient temperatures

in the enclosure. The capability does not presently exist to have

variable heat transfer coefficients although this can be handled as

a restart option (it is rather awkward, though).

3.3 DEBUG INFORMATION

Conditions - None

Number of cards to be read - 1.

Cols . Format Item(s) Note(s)

1-50 1015 IDBUG(IO) - debug control array. If 1

IDBUG(I) is different from zero, the

following information is printed.

I 1 K and matrices, JS* vector

1=2- decomposed C* matrix

I = 3 - nodal association data

I = 4 - element assembly matrices

1 = 5- interior data written

I = 6 - C* matrix written from equation

solver subroutine

I = 7 - JJ* vector written from equation

solver subroutine

1 = 8 -

1 = 9-

I = 10 -

Note(s): 1. Normally, all values should be set to zero (0), as this information

was required during the debug phase which is presently complete.

3.4 HEAT BALANCE CHECKS

Conditions - IEB ^ 0

Number of cards to be read - 1.

A-18
Cols Format Item(s) Note(s)

1-25 515 ICHECK(5) - energy balance control

array. If ICHECK(I) ^ 0, the following steady-state energy balance checks

are performed and printed, (see figure

2.5

for these heat balance regions).

1=1- region 5-6-4U-4L-1L-2L

1=2- region IC

1=3- region 1R-2R-7L-7U-8-9

1=4- region 2 C 1=5- region 3

3.5

HORIZONTAL NODE NUMBER SPECIFICATION

Conditions - None

Number of cards to be read - 1.

Note(s): 1.

2 .

3.

Cols

.

Format Item(s) Note(s)

1-10 E10.0 NY5 1,2

11-20 E10.0 NY4 3

21-30 E10.0 NY1L 3

31-40 E10.0 NYWALL 3

41-50 E10.0 NY1R 3

51-60 E10.0 NY7 3

61-70 E10.0 NY8 3

See section 2.13 for definitions of the nodal specifications.

This number represents the (number of grid spaces - 1) in the

region.

This number represents the number of grid spaces in the region

3.6

HORIZONTAL COORDINATE SPECIFICATION

Conditions - None

Number of cards to be read - 1.

Cols

.

Format Item(s) Note(s)

1-10 E10.0 Y5(m) 1

11-20 E10.0 Y4(m)

21-30 E10.0 YlL(m)

31-40 E10.0 YWALL(m)

41-50 E10.0 YlR(m)

A-19
Cols Format Item(e) Note(s)

(cont.)

51-60 E10.0 Y7(m)

61-70 E10.0 Y8(m)

Note(s): 1. The reader is referred to section 2.12 for definitions of these

coordinates (see figure 2.1).

3.7

HORIZONTAL NODE NUMBER SPECIFICATION

Condition - None

Number of cards to be read - 1.

Cols . Format Item(s) Note(s)

I-10 E10.0

II-20 E10.0

21-30 ELO.O

NZPOLY 1,2

NZFRP 3

NZWALL 3

Note(s): 1. The reader is referred to section 2.12 for definitions of the

node numbers.

2. This number represents the (number of grid spaces -1) in the

region.

3. This number represents the number of grid spaces in the region.

3.8

VERTICAL COORDINATE SYSTEM

Conditions - None

Number of cards to be read - 1.

Cols. Format

1-10

11-20

21-30

E10.0

E10.0

E10.0

Item(s)

ZPOLY(m)

ZFRP(m)

ZWALL(m)

Note(s)

Note(s): 1. The reader is referred to section 2.12 for definitions of the

coordinates (see figure 2.1).

3.9

HEAT TRANSFER COEFFICIENT INPUT

Conditions - none

Number of cards to be read - 2.

A-20
Cols

.

Format Item(s) Note(s)

(First Card)

1-80 8E10.0 HTC(I), heat transfer coefficients,

input sequentially by regions,

I(W-m"2 *K_1

)

.

1

(Second Card)

1-80 8E10.0 TAMB(I), ambient temperatures, input

sequentially by region, I (K).

1

Note(s): 1. The reader is referred to section 2.14 for definitions of the

heat transfer coefficient regions (see figure 2.4).

3.10 MATERIAL INPUT

Conditions - none

Number of cards to be read - 3.

Cols

.

Format Item(s)

(First Card)

1-50 5E10.0 RKAY(I), thermal conductivity, input

by material region, I (W»m“^*K~^)

(Second Card)

1-50 5E10.0 RHO(I), densities, input by material

region, I (kg»m-

^)

Note(s)

1

1

(Third Card)

1-50 5E10.0 CP(I), specific heats, input by 1,2

material region, I ( J*kg--*- # K“^)

Note(s): 1. The reader is referred to section 2.14 for definitions of

material regions.

2. Note that the units must be input as Joules and not

kiloj oules

.

3.11 TEMPERATURES BOUNDARY CONDITION INPUT

Conditions - NUMTBC j* 0

Numbers of cards to be read - NUMTBC/5.

A-21
Cols

.

Format Item(s) Note(s)

1-75 5(15,

E10.0)

(NODEST(I) , TBC(I) , 1-1 , NUMTBC)

(the node number at which the

1,2,3

temperature is applied, and the

applied temperature (K))

Note(s): 1. Temperatures may be applied at any nodes whether they are on the

boundary or not.

2. Node numbers vary sequentially along rows (by column) starting

with 1 at the origin. The node numbers are incremented by one

in all cases. The sequence varies by row through the entire

base. The next sequence starts at the lower left-hand of the

wall specimen and is incremented by one along the rows (again by

column) through the wall.

3. An example of nodal numbering is given in chapter 4.

3.12 TIMING INPUT

Conditions - ITRAN # 0

Number of cards to be read - 1.

Cols. Format Item(s) Note(s)

1-10 E10.0 TIM0 - initial time(s)

11-20 E10.0 TIMF - final time(s)

21-30 E10.0 DELTA - time step(s) 1

31-40 E10.0 ALFA - time integration parameter,

0 < ALFA _< 1

2

Note(s): 1. The number of time steps to be taken will be c (TIMF - TIM0)/DELTA.

2. The reader is referred to sections 2.6 and 2.7 for recommended

values of ALFA.

3.13 INITIAL TEMPRATURE INPUT

Conditions - NUMIT ± 0

Number of cards to be read - 1 or NUMNP/5.

Cols . Format Item(s) Note(s)

(Case 1, NUMIT = -1)

1-10 E10.0 TEMP0 - uniform initial temperature 1,2

(K), (1 card only read)

A-22
Cols. Format Item( s) Note( s)

(cont.)

1-75

(Case 2, NUMIT = 1)

5(15, (NODE, T(NODE) , 1=1, NUMNP) , node

E10.0) number and initial temperature of

that node (K), (NUMNP/5 cards read)

2

Note(s)

:

1. Every node point temperature is initially set ot TEMP0.

2. For reading initial temperatures from the disk file FTEMP.INP,

NUMIT should be set to 2. For writing final temperatures into

the disk file FTEMP . OUT, NUMIT should be set to -2. This file

may be renamed to FTEMP.INP and subsequently read. (See

sections 2.10 and 3.2.)

3.14 BASE PLOT DATA

Conditions - IPLT = 1

Numbers of cards to be read - 2.

Cols. Format Item( s) Note( s)

I-5 15

6-10 15

II-15 15

16-20 15

21-25 15

Card 1 - Plot Control

NUMCON - number of contour plots to be 1

plotted 0 _< NUMCON <_ 5

NUMHIS - number of time history plots 1

to be plotted 0 <_ NUMHIS < 5

IOPT - plot column option 2 IOPT = 0 - narrow plots (80 columns)

IOPT = 1 - wide plots (126 columns)

NOUT - plot output file specifier 3

NOUT = 2 - plots written to file

IMSLS. OUT

NOUT = 5 - plots written to output

file FLANK. OUT

NSTEPC - time step interval for contour 4

plots (5 time curves are accumulated

before the plot is written, this is

5 * NSTEPC time steps)

A-23
Note(s): 1 A maximum of 5 curves is permitted

2.

The IOPT 0 plots are especially convenient for tracing on

8 1/2" x 11" paper.

3. Reading an NOUT different from 2 or 5 will result in an error termination.

4. The contour plots are written every 5 time stpes or at the end

of the problem, whichever comes first.

Cols

.

Format Item(s) Note(s)

Card 2 - Range Data

1-10 E10.0 TMIN - the minimum temperature

for the vertical plot axis

value 1,2

11-20 E10.0 TMAX - the maximum temperature value 2,3

for the vertical plot axis

Note(s): 1. All plotted values should be above TMIN.

2. It is a good idea to allow some leeway in picking the range for

the plots as it is convenient to have some space between the

plotted values and the axes.

3. All plotted values should be less than TMAX.

3.15 CONTOUR PLOT DATA

Conditions - NUMCON ^ 0

Number of cards to be read - 5 x (NUMCON)

.

Contour plots are "snapshots" of the temperature profile along specified node

points at an instant in time. Five curves of temperature versus node point

number are plotted at 5 specified times (unless the last time step is taken, in

which case the plots are done before termination), on any one curve. A maximum

of 5 plots are allowed.

The following block of input items is required for each plot to be drawn. The

blocks should be input sequentially, until NUMCON of them have been input.

A-24
Cols

.

Format Card No. Item(s) Note(s)

1-72 72A1 1 ITITC(I) - plot title

1-36 36A1 2 ITITXC(I) - x-axis title

1-36 36A1 3 ITITYC(I) - y-axis title

1-5 15 4 NC(I) - number of nodes

against which temperatures

are to be plotted, 2 < NC(I)

< 20

1

1-80 1615 5 NODESC(I), the node numbers of

the temperatures to be plotted

2,3

Note(s)

:

1. A maximum of 20 nodes is permitted.

2. Actually, 2 cards may be required for this input in which case

6 x NUMCON cards must be input.

3. The nodes may be input in any order, i.e., the nodes need not be

adjacent.

3.16 TIME HISTORY PLOT INPUT DATA

Conditions - NUMHIS i 0

Number of cards to be read - 5 x (NUMHIS)

Time history plots are plots of the node temperature versus time. The abcissa

is fixed as a linear variation of time from the initial time to the final time.

Up to 5 node numbers may be included on each plot. A maximum of 5 plots is

allowed. Unlike the contour plots, only 5 plots are written to the plot file

(the contour plots allow up to 5 plots to be written every 5 x (NSTEPC) time

steps) .

The following block of input items is required for each plot to be drawn. The

blocks should be input sequentially until NUMHIS of them have been input.

Cols

.

Format Card No. Item(s) Note(s)

1-72 72A1 1 ITITH(I) - plot title

1-36 36A1 2 ITITXH(I) - x-axis title

1-36 36A1 3 ITITYH(I) - y-axis title

1-5 15 4 MH(I) - number of node points

to be included on this plot

1 £ MH(I) £ 5

1

1-80 1615 5 NODESH(I) - the node numbers

of the temperature to be

2, 3, 4,

5

plotted

A-25
te(s): 1. A maximum of 5 nodes is allowed.

2. The nodes may be input in any order.

3. Only 1 plot for each block of nodes is allowed for the entire

time interval.

4. The number of stored values is internaly limited to 80, maximum.

If the total number of time steps divided by 80 is not an integer

number, less data are stored (e.g., for 100 time steps, 50 data

points are stored).

5. The starting time temperatures and the ending time temperatures

are automatically included in the plots.

3.17 TRANSIENT DATA INPUT

Conditions - I TRAN - 1

2 sets of data to be read for each time step.

Cols. Format

1-80 8E10.0

1-80 8E10.0

Item(s) Note(s)

First Data Set - NUMTBC/8 cards to be read

(TBC(I), 1*1, NUMTBC) - temperature 1,2

boundary condition temperatures, input

sequentially (K).

Second Data Set - 1 card

(TAMBP(I), 1=1, 8) - ambient temperatures 2,3

at the end of the time step (K).

Note(s): 1. The nodes at which the boundary temperatures were set (NODEST(I))

remain unchanged.

2. Complete sets of these data must be input at each time step,

regardless of whether or not the values change.

3. The heat transfer coefficients (HTC(I)) remain unchanged (at

their initially input values).

A-26
4. RUNNING THE PROGRAM

4.1 INTERACTIVE EXECUTION

A Command Substitution System (CSS) file, FLANK. CSS, is available to run the

program interactively on the Perkin-Elmer 3230. The file loads the executable

task code, FLANK. TSK, into memory, allocates disk space for output files,

assigns local unit numbers, and prints the timing log at the user’s terminal

screen at the end of execution. Interactive execution is accomplished by typing

FLANK

followed by a (carriage) return.

The following files must be available for the program to read from:

(1) FLANK. INP - the input file as described in chapter 3.

(2) FTEMP.INP - the file of initial temperatures in a 5E14.0 format (this

file may contain only 1 blank line if no reads are required, i.e.,

NUMIT + 2.

The following output files are generated:

(1) FLANK. OUT - the output file with the echo of the input data, the

temperature states, the energy balance information, and the debug

information.

(2) FTEMP.OUT - the output file of the final temperatures written in a 5E14.6 format.

(3) IMSLS.OUT - the output file containing any IMSL debug information.

Should difficulty be encountered in running the computer code, the user should

locate and inspect the files FLANK. CSS, FLANK. INP, FTEMP.INP, and FLANK. TSK.

Usually, the problem will be that one of these files is missing, or the data

are input incorrectly in FLANK. INP. If problems are still encountered, the user should try a data check only (see section 3.2).

4.2 ERROR TERMINATION

The program internally checks for a variety of input errors. If any input errors are detected that are of such a nature as to prohibit the execution of

the program, an error termination is invoked. Then, a descriptive message is

printed in the file FLANK. OUT, together with a line denoting:

" ERROR TERMINATION ”

A-27
One of the most common mistakes is easy to decipher. This mistake occurs when

the user has input more node points than the allowed storage locations. The

following message is then printed:

"NO. OF NODE POINTS = " N1

"EXCEEDS DIMED VALUE = " N2

"STOP"

where Nl is the number of node points that the user has specified, N2 is the

maximum number allowable, and DIMED is an abbreviation for "dimensioned".

Should this occur, the user has only 2 choices: (1) reduce Nl to be N2 , or

(2) edit the FORTRAN source code FLANK. FTN to increase N2 (NNODES in FLANK. FTN)

and reload the program using the LINK33.CSS routine on disk.

4.3 USEFUL UTILITY ROUTINES

Some useful system utility routines that exist on the NBS Perkin-Elmer 3230

computer will be described in this section. All lines must be input exactly

as shown, followed by a carriage return. All lines contain parameters that

will also be described. Optional parameters will be included in parenthesis.

The parameter PR: is the system definition for the line printer, and when

included, will automatically send the output to the line printer. The utility

routines are enumerated, and the purpose of each routine is given.

4.3.1 The Rename Routine

Command Line: REN fnl .exl ,fn2 .ex2

Parameters: fnl. exl - file name fnl with extension . exl

fn2.ex2 - file name fn2 with extension. ex2 (fn2.ex2 must not

exist on disk).

Purpose: REN renames a file. This may be useful when a file is to be printed

(e.g., FLANK. OUT) for when a file is in the print queue, it cannot be deleted or rewritten. Also, the file of output temperatures may be renamed to be the

input temperatures for a subsequent run.

Examples: REN FTEMP .OUT, FTEMP.INP

REN FLANK. OUT, FI. OUT

4.3.2 The Delete Routine

Command Line: DE fnl. exl

Parameters: fnl. exl - the file fnl. exl is deleted from disk

Examples: DE FTEMP.INP

A-28
4.3.3 The Print Routine

Command Line: PRINTfnl.exl

Usage: PRINTF1.0UT

4.3.4 The Compile Routine F0RT7C

Command Line: F0RT7C fnl, (PR:), (XREF) , 10, (lib.exl)

Parameters: fnl - the file of FORTRAN source code fnl.FTN (note the extension

is not given on the command line, and it must be FTN) is compiled,

file fnl. OBJ (the object code) is generated and written to disk.

PR: - the compiler listing is sent to the printer

XREF - the cross reference map is sent to the printer

lib.exl - the library file lib with extension, ,exl=.LIB, is loaded

also. (For compiling FLANK. FTN, this space should be left blank.)

Examples: F0RT7C FLANK, , , 10

F0RT7C FLANK, , XREF, 10,

4.3.5 The Compile/Library Load Routine F0RT7CR

Command Line: F0RT7CR fnl, (PR:), 10

Parameters: fnl - the FORTRAN source code file fnl.FTN is complieted

PR: - the compiler listing is sent to the printer.

Purpose: The object code (compiler output)

fnl. OBJ is generated from the compiler. If this routine exists in

the library FLANK. LIB, it is deleted. The routine fnl. OBJ is then

written into the object code library FLANK. LIB. This object code

library is loaded together with FLANK. OBJ, the IMSLS library, and

the system FORTRAN library by the utility routine LINK33.CSS.

Example: F0RT7CR EQS0L2, , 10

F0RT7CR EQS0L2, PR: ,10

4.3.6 The Linker LINK33.CSS

Command Line: LINK33 FLANK, (PR:), FLANK. LIB

A-29
Parameters: FLANK - the main routine object code file FLANK. OBJ is loaded

PR: - a cross reference map of the load is sent to the printer

FLANK. LIB - the object code files of all of the subroutines are loaded from the object code library FLANK. LIB.

Purpose: This routine loads all of the object code, together with the IMSLS

library and the system FORTRAN library to generate the executable task code

FLANK. CSS.

Usage: LINK33 FLANK, , FLANK. LIB

LINK33 FLANK, PR:, FLANK. LIB

4.3.7 The DELELIB Utility Routine

Command Line: DELELIB f nl , fn2.LIB

Parameters: fnl - the object code file fnl.OBJ is deleted from the object

library fn2.LIB

fn2 - the object code library fn2.LIB

Usage: DELELIB FLANK, FLANK. LIB - this removes the routine FLANK. OBJ from

FLANK. LIB. Note that this file should not be in the library as it gets loaded

together with the library at link time via LINK33.

4.4 CRT MESSAGES

When executing the program, various error messages and other messages will

appear on the CRT terminal. If an error does occur, the user must then cancel

the job using the command

CA

The following examples are listed, together with an explanation.

4.4.1 Input Data Incorrect

If the input data are not in the correct format or if data are missing, the

following sequence is likely to be seen on the screen.

* FLANK - program initiation step

- IN S0LVE2 - AT BEGINNING

- IN S0LVE2 - CALLING ASEMB2

- IN S0LVE2 - ASEMB2 DONE, CALLING FORMB

- NON-EXISTENT SEGMENT ERROR (PST) AT B4C80

- MEMORY FAULT ADDRESS = FFFD30

TASK PAUSED

*

A-30
This error message resulted because the number of temperature boundary

conditions NUMTBC was 3 but only two temperatures were actually supplied to

the input file. The only corrective action necessary is to edit FLANK. INP and

change NUMTBC to 2 or else supply a third node specification and temperature

for a boundary point.

4.4.2 Underflow Errors

Often with specified temperature boundary conditions, the message

- FLOATING POINT - UNDERFLOW ERROR AT D70D2

- NEXT INSTRUCTION AT D70D8

will occur. The D70D2 are the hexadecimal memory locations for the particular

subroutine in which the underflow occurred. These errors can be ignored on the

P-E since penalty functions are being used and division by a very large number

occurs. The PE simply sets the result to zero.

4.4.3 Overflow Error

The determinant of the £ or K matrix will often be greater than 10^6 if a large

number of nodes is used. The typical value of a diagonal element is about 2 so that if there are 246 nodes, 2^46 = 10^4. This determinant value is computed

but never used so the error may be ignored. The CRT screen will display the

following message

ERR 517 (0C 94A6):

RXXR: EXPONENTIAL OVERFLOW

IN SOLVE 2 - DONE DECOMPOSITION, BEGIN TIME LOOP

This message can be ignored.

4.4.4 Program Operating Messages

When the program is initiated, the progress through the code is indicated on

the screen. The following sequence will occur on the CRT indicating the program

is progressing in a satisfactory manner.

* FLANK

- IN SOLVE 2 - AT BEGINNING

- IN SOLVE 2 - CALLING ASEMB2

- IN SOLVE 2 - ASEMB2 DONE, CALLING FORM B - IN SOLVE 2 - DONE FORM B, CALLING APPTBC

- DONE APPTBC, CALLING ESQSOL2(0)

- FLOATING POINT - UNDERFLOW ERROR AT D70D2

- NEXT INSTRUCTION AT D70D8

This indicates that EQS0L2 is being executed which begins the matrix decomposition.

The underflows in this section can be disregarded.

A-31
4.4.5 Time Step Indication

After the matrix has been decomposed and the temperature distributions at the

various nodes are going to be calculated , the message

- BEGIN TIME STEF NO. 4 AT TIME 0.27000E + 04

occurs
5. EXAMPLE PROBLEMS

The sample problems that follow are provided so that the user can see the input

data as well as the output. In both cases, the input file is displayed but

the complete output is displayed only for the static case since the output for

the dynamic problem consists of about 200 computer pages. The output is obtained

by typing PRINT FLANK. OUT on the PE. The output file - name is FLANK. OUT.

5.1 A STATIC PROBLEM

The input for this is shown as indicated on page A-90. Note that in the second

line of FLANK. INP, the next to the last element is 0 (which sets ITRAN to 0).

The output indicates that the net metering chamber energy is -114.7186 watts.

The minus sign indicates that heat is flowing into the chamber which is reasonable since the MC temperature is 38° C whereas the EC temperature is 65° C. The

heat transferred through the 1-D wall is 110.7266 watts. The two different

frame coefficients are 2.62 percent and 3.48 percent.

5.2 A DYNAMIC PROBLEM

The input for the dynamic problem is also shown on page A-100. Notice that

NUMIT.is - 1 indicating that a uniform initial temperature distribution given

by TEMP0 is used. Line 16 provides this initial temperature (25°C). IPR is

set equal to 4 so that the temperatures are printed out every 4th time step or

in this case each hour. IQPR is the energy tally printout which is provided

every two hours in this example since the time step is 15 minutes. After 24

hours, the frame coefficients were 0.48 percent and 0.95 percent.

A-33
6 REFERENCES

1. Lavine, A.G., Rucker, J.L., and Wilkes, K.E., "Flanking Loss Calibration

for a Calibrated Hot Box," Owens-Corning Fiberglass, Technical Center,

Granville, Ohio, (1981 or 1982).

2. Nakamura, S., "Computational Methods in Engineering and Science," Published

by John Wiley & Sons, 1977, p. 143.

A- 34
7. CODE LISTING

A listing of the code is provided so that definitions of various terms etc.,

may be obtained directly from the FORTRAN statements if desired. The code

consists of a main section which controls the sequence of steps through the

code which is determined by the input data.

The line numbers, right arrows and the message "! RANGE 2 NOT FOUND" are output from the editor and may be ignored.

A-35

o

n

o

ooo

ooo

C FLANK. FTN

C FLANKING LOSS CALCULATION FOR CALIBRATED HOT BOX

C IMPLICIT TIME INTEGRATION SCHEME, ALFA IS GREATER THAN ZERO

C AND LESS THAN OR EQUAL TO ONE LINEAR NODAL EQUATIONS, ALL COEFFICIENTS ARE CONSTANT

CODED BY ONEGA AND BURNS ON 7/6/82

MATERIAL PROPERTIES: MATERIAL NUMBER

12

34

5

6

MATERIAL TYPE

POLYURETHANE FOAM

FRP NEOPRENE RUBER GASKET

INFLATABLE RUBBER GASKET

WALL MATERIAL

NULL MATERIAL

COMMON /HEADG/ HEAD (16), TIMEM , DATEM

COMMON/CONTRL/NUMNP, NUMIT , NUMNPL , NUMTBC , IPR , IDATA , TEMPO

COMMON /NUMZ/ NZPOLY, NZFRP, NZWALL

COMMON /NUMY/ NY5,NY4,NY1L,NYWALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTY4,NTY1L,NTYWAL,NTY1R,NTY7,NCOL

COMMON /GRIDBZ/ NTZFRP, NTZWAL COMMON /ZCOOR/ ZPOLY.ZFRP.ZWALL

COMMON /YCOOR/ Y5 ,Y4.Y1L, YWALL, Y1R.Y7.Y8

COMMON /MATS/ RKAY (6) ,RHO(6) ,CP(6)

COMMON /HTCOEF/ HTC(8) ,TAMB(8) ,TAMBP(8)

COMMON /TEMPBC/ NODEST( 100) ,TBC( 100)

COMMON /TIM/ TIMO, TIMF, DELTA, ALFA, IRFINE

COMMON /DBUG/ IDBUG(IO)

COMMON /HEAT/ Q(21 ) ,IEB, IQPR,ICHECK(5)

COMMON /FIXT/ ITRAN COMMON /PLOTS/ IPLT.NUMCON,NUMHIS,IOPT,NOUT,RANGE(4),NSTEPH,

• ICHAR (10), NSTEPC

COMMON /PCONT/ ITITC(5,72) , ITITXC ( 5 , 36 ) ,ITITYC(5,36) ,NC(5)

,

• NODESC (5,20)

COMMON /PHIST/ ITITH(5,72) , ITITXH ( 5 , 36 ) ,ITITYH(5,36) ,MH(5)

,

• NODESH ( 5 , 5 )

, N , XMIN , XMAX

C DIMENSION C(246 ,246) ,A(246 ,246) ,UL(246 ,246)

DIMENSION B ( 246 )

, T ( 246 )

, DT ( 246 )

, EQUIL ( 246 ) , IPVT ( 246 )

,

• DX( 246 ) , RES( 246 )

, BTEMP ( 246

)

NNODES=246

DATA INPUT AND ECHO

CALL INPUT2 (C, A, UL,B,T,DT, EQUIL, IPVT, DX, RES, NNODES.IDGT) IF (IDATA. NE.O) CALL ADIOS (1)

SOLUTION PHASE

CALL S0LVE2 (C, A, UL,B, BTEMP, T.DT, EQUIL, IPVT. DX, RES, NNODES,

• IDGT)

WRITE TEMPERATURE FILE

IF (NUMIT. LE. -2) WRITE (3,1001) (T(J) . J=1 ,NUMNP)

C CALL ADIOS (1) C

1001 FORMAT (5E14.6) C

END

A-36
cCIIIIIIIIIIHHIHIIIHHIHIIIIHHIIINIIIIIIIIIIIIIIIIIIIIIIIHII

c SUBROUTINE ADIOS (N) CCHIIIMIIIIIIIIIIIIIIIIIIIIMIHHMIIIHIHIIHIIIIIIIIIIIIIIIIII

C* SOBRODTINE TO PERFORM EUT FUNCTIONS FOR PROGRAM

CIIIIIIIIIHIIIIIHIIIIIIIIIIIIIMIMHIimmilllllllllHIIIIHII

c

CALL HEADER

1

IF (N.NE.2) GO TO 100 WRITE (6,3000) GO TO 200

100 WRITE (6,3001)

200 CONTINUE

C

STOP c

3000 FORMAT (///20X,» ERROR TERMINATION *//)

3001 FORMAT (///20X, NORMAL TERMINATION

• ///) C

END

C SOBRODTINE HEADER 1 CCIHIIHIItIHIIItlllllHHIimillHIIIIHIIHIIIIIHIIHIItllllll

C* SUBROUTINE TO PERFORM PAGING AND WRITE HEADER • CMIIIimHIHIIIHIIIIIIIIHHIIIIIIIIIIIIHIIIIIIIIIHIIIHIIIII

c

COMMON /HEADG/ HEAD( 16) .TIMEM.DATEM

DATA STAR/4H****/

C WRITE (6,2000) (STAR, Id, 20) WRITE (6,2001) (HEAD(I),Id,l6)

WRITE (6,2002) (STAR, Id ,20) C RETURN

C

2000 FORMAT (1H1//20A4/)

2001 FORMAT (16A4)

2002 FORMAT (/20A4//)

C

END

A-37
ooooo

ooooo

o

ooooo

ooo

SUBROUTINE INPUT2 (C, A,UL,B,T,DT,EQUIL,IPVT,DX, RES, N.IDGT)

SUBROUTINE TO PREFORM DATA INPUT, INITIALIZATION AND ECHO •

COMMON/CONTRL/NUMNP, NUMIT , NUMNPL , NUMTBC , IPR , IDATA , TEMPO

COMMON /HEADG/ HEAD (16) ,TIMEM,DATEM

COMMON /NUMZ/ NZPOLY , NZFRP, NZWALL

COMMON /NUMY/ NY5,NY4,NY1L,NYWALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTYH,NTY1L,NTYWAL,NTY1R,NTY7,NCOL

COMMON /GRIDBZ/ NTZFRP, NTZWAL COMMON /ZCOOR/ ZPOLY , ZFRP , ZWALL COMMON /YCOOR/ Y5,Y4,Y1L,YWALL,Y1R,Y7,Y8

COMMON /MATS/ RKAY ( 6 )

, RHO ( 6 ), CP ( 6

)

COMMON /HTCOEF/ HTC(8) ,TAMB(8) ,TAMBP(8)

COMMON /TEMPBC/ NODEST( 100) ,TBC( 100)

COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

COMMON /DBUG/ IDBUG(IO)

COMMON /HEAT/ Q(21 ) ,IEB,IQPR,ICHECK(5)

COMMON /FIXT/ ITRAN

COMMON /PLOTS/ IPLT , NUMCON , NUMHIS , IOPT , NOUT , RANGE ( 4 )

, NSTEPH

,

• ICHAR(IO), NSTEPC

DIMENSION C(N,N),A(N,N),UL(N,N)

DIMENSION B(N) ,T(N) ,DT(N) ,EQUIL(N) ,IPVT(N) ,DX(N) ,RES(N)

READ (5,1000) (HEAD(I),I=1,16)

CALL GE0M2 (IDGT)

CALL GRIDBY

CALL GRIDBZ

CALL COEF

DO 10 K=1 ,8 10 TAMBP(K)=TAMB(K)

CALL MAT

IF (NUMTBC. NE.O) CALL TEMPB

TRANSIENT INPUT

DELTAs 1. ALFA=1.

TIM0=0. TIMF=DELTA

IF (ITRAN. NE.O) CALL TIMES

C IF (NUMIT. NE.O) CALL INIT (T.N, ITRAN)

C IF (IPLT. NE. 0) CALL PLOTIN (ITRAN)

C

IF (NUMNP.GT.N )GO TO 100 CALL HEADER

WRITE (6,2000) RETURN

, 100 WRITE (6,3000) NUMNP.N

CALL ADIOS (2)

C 1000 FORMAT ( 16AH)

2000 FORMAT (//20X," INPUT PHASE COMPLETE "//)

3000 FORMAT (////20X, "NO. OF NODES POINTS =",110//

• / 20X, "EXCEEDS DIMED VALUE =",110//

• / 20X, "STOP .") C END A-38
c

c SUBROUTINE CE0H2 (IDCT)

CC«» c»CM

c

SUBROUTINE TO INPUT GRID

COt«ON/CONTRL/NUMNP , NUKI T , NUMN PL , NUKTBC , I PR , IDATA , TEMPO

COfMON /NUMZ/ NZPOLT , NZFRP , NZWALL COMMON /NUMT/ NTS , NT4 , NY 1L , NTWALL , NY 1R, WY7 , NY8 COMMON /GRIDBT/ NTY4 , NTY 1L, NTYWAL, NTT 1R, NTY7 , NCOL COMMON /GRIDBZ/ NTZFRP, NTZWAL

COMMON /ZCOOR/ ZPOLT.ZFRP.ZWALL

COMMON /TCOOR/ T5.T4.T1L,TWALL,T1R,T7,T8

COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

COMMON /DBUG/ IDBUGOO)

C0»*10N /HEAT/ Q( 21 ), IEB, IQPR , ICHECK (5) COMMON /FIIT/ ITRAN COmON /PLOTS/ IPLT, NUMCON, NUMHIS, IOPT.NOUT , RANGE ( A ) , NSTEPH

,

• ICHAR ( 10) , NSTEPC

READ (5,1000) NUMTBC, NUMIT, IDATA, IPR, IRFINE, IDCT, IEB, IQPR,

• ITRAN, I PL

T

READ (5,1000) (IDBUG(I),I*1,10)

IF (IEB. HE. 0) READ (5.1000) (ICHECKI ) ,1* 1 ,5) READ (5.1000) NT5 ,NT4 , NT 1L, NTHALL, NT 1R,NT7,HT8

READ(5, 1001 ) Y5,Y4,Y1L, YWALL,T1R,T7,Y8

READ (5, 1000) NZPOLT, NZFRP, NZWALL

READ (5, 1001) ZPOLT.ZFRP.ZWALL

CALL HEADER 1 WRITE (6,2000) NUMTBC, NUMIT, IDATA, IPR, IRFINE, IDGT, IEB, IQPR,

• ITRAN, IPLT WRITE (6,3000) (IDBUG(I), 1x1,10)

IF (IEB. NE. 0) WRITE (6,3001) (ICHECEd ),Ix1 ,5) CALL HEADER 1 WRITE (6,4000) WRITE (6,2001) NT5 , NT4 , NT 1L, NTWALL, NT 1R, NT7 , NT8 WRITE (6,2002) T5,T4,T1L,TWALL,T1R,I7.I8

3000

3001

4000

2001

2002

2003

2004

WRITE (6,2003) NZPOLT, NZFRP, NZWALL WRITE (6,2004) ZPOLT.ZFRP.ZWALL

FORMAT (1615) FORMAT (8E10.0) FORMAT (//

/ 20Z , "C ONTROL DAT A"//

/ 20X, "NUKTBC «",I5/

/ 20X, "NUMIT •",15/

/ 201, "IDATA «",I5/

/ 20X, "IPR *",15/

/ 20X, "IRFINE *",15/

/ 20X, "IDGT *",15/

/ 20X, "IEB *",15/

/ 20X, "IQPR «",I5/

/ 20X, "ITRAN «*,I5/

/ 20X, "IPLT *",I5///) FORMAT (//

/ 20X, "DEBUG PARAMETERS"//

/ 15X," 1 (WRITE A.C MATRICES, B VECTOR) x",I5/

/ 15X," 2 (WRITE DECOMPOSED C MATRIX ) «",I5/

/ 15X," 3 (WRITE NODAL ASSOC. ) *",15/

/ 15X, 4 (WRITE ELEMENT MATRICES ) *",15/

/ 15X," 5 (WRITE INTERIOR DATA ) -M5/

/ 15X, " 6 (WRITE C BEFORE DECOMP. ) *",15/

/ 15X, 7 (WRITE B VECTOR ) "",15/

/ 15X," 8 ( ) *",15/

/ 15X, " 9 ( ) *",15/

/ 15X, • 10( ) *",I5//) FORMAT (//20X, "E NERGT BALAN C E S"/

/ 15X, "ICHECE ( 1

)

*",15/

/ 15X,"ICHECX(2) *",15/

/ 15X,"ICHECI(3) *",15/

/ 15X,"ICHECK(4) *",15/

/ 15X, "ICHECKI 5) *’,I5) FORMAT (//20X, "G RID SPECIFI C A T I 0 N"//) PORMAT (//20X , "I NUMBERS FOLLOW"

/ 20X, "NT5 *",15/

/ 201, "NY

4

*",15/

/ 20X, "NT 1L *",15/

/ 20X, "NTWALL *",15/

/ 20X , "NT 1R *",15/

/ 20X,"NT7 *",15/

/ 20X,"NT8 «",I5//) FORMAT (/20X,"T'S FOLLOW"/

/ 20X , "15 *",E1 1. 4/

/ 20X, "T4 *",E11.4/

/ 20X, "TIL *",E11.4/

/ 201, "TWALL «",E11.4/

/ 20X, "I1R «",E11.'4/

/ 20X, "T7 «",E11.4/

/ 20X, "T8 «",B11.4//) FORMAT (/20X, *Z NUMBERS FOLLOW"/

/ 20X, "NZPOLT ",15/

/ 20X, "NZFRP ",15/

/ 20X, "NZWALL •*,I5//) FORMAT (/20X,"Z'S FOLLOW"/

/ 20X, "ZPOLI *",E1 1. 4/

/ 20X, "ZFRP A — "3Q

/ 20X , "ZWALL

*",E1 1. 4/ * ,E1 1. 4//)

RETURN
o

o

o

o

o

o

o

SUBROUTINE GRIDBY

SUBROUTINE TO SPECIFY ALL THE Y BOURDARY NODE NUMBERS •

COMMON /NUMZ/ NZPOLY,NZFRP,NZWALL

COMMON /NUMY/ NY5,NYit,NY1L,NYWALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTY4,NTY1L,NTYWAL,NTY1R,NTY7,NCOL

COMMON /GRIDBZ/ NTZFRP, NTZWAL

COMMON /ZCOOR/ ZPOLY, ZFRP, ZWALL COMMON /YCOOR/ Y5, YH,Y1L,YWALL,Y1R,Y7,Y8

C NTY4=NY5+NY4

NTY1L=NTYH+NY1L

NTYWAL=NTY 1L+NYWALL

NTY1R=NTYWAL+NY 1R NTY7 =NTY 1 R+NY7

NCOL=NTY7+NY8

C RETURN

C END

CIIIIIIHIIHIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIHIIIMfllllHIII

C SUBROUTINE GRIDBZ

CCIIHKIlKimHHIIIIIIIlfMHIimitHllllfHHIHIIHIIIII

C« SUBROUTINE TO SPECIFY BOUNDARY NODES IN THE Z DIRECTION

C COMMON/CONTRL/NUMNP , NUMIT , NUMNPL , NUMTBC , IPR , IDATA , TEMPO

COMMON /NUMZ/ NZPOLY , NZFRP , NZWALL

COMMON /NUMY/ NY5,NY4,NY1L,NYVALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTY^NTY1L,NTYWAL,NTY1R,NTY7,NCOL

COMMON /GRIDBZ/ NTZFRP, NTZWAL COMMON /ZCOOR/ ZPOLY, ZFRP, ZWALL

COMMON /YCOOR/ Y5,Y4,Y1L,YWALL,Y1R,Y7,Y8

C NTZFRP=NZPOLY+NZFRP

NTZWAL =NTZFRP+NZWALL NUMNP=NCOL#NTZFRP+(NYWALL+1 )*NZWALL

NUMNPL=NTZFRP*NCOL

C CALL HEADER 1 WRITE (6,2000) WRITE (6,2001) NTYJ»,NTY1L,NTYWAL,NTY1R,NTY7,NCOL WRITE (6,2002) NTZFRP, NTZWAL WRITE (6,2003) NUMNP, NUMNPL

C RETURN

C2000 FORMAT ( • / 25X,"C UMULATIVE NODE NUMBER S"//)

2001 FORMAT (/20X,"Y - DIRECTION NODE NUMBERS"/

• / 20X, "NTY4 =",15/

• / 20X,"NTY1L =",15/

• / 20X, "NTYWAL =",15/

• / 20X, "NTY1R =",15/

• / 20X, "NTY7 =",15/

• / 20X, "NCOL =",I5)

2002 FORMAT (/20X,"Z - DIRECTION NODE NUMBERS"/

• / 20X, "NTZFRP =",15/

• / 20X, "NTZWAL =",I5/)

2003 FORMAT (/20X, "TOTAL NUMBER OF NODES (EQUATIONS) =",110/

• / 20X, "LOWER NUMBER OF NODE POINTS =",I10///)

C END

A-40
cc

cC*1 c«

C*' c

c

c

cc

c

c

50

SUBROUTINE INIT (T.N.ITRAN)

SUBROUTINE TO READ IN AND ECHO INITIAL TEMPERATURES

COMMON /CONTRL/NUMN P , NUMIT , NUMN PL , NUMTBC , I PR , IDATA , TEMPO

DIMENSION T(N)

IF (ITRAN.EQ. 0) GO TO 50

IF (NUMIT. GT.O) GO TO 200 READ (5,1001) TEMPO

CORRECT FOR STEADY CALCS.

IF (ITRAN.EQ. 0) TEMP0=0.

DO 100 1=1,

N

100 T(I)=TEMPO

C

CALL HEADER 1 WRITE (6,2000) TEMPO C

RETURN C

200 CONTINUE

IF (NUMIT. GE. 2) GO TO 300 READ (5.1003) (NODE, T(NODE), 1=1 .NUMNP) GO TO 400

300 READ (1,1002) (T(I), 1=1, NUMNP) C

400 CALL HEADER 1 WRITE (6,2000) WRITE (6,2001) (T(N0DE),N0DE=1, NUMNP)

C

RETURN C

1001 FORMAT (8E10.0)

1002 FORMAT (5E14.0)

1003 FORMAT (5(15, E10. 0)

)

2000 FORMAT (///20X,*I NITIAL TEMPERATURE S*//

• / 15X, "INITIAL TEMPERATURE =*,E14.5)

2001 FORMAT (8(E10. 0,4X)

)

C

END CHKIIMHHIHIHIIIIIIIHIIIIIIIIimailllllllllllllllllll

C

SUBROUTINE MAT

CClt*tl»HIIII«llllllllltllllllllf«IIIIIHI««*HHHIIIIIHt*<

C» SUBROUTINE TO INPUT AND ECHO MATERIAL PROPERTIES CtHHillHIlHHIlHIIHIHHIHKKHIIIIHIHIHHHIIIII

C

COMMON /MATS/ RKAY ( 6 ), RHO ( 6 ), CP ( 6

)

C

READ (5,1001) (RKAY(I ) ,1=1,5) READ (5.1001) (RHO(I),I=1,5)

READ (5,1001) (CPU), 1 = 1, 5) RKAY(6)=0.

RH0(6)=0.

CP(6)=0. C

CALL HEADER 1 WRITE (6,2000)

DO 100 J=1,6

100 WRITE (6,2001) J,RKAY(J),RHO(J),CP(J)

C

DO 200 J=1 ,5 200 RHO( J)=RH0(J )*CP( J)

C

RETURN C

1001 FORMAT (8E10.0)

2000 FORMAT ( //25X,*M ATERIAL PROPERTIES*

• /// 17X,*I*,10X,*RXAY»,13X,*RH0*,14X,*CP*//)

2001 FORMAT (15X,I5,3(5X,E11.4))

C END A-41
on

oo

o

o

o

o

oo

o

o

o

o

o

I C SUBROUTINE TEMPB

SUBROUTINE TO READ TEMPERATURE BCS

COMMON/CONTRL/NUMNP , NUMIT, NUMNPL , NUMTBC , IPR , IDATA , TEMPO

COMMON /TEMPBC/ NODEST( 1 00) ,TBC( 100)

READ (5,1003) (NODEST(I),TBC(I), 1=1, NUMTBC)

CALL HEADER 1 WRITE (6,2000) NUMTBC WRITE (6,2001) (NODEST(I),TBC(I), 1=1, NUMTBC)

RETURN

1003 FORMAT (5(I5,E10. 0)

)

2000 FORMAT (1H1

• / 25X, "T EMPERATURE BOUNDARY C 0 N D"

• / 40X, "NUMBER OF TEMPERATURE BCS =",I5//)

2001 FORMAT (5( 3X, 15 ,2X,E10. 3)

)

END

SUBROUTINE TIMES

C« SUBROUTINE TO INPUT AND ECHO TIMING INFORMATION • CttHIIHmitHIHfllllllHtlfHHHIimillllllllHIlHHIIIIIIIH

C COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

C READ (5,1001) TIMO.TIMF, DELTA, ALFA

C CALL HEADER 1 WRITE (6,2000) WRITE (6,2001) TIMO.TIMF, DELTA, ALFA

C RETURN

1001 FORMAT (8E10.0)

2000 FORMAT (//25X,"T I M I N G 2001 FORMAT (/20X, "TIMES, ETC."

• / 20X, "TIMO

• / 20X, "TIMF

• / 20X, "DELTA

• / 20X, "ALFA

PARAMETER S"//)

=",E1 1. 4/

=",E1 1. 4/

=",E11.4/

=",E1 1. 4/)

END

A- A 2
cCHIIHIHHmiHHimHIIIHHmHmiMHIlimilllHHIMHHa

c SUBROUTINE S0LVE2 (0, A , UL, B, BTEMP.T, DT, EQUIL, IPVT.DX, RES, N,

• IDOT)

CCHHIHIIIHIHHIlHimilHHUHIIHIIIIIIHIHUHmiimMIH

C» SUBROUTINE TO PERFORM TIKE INTEGRATION AND EBS • CHHIHHtllimmiHHHHHHHHHHHMHMIlHHIUmilUm

C COMMON /CONTRL/NDMNP , NOKIT , NUKNPL , NOKTBC , I PR , IDATA , TEMPO comoN

COMMON

COt«ON

COWiON

COMMON

COMMON

COMMON

COMMON COMMON

COMMON

COMMON

COMMON COMMON

COHMON

C01*K)N

tCOMtON

•C0M40N

• NODESH ( 5 , 5 ) , NH , ZMIN , ZMAZ C

DIMENSION C(N,N),A(N,N),UL(N,N)

DIMENSION B(N) , BTEMP(N) ,T(N) ,DT(N) ,EQOIL(N) ,IPVT(N) ,DX(N)

• RES(N) C

DATA BIG /1.E30/ WRITE (7,7000) 7000 FORMAT ( IN SOLVE2 - AT BEGINNING")

C

C ASSEMBLE A AND C MATRICES AND THE B VECTOR

C

DO 100 J=1,N

BTEMP(J)=0.0

B(J)sO. DT(J)=0. DO 100 1=1,

C(I,J)=0.

100 A(I,J)=0. WRITE (7,7001)

7001 FORMAT ( IN S0LVE2 - CALLING ASEMB2")

CALL ASEMB2 (C,A,N) WRITE (7,7002) 7002 FORMAT (* IN S0LVE2 - ASEMB2 DONE, CALLING FORMB")

CALL FORMB (0,B,BTEMP,N)

C

C FORM AND REDUCE SOLUTION MATRIX C

CDIVsI.

IF (ITRAN.EQ.O) CDIV=1.E30

DO 200 J = 1 ,NUMNP

DO 200 1=1 ,NUMNP 200 C(I,J)=C(I,J)/CDIV A(I, J) C

C ADD TEMP BCS FOR STEADY STATE PROBLEM

C IF (NUMTBC.EQ. 0) GO TO 220

DO 210 1=1 ,NUMTBC JeNODEST(I)

210 C(J,J)=BIG WRITE (7,7003) 7003 FORMAT (• IN SOLVE2 - DONE FORMB, CALLING APPTBC")

IF (ITRAN.EQ.O) CALL APPTBC (NUHTBC,ITRAN,BIG,B,T,N)

220 CONTINUE

C

C WRITE OUT DBUG MATRICES C

IF (IDBOG(I).EQ.O) 00 TO 260 WRITE (6,3000)

DO 251 1=1 ,NUMNF WRITE (6,4001) I 251 WRITE (6,4000) (A(I, J) , J=1 ,NUMNP) WRITE (6,3001)

DO 255 I*1,HUMNP WRITE (6,4001) I 255 WRITE (6,4000) (C(I, J) , J=1,NUMNP) WRITE (6,3002) WRITE (6,4000) (B(J),J=1 .NUMNP) IF (IDBUG(I).LT.O) CALL ADIOS (1) 260 CONTINUE

C

/HEADG/ HEAD( 16 ) ,TIMEM, DATEM

/NUM2/ NZPOLI, NZFRP, NZWALL

/NUMY/ NY5,NY4,NY1L,NYWALL,NY1R,NY7,NY6

/GRIDBY/ NTY4 ,NTY 1L, NTYWAL, NTY 1R, NTY7,NCOL

/GRIDBZ/ NTZFRP, NTZWAL

/ZCOOR/ ZPOLY, ZFRP, ZWALL /YCOOR/ Y5, Y4, Y1L, YWALL, Y1R, Y7, Y8 /MATS/ RKAY (6) , RHO(6) , CP ( 6

/HTCOEF/ HTC(8) ,TAMB(8) ,TAMBP( 8) /TEMPBC/ NODEST (100), TBC (100) /TIM/ TIMO , TIMF, DELTA , ALFA , IRFINE

/DBUG/ IDBUG(IO)

/HEAT/ Q(21 ), IEB, IQPR, ICHECK(5) /FIXT/ ITRAN

/PLOTS/ IPLT , NUMCON , NUHHIS , IOPT , NOOT , RANGE ( 4 ), NSTEPH

ICHAR ( 10) ,NSTEPC /PCONT/ ITITC(5,72),ITITXC(5,36),ITITYC(5,36),NC(5),

NODESC(5,20)

/PHIST/ ITITH(5,72),ITITXH(5,36),ITITYH(5,36),MH(5),

A-43
WRITE (7, 7004) 700* FORMAT (» IN S0LVE2 - DONE APTBC, CALLING EQS0L2(0)")

CALL EQS0L2 (O,C,A,0L,B,T,DT,EQUIL,IPVT,DX,RES,N,IDGT,

• TIME)

c

IF (IDBUG(2).LT.O) CALL ADIOS (1) C

C WRITE INITIAL STATE INTO PLOT FILE

C IF (IPLT.EQ. 0) GO TO 270

IF (NUMCON.NE. 0) CALL PLOTC (O.TIME.T.N)

IF (NUMHIS.NE.O) CALL PLOTH (0,TIME,T,N)

C

C ICOONT - COUNTER FOR PRINTING

C KOUNT - COUNTER FOR ENERGY BALANCE WRITES

C KOUNTC - COUNTER FOR CONTOUR PLOTS

C KOUNTH - COUNTER FOR TIME HISTORY PLOTS

C270 ISTEP=0 ICOUNT=0 KOUNTrO KOUNTC=0

KOUNTH=0 TIME=TIMO

C

C PRINT INITIAL TEMPERATURES

C CALL PRINT2 (TIME,T,N,NUMNP)

C

C TIME STEP BEGINS

C WRITE (7,7005) 7005 FORMAT (• IN S0LVE2 - DONE DECOMPOSITION , BEGIN TIME LOOP")

300 ISTEP=ISTEP+1 WRITE (7,7006) ISTEP.TIME

7006 FORMAT (" BEGIN TIME STEP NO. ",I5," AT TIME »,E12.5) IC0UNT=IC0UNT+1 K0UNT=K0UNT+1 ALTIM=TIME+ALFA "DELTA TIME=TIME+DELTA

C

C FORM RHS - B VECTOR

C

IF (ITRAN.GE.O) GO TO 350 READ (5,1001) (TBC(I) ,1=1 .NUMTBC

)

READ (5,1001) (TAMBP(I ) ,1=1 ,8) 350 IF (ITRAN.NE.O) CALL FORMB (ITRAN, B, BTEMP, N) IF (ITRAN.NE.O) CALL L0DE2 (B,T, A,N, ALFA) IF (ITRAN.NE.O) CALL APPTBC (NUMTBC, ITRAN. BIG, B,T,N)

C

C SOLUTION PHASE

C CALL EQS0L2 ( 1 ,C, A,UL,B,T,DT,EQUIL, IPVT,DX,RES,N,IDGT,

• TIME)

C

C UPDATE TEMPERATURES

C

CALL UPDATE (T.DT.N)

C

C PRINTING

C

IF (ICOUNT.LT. IPR) GO TO *00 CALL PRINT2 ( TIME, T,N,NUMNP) ICOUNT=0

C

C PERFORM ENERGY BALANCES AND PRINTING THEREOF

C*00 IF (IEB.EQ. 0) GO TO 500 CALL QSEMB2 (T.N) CALL QS0LV2 (ISTEP, KOUNT, T.N, TIME) e

C PLOTTING

C500 IF (IPLT.EQ. 0) GO TO 600

C

C CONTOURS

c

IF (KOUNTC. LT.NSTEPC) GO TO 550

IF (NUMCON.NE. 0) CALL PLOTC (0, TIME, T.N) KOUNTC=0

C

C TIME HISTORIES C

550 IF (KOUNTH. LT.BSTEPH) GO TO 600

IF (NUMHIS.NE.O) CALL PLOTH (0, TIME, T.N) KOUNTH=0

CC TEST FOR END OF TIME LOOP

C A-44
600 IF (TIME. LT. TIMF) 00 TO 300 C

C FINAL PLOTS c

IF (NUMCON. HE. 0) CALL PLOTC (I.TIME.T.N)

IF (NUMHIS. NE. 0) CALL PLOTH (0,TIME,T,N)

IF (NUMHIS. NE. 0) CALL PLOTH (1,TIME,T,N) c

RETURN C

1001 PORMAT (8E10.0) 3000 FORMAT (////

• / 15X,*A MATRIX FOLLOWS*//)

3001 FORMAT (////

• / 15X, *C MATRIX FOLLOWS*//)

3002 FORMAT (////

• / 15X, *B VECTOR FOLLOWS*//)

000 FORMAT (10( IX, Ell. A)) 001 FORMAT ( // 25X,*R0W NUMBER -*I5//) c

END

A-45
cCMIIHIIHIIIIIIIHHIIIII<IHHIIHI<IIII<IIIIHIIIIMIIH

c

SUBROUTINE ACONV (N,M,DY,DZ,MM,A0,A1,A2,A3,A4,C0)

CC<HIIHHIHHIIHHHH<HIH< HHHHI<HHI<H<

C‘ SUBROUTINE TO DETERMINE A(IJ) FOR CONVECTION BOUNDARY C<HHHHHHHHHHHHH<HHHHHHH<HHHHHHHH

c

COMMON /MATS/ RKAY(6) , RHO(6) ,CP(6) COMMON /HTCOEF/ HTC(8) ,TAMB( 8) ,TAMBP(8) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

C CO=(RHO(MM)‘DY‘DZ / (2. ‘DELTA))

A01=HTC(M)‘DY

A 1 =-RKAY (MM )‘DZ/ ( 2. ‘DT

)

A2=A1 A3=-RKAY(MM)‘DY/DZ

A4sA3

C GO TO (1,2, 3, 4) ,N

1 A3=A2

A2=A4

A4sA1

A 1=0. GO TO 10 C

2 A1=A3

A3=A2

A4=A3

A2=0. GO TO 10 C

3 A3=0. GO TO 10 c

4 A4=0.

10 A0=-(AHA2+A3+A4)+A01

C

RETURN

C END

C“

c

cc*«

c«c»‘

c

c

c

c

SUBROUTINE AIC1 (MA,MB,MC,MD,DX1,DX2,DX3,DX4,

• A0,A1,A2,A3,A4,C0)

SUBROUTINE TO CALCULATE A(I,J) FOR INTERIOR CORNERS

COMMON /MATS/ RKAY(6) ,RH0(6) ,CP(6) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

C0=(RH0(MA)‘DX2*DX4 + RH0(MB)»DX4»DX1 * RH0(MC)»DX1»DX3

• RHO ( MD ) *DX 3*DX2 ) / (4. ‘DELTA)

A0=(RKAY(MA)‘(DX2/DX4+DX4/DX2) + RKAY(MB)‘(DX1/DX4+DX4/DX1)

• RKAY(MC)‘(DX3/DX1+DX1/DX3) RKAY(MD)‘(DX2/DX3+DX3/DX2)

)

• / 2. A1=-(RKAY (MB)‘DX4/DX1 + RKAY(MC)‘DX3/DX1) / 2. A2=- ( RKA Y ( MD ) *DX 3/DX2 + RKAY(MA)‘DX4/DX2) / 2. A3=-(RKAY(MC)‘DX1/DX3 + RKAY(MD)*DX2/DX3) / 2. A4=-( RKAY (MA )‘DX2/DX4 + RKAY(MB)‘DX1/DX4) / 2.

RETURN

END

A-46
eCNIIINIIIIIIIIIHIIHIIKIIIHHIIHIIIIMIIKIHIIIHIIIIIII

c SOBROOTINE AIC2 ( N 1 , N2 , WH 1

, NH2 , MA , KB , MC , KD

• DX1,DX2,DX3,DXA,A0,A1,A2,A3,AA,C0) ccmmiiimiiiiimiimmiHHiHimmimHimiHiimi

C‘ SOBROOTINE TO CALCOLATE A(IJ) POR REENTRANT CORNERS C“‘“

c

COMMON /MATS/ RXAK6) ,RH0(6) ,CP(6) COMMON /HTCOEP/ HTC (8) ,TAMB( 8) ,TAMBP(8)

C0ANA0N /TIM/ TINO,TIMF,DELTA p ALFA,IRPINE

C C0z(RH0(MA )*DXA*DX2 RH0(MB)‘DI1‘DXA *RH0(MC )‘DX 1‘DX3

• RH0(MD)‘DX3*DI2) / (A. ‘DELTA) DDXIaDXI

DDX2*DXA

IP (N1.EQ. 2) DDXUDX2

A0r0.5‘HTC(NH1)‘DDX1

A01«0.5‘HTC(NH2)‘DDX2

A1z-(RKA?(MB)*DXA/DX1 RIAT (MC)«DX 3/0X1) / 2. A2*-(RKAT(MD)»DX3/DX2 RKAT (MA )*DXA/DX2) / 2. A3*-(RXAT(MC)»DX 1/0X3 HKAT(MD)*DX2/DX3) / 2. AA*-(RKAT(MA)‘DX2/DXA RKAT(MB)‘DX1/DXA) / 2. C

A Or- ( A 1 A2*A 3+A A ) AO+A0

SOBROOTINE AIC3 (N,NH,NA,MB,DI,DXA,DXB,A0,A1 ,A2,A3,AA,C0) C

C‘ SOBROOTINE TO COMPOTE A(IJ)P0R COHNER NODE WITH CONVECTION • C‘ ON ONE SIDE • CHMKIIIimiHimtHIMMHIIIMIIItHIIIHIIIHHIIIMIMIIMIt

c

COAMOH /MATS/ RXAK6) ,RH0(6) ,CP(6)

COANION /HTCOEF/ HTC(8),TAMB(8),TAMBP(8)

C0AN10N /TIM/ TIMO , TIMF , DELTA , ALFA , IRPIHE C

C0r( RH0(MA )*DXA RHO(MB)*DXB)*DX/(A. ‘DELTA)

A01 rHTC (NH )‘0. 5‘ (DIA^DXB)

A 1 »-( RIAT (MA )‘DIA RKAT (MB) ‘DIB) / (2.‘DX) A2<A1 A3=-RXAT (MA )‘DX/(2. ‘DU ) AA«-RLAT (MB)‘DX/(2. ‘DIB)

00

TO (1,2, 3. A) ,N

1 AlrO. 00 TO 10 c

2 A2rA3 A3«AA AA«A2

A2r0. 00 TO 10 C

3 AlsAA AAcA2

A2xA3

A3*0. C

00 TO 10

A A1rA3

A3=A2

A2rAA

AA»0. C

10 AOa- ( A 1 +L2+t 3*A A ) +A0

C

RETURN C

END A-47
cC“

c

cc“

c»

c«CM

c

c

c

c

c

c

c

c

SUBROUTINE AIC4 (II ,N2,NH1 ,HH2,M,DI1,DX2,

• AO,A1,A2,A3,A4,CO)

SUBROUTINE TO CALCULATE A(I,J) FOR EXTERIOR CORNER WITH HEAT TRANSFER ON TWO SIDES

COMMON /MATS/ RXAT(6) , RH0(6) ,CP(6

)

COMMON /HTCOEF/ HTC(8) ,TAMB( 8) ,TAHBP(8) COMMON /TIM/ TIMO.TIMF .DELTA, ALFA, IBFINE

CO=RHO(M)*DX1 ‘DX2/ ( 4. ‘DELTA

)

A1=-RKAY(M)‘DX2/(2. «DX1) A2=A1 A3=-RKAY(M)‘DX1 / (2. *DX2) A4=A3

A0=0. 5*DX2*HTC ( NH 1

)

A01s0.5‘DX1‘HTC(NH2)

IF (N1.NE. 1) GO TO 2 IF (H2.NE.3) GO TO 1

A1=0.

A3=0. GO TO 10

1 AlsO.

A4s0. GO TO 10

2 IF (N2.NE.3) GO TO 3 A2s0.

A3=0. GO TO 10

3 A2s0.

A4s0.

10 A0sAO+A01-(A1+A2+A3+A4

)

RETURN

END

SUBROUTINE AINF (»«,MA,MB,DX,DXA,DIB t A0,A1 , A2,A3,A4,C0)

CCillimiltlltflKHIIItlHIItUlHIfllHIHHIlHIHHIIMII

C‘ SUBROUTINE TO CALCULATE A(IJ) AT THE INTERFACES

C

COMMON /MATS/ RKAY(6) ,RH0(6) ,CP(6) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA .IRFINE

C COs(RHO(MA)*DXA RHO(MB)‘DXB)‘DX / (2. ‘DELTA)

A 1 s- ( REAY (MA ) *DXA/ ( 2. *DX ) + RKAY ( MB ) ‘DXB/ ( 2. ‘DI )

)

A2sA1 A3=-REAY(MB)‘DX/DXB

A4s-RXAI(MA)‘DI/DXA

AOs- ( A 1 +A2+A 3+A 4

)

C

IF (MLEQ. 1) RETURN

A1sA4

A4sA2

A2sA3 A3>A4 C

RETURN C

END

48

t
e

SUBROUTINE AINT (DY,DZ,M, A0,A1 ,A2,A3,AI|,C0)

C

C* ROUTINE TO CALCULATE A(IJ) FOR GENERAL INTERIOR NODE • CHHimimiHUIIHMHIHimHHHimiimHHIlHIHIMIHIII

c

COMMON /DBUG/ IDBOG(IO)

C COMMON /MATS/ REAY(6) ,RH0(6) ,CP(6)

C0AN10N /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE C

CO* ( RHO ( M )»DY»DZ/DELTA

)

AObRKAY (M)*2. •(dy/dz*dz/dt)

A 1 =-RIAY (M)*DZ/DY A2.A1 A3*-RKAY(M)»DY/DZ

A4cA3 C

IF (IDB0G(5).EQ. 0) RETURN WRITE (6,3000) M,REAY(M),C0,A0,A1,A2,A3,AA

3000 FORMAT (5X,"M,E,C0,A(5) »*,I5,7(2X,E12.5)) c

END

C

c

cCM*

c«

cc

c

100 c

c

c 3000

c

SUBROUTINE AMATRZ (A, N, NN, AN

)

SUBROUTINE TO PERFORM ALGEBRAIC ASSEMBLY OF THE A MATRIX

COMMON /DBUG/ IDBUG(IO)

DIMENSION NN(5) ,AN(5) DIMENSION A(N,N)

DO 100 K*1 ,5 A(NN(1),NN(K))*A(NN(1) ,NN(E ) )*AN(E)

IF (IDBUG(A).EQ. 0) RETURN WRITE (6,3000) (NN(L),L=1,5),(AN(L),L*1,5)

RETURN

FORMAT (51, "IN AMATRX - NN(5),AN(5) «",2I,5I5,5(2I,B12.5))

m>

i

cc*» c

c

C»«' c»

c

c

c

c

100 c

c

SUBROUTINE APPTBC ( NUMTBC , ITRAN , BIG , B , T , NN

)

SUBROUITNE TO APPLY TEMPERATURE BCS VIA PENALTY FUNCTION

COMMON /TEMPBC/ NODESTt 100) ,TBC( 100)

DIMENSION B(NN) ,T(NN)

IF (NUMTBC. EQ.O) RETURN

DO 100 1*1, NUMTBC J*N0DEST(I)

B(J)=BIG»(TBC(I)-T(J))

RETURN m>

A-49
ce

c

SUBROUTINE EQSOL2 (K , C, A, UL,B, T, DT, EQUIL, IPVT.DX, RES, N, IDGT,

• TIME)

CCHimHIIIIIIIIIMHIlHHIIIIIIHIIIIIIIHIIHIHMHIIMIIHIII

C SUBROUTINE TO PERPOHM MATRIX CALCULATIONS CMIIHIIHIIHIIIHIIIHIHHIIIIIHIimHIHHIIIHHIIIHHIIII

c

COMMON /CONTRL/ NUMNP, NUMIT, NUMNPL, NUMTBC, IPR, IDATA, TEMPO COMMON /TIM/ TIMO.TIMF, DELTA , ALFA, IRFINE

COMMON /DBUG/ IDBUG(IO) C

DIMENSION C(N,N),A(N,N),UL(N,N)

DIMENSION B(N) ,T(N) ,DT(N) ,BQUIL(N) ,IPVT(N) ,DX(N) ,RES(N) C

IF (X) 200,100,200

c

C LU DECOMPOSITION

C

100 CONTINUE

IA=N

NN= NUMNP

IF (IDBUG(6).EQ.O) GO TO 120 WRITE (6,3500) DO 110 1=1, NUMNP WRITE (6,3501) I 110 WRITE (6,3502) (C(I, J) ,J=1 , NUMNP)

IF (IDBUG(6).LT. 0) CALL ADIOS (1) 120 CONTINUE

CALL LUDATF (C, UL,NN,IA, IDGT, D1,D2,IPVT, EQUIL, WA, IER) D=D1*2.

#D2 CALL HEADER 1 WRITE (6,2000) D,WA,NN,IA,IER

IF (IER.EQ.O) GO TO 150 WRITE (6,3000) IER

IF (IER. GE. 128) CALL ADIOS (2) 150 IF ( IDBUG ( 2 ) . EQ. 0 ) RETURN

CALL HEADER 1 WRITE (6,4000)

C DO 170 1=1, NUMNP WRITE (6,3501) I 170 WRITE (6,3502) (C(I,J),J=1, NUMNP)

IF (IDBUG( 2). LT. 0) CALL ADIOS (1) C

RETURN C

C SOLUTION PHASE

C 200 CONTINUE

CALL LUELMF (UL,B,IPVT,NN,IA,DT)

C

IF (IDBUG(7).EQ. 0) GO TO 300

C

C WRITE SOLUTION VECTOR

C WRITE (6,5000) TIME WRITE (6,3502) (B( J) ,J=1 .NUMNP) 300 CONTINUE

C

C REFINEMENT C

IF (IRFINE. EQ.O) RETURN

C CALL LUREFF (C,B,UL,IPVT,NN,IA,DT,IDGT,RES,DX,IEH)

IF (IER. EO. 0) RETURN WRITE (6,3001) IER CALL ADIOS (2) C

2000 FORMAT (//

• / 15X,"E QUATION DAT A"//

• / 15X , "DETERMINANT =",E14.7/

• / 15X, "ACCURACY TEST PARAMETER =",E14.7/

• / 15X, "NUMBER OF EQUATIONS =",110/

• / 15X, "ROW SIZE IN CALLING ROUTINE =",110/

• / 15X, "ERROR PARAMETER ( 129-SIN. , 34-ACC. ) =",I10///)

3000 FORMAT (//

• / 15X, " ERROR FROM EQUATION SOLVER /

• / 20X, "DECOMPOSITION PHASE"/

• / 20X, "IER =",I10//)

3001 FORMAT (//

• / 15X, " ERROR FROM EQUATION SOLVER •/

• / 20X, "REFINEMENT PHASE"/

• / 20X, "IER =",I10//)

3500 FORMAT (//

• 15X,"IN EQSOL2 C MATRIX BEFORE DECOMPOSITION FOLLOWS"//

• ) 3501 FORMAT (//25X,"R OW NUMBER =",I10//)

3502 FORMAT ( 10(2X,E1 1. 4)

)

4000 FORMAT (//15X,"IN EQSOL2 «/

• / 15X, "DECOMPOSED CAPACITANCE MATRIX FOLLOWS"///)

5000 FORMAT (//

• / 25X, " IN EQSOL2 AT TIME =",E14.7//

• / 30X,»R H S V E C T 0 R"//)

C A-50

END
ClllltHIIHKItlllHIlHlllltHtHHIIIIHIimiHIIIHIHHimill

c SUBROUTINE GETNOD (I,J,N) CCHIIIIHIIIimillHIIIIIHIHIIIimHHHHIIIIHIIIIIIIHIIIIIII

C» SUBROUTINE TO MAP FROM I,J TO NODE NUMBER N(5) • CIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIHIIIIIIIIIIHIIIIIIIIHIIIIIIIIII

C

COHMON/CONTRL/NUMNP , NUMIT , NUMNPL , NUMTBC , IPR . IDATA , TEMPO COMMON /NUMZ/ NZ POLY, NZFRP, NZWALL COMMON /NUMY/ NY5,NY4,NY1L,NYWALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTY , NTY 1L , NTYWAL , NTY 1 R , NTY7 , NCOL COMMON /GRIDBZ/ NTZFRP, NTZWAL COMMON /ZCOOR/ ZPOLY, ZFRP, ZWALL COMMON /YCOOR/ Y5,Y4,Y1L,YWALL,Y1R,Y7,Y8

COMMON /DBUG/ IDBUG(IO) C

DIMENSION N(5)

C

IF (I.GT.NTZFRP+1) GO TO 300 IF (I. GE. NTZFRP) GO TO 200 C

C LOVER REGION

C N(1)=(I-1)*NCOL«J

N(2)=N( 1 )-1

N(3)=N(1)f1

N(4)=N( 1)-NCOL

N(5)=N(1)+NCOL

N ( 2 ) =MAX 0 ( N ( 2 ), 1

)

N(4)=MAX0(N(4) , 1

)

N(3)=MIN0(N(3), NUMNPL)

N(5)=MIN0(N(5), NUMNPL)

C GO TO 400 C

C MIDDLE REGION

C 200 IF (I. GT. NTZFRP) GO TO 250

C

C ON TOP OF BASE

C N(1)=(I-1)«NCOL+J

N(2)rN(1)-1

N(3)=N( 1 )+1

N(4)=N( 1 )-NCOL

N ( 5 ) =NUMNPL+ ( J-NTY 1 L+ 1

)

N( 5) =MIN0(N(5) , NUMNP

)

GO TO 400 C

C IN UPPER WALL

C 250 N( 1 )=NUMNPL+ ( J-NTY1L+1

)

N(2)=N( 1 )-1

N(3)=N( 1 )+1

N(4)=NUMNPL-( NCOL-J

)

N ( 5 ) = N ( 1 ) + ( NTYWAL-NTY 1 L* 1

)

GO TO 400 C

C UPPER REGION

C

300 N(1)=NUMNPL+(I-NTZFRP-1)»(NTYVAL-NTY1U1MJ-NTY1U1)

N( 2) =N( 1 )-1

N(3)=N( 1 )+1

N(4)=N(1)-(NTYWAL-NTY1U1)

N(5)=N( 1 )( NTYWAL-NTY 1L+ 1

)

N(3)=MINO(N(3),NUMNP)

N(5) =MIN0(N(5) ,NUMNP) C 400 IF (IDBUG(3).EQ. 0) RETURN WRITE (6,3000) I, J,N( 1 ) ,N(2) ,N(3) ,N(4) ,N(5)

C 3000 FORMAT (15X,"IN GETNOD - I,J,N(5) «",8I5)

C END A-51
c SUBROUTINE PRINT2 (TIME,T,N,NUMNP)

CCfllHIIIHIimilllllllHIIIIIIIIHHIIIIIIHIIHItlimillllHHI

C» SUBROUTINE TO PRINT TEMPERATURE STATES

CMIItMllllfllllllllllllllllllllllllllllllllllllllllllllllllKIIII

C DIMENSION T(N) , L(5) ,U(5) C

N5=NUMNP/5

NLEFT=NUMNP-5*N5

ND05=N5

IF (NLEFT.NE.O) ND05=N5+NLEFT-4

NSTART=ND05+1

NST0P=NSTART+4-NLEFT

C

CALL HEADER 1 WRITE (6,2000) TIME

C

DO 200 1=1 ,ND05 DO 100 J=1 ,5 L(J)=(J-1 )*NST0P+I

100 U(J)=T(L(J))

200 WRITE (6, 2001 )L(1),U(1),L(2),U(2),L(3),U(3),L(4) t U(4) f L(5),U(5)

C IF (NLEFT.EQ.0) RETURN

C

DO 300 I=NSTART, NSTOP

DO 350 J=1 ,4

L(J)=(J-1)«NST0P+I

350 U(J)=T(L(J))

300 WRITE (6, 2001 )L(1) I U(1),L(2),U(2),L(3),U(3),L(4),U(4)

C2000 FORMAT (///

• / 15X,"T EMPERATURE STATE"/

• / 15X, " AT TIME = ",E12. 5////)

2001 FORMAT (5(2X, I5,3X,E15. 7)

)

C RETURN

END

CHIHIIIIIIHHIIHIIHIHIIIIIIIIIII<IHIIHIIIIIIIII<IHIHIH<II

C SUBROUTINE L0DE2 (B,T,A,N, ALFA)

CC<HHHHHH<HHHHHHHHHHHHHHHHHHHHHH<HHHH

C« SUBROUTINE TO FORN RHS DURING TRANSIENT RUN • C<HHHHHHHHHHHHHHHHHHHHHHHHHHH<HHHHH<

C

DIMENSION B(N),T(N),A(N,N)

C

DO 200 1=1 ,N SUM=0. DO 100 J=1,N

100 SUM=SUM+A (I,J)*T(J)

200 B(I )=B(I )-SUM

C

RETURN

C

END

A-52
c

C'c

I

SUBROUTINE UPDATE (T,DT,N) CCHHHIIIIIIHIIIIIIIItllHIIIHIIIIIIHIIII

C» SUBROUTINE TO UPDATE TEMPERATURES CHMHHIiniHHimmiHHHHHIHHIIIt

C COMMON /HTCOEF/ HTC(8) ,TAMB(8) ,TAMBP(8) c DIMENSION T(N),DT(N) C

DO 100 1*1 ,N 100 T(I )*T(I )+DT(I

)

c

DO 200 1=1,8

200 TAMB(I )=TAHBP(I

)

c

RETURN C

END

Cc—

c

cCM

c«CM

c

c

c

c

c

c

SUBROUTINE HPT (A,C,N,I, J, A0, A1 ,A2,A3,A4, CO)

SUBROUTINE TO WRITE A POINT TO THE A MATRIX

DIMENSION NN(5),AN(5) DIMENSION A(N,N),C(N,N)

AN( 1 )=A0 AN(2)=A1

AN(3)=A2

AN ( 4 ) =A3 AN(5)*A4

CALL GETNOD (I,J,NN) CALL AMATRX (A,N,NN, AN) C(NN(1) t NN(1))*C(NN(1),NN(1) )+C0

RETURN

END

t

SUBROUTINE WMAT1 (A,C,N,I,NC1,NC2,A0,A1 ,A2,A3,A4,C0)

DIMENSION AN(5),NN(5) DIMENSION A(N,N) ,C(N,N)

AN( 1)=A0 AN(2)*A1

AN(3)«A2

AN(4)«A3

Af(5)«A4

DO 100 J=NC1,NC2

CALL GETNOD (I.J.NN) CALL AMATRX (A,N,NN, AN) C(NN( 1 ) ,NN( 1) )cC(NN( 1) , NN( 1) )+C0

100 CONTINUE

RETURN

END

A-53
ooooooo

o

n

o

oo

o

o

oooooo

ooo

SDBRODTINE WMAT2 (A,C,N,NR1 ,NR2,J,A0,A1 ,A2,A3,A4,

• CO)

SUBROUTINE TO LOGICALLY ASSEMBLE A ROW OF THE A MATRIX

DIMENSION AN(5),NN(5)

DIMENSION A(N,N) , C(N,N)

AN( 1 )=A0

AN(2)=A1

AN(3)=A2

AN(l»)=A3

AN(5)=A4

DO 100 I=NR1,NR2

CALL GETNOD (I,J,NN)

CALL AMATRX (A,N,NN,AN)

C(NN(1),NN(1))=C(NN(1),NN(1) )+C0

100 CONTINUE

RETURN

END

SUBROUTINE WRTMAT (A,C, N,NR1 ,NR2,NC1 , NC2.A0, A1 ,A2,

1 A3.A4.C0)

SUBROUTINE TO WRITE ROWS AND COLUMNS OF GEOMETRY TO

THE A MATRIX

DIMENSION NN(5),AN(5)

DIMENSION A(N,N) ,C(N,N) C

AN( 1 )=A0

AN(2)=A1

AN ( 3 ) =A2 AN(4)=A3

AN(5)=A4

C DO 100 I=NR1,NR2

DO 100 J=NC1 , NC2 CALL GETNOD (I,J,NN)

CALL AMATRX (A,N,NN,AN)

100 CONTINUE

C RETURN

C

END

A-54
non

o

o

o

ooo

o

oo

oooooo

oooooooo

$BATCH

SUBROUTINE ASEMB2 (C,A,NN)

SUBROUTINE TO FORM ELEMENTS OF THE A MATRIX

COMMON /NUMZ/ NZPOLY, NZFRP, NZWALL COMMON /NUMY/ NY5, NY4,NY1L,NYWALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTY4,NTY1L,NTYWAL,NTY1R,NTY7,NCOL

COMMON /GRIDBZ/ NTZFRP.NTZWAL

COMMON /ZCOOR/ ZPOLY , ZFRP , ZWALL COMMON /YCOOR/ Y5.Y4.Y1L, YWALL, Y1R, Y7.Y8

COMMON /DBUG/ IDBUG(IO)

DIMENSION A(NN,NN),C(NN,NN)

1. INTERIOR NODE POINTS

- REGION 5.

DY=Y5/FLOAT(NY5-1)

DZ=ZPOLY/FLOAT( NZPOLY- 1

)

NCI =2 NC2=NY5-1

NR 1=2 NR2=NZP0LY-1

CALL AINT (DY,DZ,1,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A, C, NN, NR 1 , NR2, NC 1 ,NC2, AO, A 1, A2 1 A3, A4 , CO)

- REGION 4L

DY=(Y4-Y5)/FLOAT(NY4)

DZ=ZPOLY/(NZPOLY-1

)

NC1=NY5+1

NC2=NTY4-1

CALL AINT (DY,DZ,3,A0,A1 ,A2,A3,A4,C0)

CALL WRTMAT (A,C,NN, NR1 ,NR2,NC1 ,NC2,A0,A1 , A2, A3, A4.C0)

- REGION 40

DZ= ( ZFRP-ZPOLY ) /NZFRP

NRUNZPOLY+1

NR2=NTZFRP-1

CALL AINT (DY,DZ,3,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A,C,NN,NR1 ,NR2,NC1 ,NC2,A0,A1 ,A2,A3,A4,C0)

- REGION 1L

DY= ( Y1L-Y4) /FLOAT ( NY 1L) DZ= ZPOLY/FLOAT ( NZPOLY- 1

)

NR2=NZP0LY-1

NC1=NTY4+1

NC2=NTY1L-1

NR 1=2 CALL AINT (DY.DZ, 1 ,A0, A1 ,A2,A3, A4.C0) CALL WRTMAT (A,C,NN, NR1 ,NR2,NC1 ,NC2, A0,A1 ,A2, A3, A4.C0)

- REGION 1C

A-55
ooo

oo

o

ooo

ooo

n

no

ooo

DY= ( YWALL-Y 1 L ) /FLOAT ( NYWALL

)

DZ= ZPOLY/FLOAT ( HZPOLY- 1

)

NC1=NTY1L+1

NC2=NTYWAL-1

NR 1=2 NR2=NZP0LY-1

CALL AINT (DY,DZ,1,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A, C,NN, NR1 ,NR2,NC1 ,NC2, AO, A1 , A2 f A3, AH, CO)

- REGION 1R

DY= ( Y 1 R-YWALL ) /FLOAT ( NY 1 R

)

NC 1 =NTYWAL+1

NC2=NTY1R-1

CALL AINT (DY,DZ, 1 ,A0,A1 , A2,A3» A4,C0) CALL WRTMAT (A,C,NN,NR1,NR2,NC1,NC2,A0,A1,A2,A3,A4,C0)

- REGION 7L

DY = ( Y7-Y 1 R ) /FLOAT (NY7) NC 1 =NTY 1 R+

1

NC2=NTY7-1

NR2=NZP0LY-1

CALL AINT (DY,DZ,4,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A,C,NN,NR1 ,NR2,NC1 ,NC2,A0,A1 ,A2,A3,A4,C0)

- REGION 70

DZ= ( ZFRP-ZPOLY ) /NZFRP

NR1=NZPOLY+1

NR2=NTZFRP-1

CALL AINT (DY,DZ,4,A0,A1 f A2 f A3,A4,C0) CALL WRTMAT (A,C,NN,NR1 ,NR2,NC1 ,NC2,A0,A1 ,A2, A3.A4.C0)

- REGION 8

DY= ( Y8-Y7 ) /FLOAT ( NY8

)

DZ=ZP0LY/FL0AT(NZP0LY-1

)

NC1=NTY7+1

NC2=NCOL-1

NR 1=2 NR2=NZPOLY-1

CALL AINT (DY,DZ,1,A0,A1,A2 f A3,A4,C0)

CALL WRTMAT (A, C,NN, NR1 ,NR2,NC1 ,NC2,A0,A1 ,A2,A3, A4.C0)

- REGION 6

DY=T5/FL0AT(NY5-1)

DZ= ( ZFRP-ZPOLY ) /FLOAT ( NZFRP

)

NCI =2 NC2=NY5-1

NR1=NZP0LY+1

NR2=NTZFRP-1

CALL AINT (DY,DZ,2,A0 t A1 f A2,A3,A4 t C0) CALL WRTMAT (A,C,NN,NR1 ,NR2,NC1 ,NC2,A0,A1 , A2,A3,A4,C0)

- REGION 2L

DY= (Y1L-Y4) /FLOAT (NY 1L)

NC1=NTY4+1

NC2=NTY 1L-1 CALL AINT (DY,DZ,2,A0,A1,A2 f A3,A4,C0)

CALL WRTMAT (A,C,NN,NR1 ,NR2,NC1 ,NC2, AO, A1 , A2, A3, A4 ,C0)

A- 5 6
ooo

ooo

ooo

oo

o

o

o

ooo

ooo

ooo

no - REGION 2C

DY* ( YWALL-I 1 L ) /FLOAT ( NTHALL

)

NC1=NTY1L+1

NC2=NTYHAL-1

CALL AINT (DY,DZ,2,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A, C, NN, NR1 , NR2.NC1 ,NC2, AO, A1 , A2, A3, A4.C0)

- REGION 2R

DY= ( Y 1 R-YVALL ) /FLOAT ( NY 1 R

)

NC1=NTYWAL+1

NC2=NTY1R-1

CALL AINT (DY,DZ,2,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A,C, NN, NR1 , NR2.NC1 ,NC2, AO, A1 , A2, A3, A4.C0)

- REGION 9

DY= (Y8-Y7) /FLOAT (NY8)

NC1=NTY7+1

NC2=NC0L-1 CALL AINT (DY,DZ,2,A0,A1,A2,A3,A4,C0)

CALL WRTMAT (A,C,NN,NR1,NR2,NC1,NC2,A0,A1,A2,A3,A4,C0)

- REGION 3

DY =( YWALL-Y 1 L ) /FLOAT ( NYWALL

)

DZ=( ZWALL-ZFR P ) /FLOAT ( NZWALL

)

NC1=NTY1U1

NC2=NTYHAL-1

NR1=NTZFRP+1

NR2=NTZWAL-1

CALL AINT (DY,DZ,5,AO,A1,A2,A3,A4,CO)

CALL WRTMAT (A,C,NN,NR1 ,NR2,NC1 ,NC2,A0,A1 ,A2,A3,A4,C0)

2. CONVECTION BOUNDARIES

- REGION 5B

DY=Y5/FLOAT(NY5-1)

DZ=ZP0LY/FL0AT(NZP0LY-1

)

1=1

NCI =2 NC2=NY5-1 CALL ACONV (3, 1 ,DY,DZ, 1 ,A0,A1 ,A2,A3,A4,C0) CALL WMAT1 ( A, C , NN, I, NC 1

, NC2.A0, A1 , A2, A3, A4 ,C0)

LINE 4B

DY= (Y4-Y5) /FLOAT (NY4) NC 1 =NY5+

1

NC2=NTY4-1

CALL ACONV (3, 1 ,DY,DZ,3,A0,A1 ,A2,A3, A4.C0) CALL WMAT1 (A,C, NN,I,NC1 ,NC2,A0, A1 , A2,A3,A4, CO)

- LINE 1LB

DY= ( Y1L-Y4) /FLOAT (NY 1L)

NC1=NTY4+1

NC2=NTY1L-1

CALL ACONV (3, 1 ,DY,DZ, 1 ,A0,A1 , A2, A3, A4 ,C0) CALL WMAT1 (A, C, NN, I, NCI ,NC2, A0,A1 , A2, A3, A4 , CO)

- LINE 1CB

DY= ( YWALL-Y 1L ) /FLOAT ( NYWALL

)

NC 1 =NTY 1L+1 NC2 = NTYWAL— 1 CALL ACONV (3, 1 ,DY,DZ, 1 ,A0,A1 , A2,A3, A4.C0) CALL WMAT1 (A,C,NN,I,NC1,NC2,A0,A1,A2,A3,A4,C0)

A-57
non

non

ooo

ooo

ooo

ooo

non - LINE 1RB

DY= ( Y 1 R-YVALL ) /FLOAT ( NY 1 R

)

NC1=NTYWAL+1

NC2=NTY 1R-1 CALL ACONV (3, 1 ,DY,DZ, 1 ,A0,A1 , A2,A3, A4.C0) CALL UMAT1 (A,C,NN,I,NC1,NC2,A0,A1,A2 f A3,A4,C0)

- LINE 7B

DY=(Y7-Y1R)/FL0AT(NY7)

DZ=ZP0LY/(NZP0LY-1

)

NCUNTY1R+1

NC2=NTY7-1

CALL ACONV (3, 1 ,DY,DZ,4, A0,A1 ,A2,A3,A4,C0) CALL WMAT1 (A,C,NN,I,NC1,NC2,A0,A1,A2,A3,A4,C0)

- LINE SB

DY= ( Y8-Y7 ) /FLOAT (NY8)

DZ=ZPOLY/FLOAT(NZPOLY-1

)

NCUNTY7+1

NC2=NC0L-1

CALL ACONV (3, 1 ,DY,DZ, 1 ,A0,A1 ,A2,A3,A4,C0)

CALL WMAT1 (A,C,NN,I,NC1,NC2,AO,A1,A2,A3,A4,CO)

- LINE 5L

DY=Y5/FLOAT(NY5-1)

DZ=ZPOLY/FLOAT(NZPOLY-1

)

Jsl NR 1=2 NR2=NZP0LY-1

CALL ACONV ( 1,8,DZ,DY, 1 ,A0,A1 ,A2,A3,A4,C0) CALL WMAT2 (A, C,NN f NR1 ,NR2, J, A0,A1 ,A2, A3, A4 ,C0)

- LINE 6L

DZ= ( ZFRP-ZPOLY ) /FLOAT (NZFRP

)

NR1=NZP0LY+1

NR2=NTZFRP-1

CALL ACONV ( 1 ,8,DZ,DY,2,A0,A1 ,A2,A3, A4.C0)

CALL UMAT2 (A,C, NN,NR1 , NR2,J,A0,A1 ,A2, A3.A4.C0)

- LINE 6T

DY=Y5/FLOAT(NY5-1)

DZ= ( ZFRP-ZPOLY ) /FLOAT (NZFRP

)

IsNTZFRP

NCI =2 NC2=NY5-1

CALL ACONV (4,7,DY,DZ,2,A0,A1 ,A2,A3,A4, CO) CALL WHAT1 (A, C, NN, I, NCI ,NC2 , AO, A1 , A2, A3, A4 , CO)

- LINE 4T

DY= (Y4-Y5) /FLOAT (NY4)

NC1=NY5+1

NC2=NTY4-1

CALL ACONV (4,7,DY,DZ,3,A0,A1 , A2,A3,A4,C0)

CALL WMAT1 (A,C,NN, I, NCI ,NC2,A0,A1 ,A2,A3,A4,C0)

A-58
V

C - LINE 2LT

C

DY= ( Y1L-Y4) /FLOAT ( NY 1L) NC1=NTY4«-1 NC2:NTY1L-1

CALL ACONV (4,7,DY,DZ,2,A0,A1 ,A2,A3,A4,C0) CALL WMAT1 (A, C, NN,I, NCI ,NC2 ,AO, A1 , A2 , A3, A4.C0) C

C - LINE 3L

C

DY : ( YVALL-Y 1 L ) /FLOAT ( NYWALL

)

DZ = ( ZVALL-ZFRP ) /FLOAT ( NZWALL

)

NR1=NTZFRP+1

NR2=NTZWAL-1

J=NTY1L

CALL ACONV ( 1 , 6 ,DZ,DY,5, AO, A1 ,A2,A3, A4 , CO) CALL WMAT2 (A,C,NN,NR1,NR2,J. A0,A1,A2,A3,A4,C0)

C

C - LINE 3T C IsNTZWAL NC1=NTY1U1

NC2=NTY¥AL-1

CALL ACONV (4,5,DT,DZ,5,A0,A1,A2,A3,A4,C0)

CALL WMAT1 ( A, C, NN, I, NC 1, NC2 , AO, A1 , A2 , A3 , A4 , CO)

C

C - LINE 3R C

NR1=NTZFRP+1

NR2=NTZWAL-1 J=NTYNAL

CALL ACONV (2,4,DZ,DY,5,A0,A1,A2,A3,A4,C0)

CALL WMAT2 (A, C,NN, NR1 ,NR2 , J. AO,A1 , A2,A3,A4,C0)

C

C - LINE 2RT C

DY: ( Y 1 R-YWALL ) /FLOAT ( NY 1 R

)

DZ: ( ZFRP-ZPOLY ) /FLOAT ( NZFRP

)

IrNTZFRP

NC1:NTYVALf1

NC2:NTY1R-1

CALL ACONV (4,3,DY,DZ,2,AO,A1 ,A2,A3,A4,CC) CALL WMAT1 (A,C,NN,I,NC1,NC2,A0,A1,A2,A3,A4,C0)

C

C - LINE 7T C

DY: (Y7-Y1R) /FLOAT (NY7)

NC1:NTY1R+1

NC2:NTY7-1

CALL ACONV (4,3,DY,DZ,4,A0,A1, A2,A3,A4,C0)

CALL WMAT1 (A,C,NN,I,NC1 ,NC2,A0,A1 ,A2,A3,A4,C0)

C

C - LINE 9T

C

DY: ( Y8-Y7 ) /FLOAT (NY8)

NC1:NTY7+1

NC2=NCOL-1 CALL ACONV (4,3,DY,DZ,2,A0,A1,A2,A3,A4,C0)

CALL WMAT1 (A,C,NN,I,NC1,NC2,A0,A1 ,A2,A3,A4,CO)

. C

C - LINE 9R C

J:NCOL

NR 1:NZPOLYPI NR2:NTZFRP-1

CALL ACONV (2,2,DZ,DY,2,A0,A1 ,A2,A3,A4,C0) CALL WMAT2 (A,C, NN, NR1 ,NR2, J. A0,A1 ,A2,A3, A4.C0)

C

A-59
ooo

ooo

ooo

n

o

o

ooo

oooooo

or - LINt OK

DZ=ZP0LY/FL0AT(NZP0LY-1

)

NR 1=2 NR2=NZP0LY-1

CALL ACONV (2,2,DZ,DY, 1 ,A0,A1 ,A2,A3,A4,C0)

CALL VHAT2 (A, C,NN, NR1 ,NR2, J, AO, A 1, A2, A3, A4, CO)

INTERFACE CONDITIONS

LINE 5-6

DX=(Y5)/FLOAT(NY5-1)

DIA= ( ZFRP-ZPOLY ) /FLOAT ( NZFRP

)

DXB= ZPOLY/FLOAT(NZPOLY- 1

)

IrNZPOLY

NCI =2 NC2=NY5-1 HA=2 MB=1 CALL AINF ( 1 ,MA,MB,DX,DXA,DXB, A0,A1 ,A2,A3, A4,C0) CALL WMAT1 (A, C,NN f I, NCI , NC2, AO, A1 ,A2,A3 f A4 ,C0)

INTERFACE 4L4U

DX=(Y4-Y5)/NY4

NC 1=NY5+1

NC2=NTY4-1

MA=3 MB=3 CALL AINF (1,MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WMAT1 (A,C,NN,I,NC1 ,NC2,A0,A1 ,A2,A3, A4.C0)

- INTERFACE 12L

NC1=NTYM+1

NC2=NTY1L-1

MA=2

MB= 1 DX=( Y1L-Y4) /FLOAT ( NY 1L) CALL AINF ( 1 ,MA,MB,DX,DXA,DXB, A0,A1 ,A2,A3, A4.C0) CALL VMAT1 (A, C,NN,I,NC1 ,NC2,A0,A1 , A2,A3,A4,C0)

- INTERFACE 12C

NC1=NTY1L+1

NC2=NTYWAL-1

DX = ( YWALL-Y 1 L ) /FLOAT ( NYWALL

)

CALL AINF (1,MA f MB,DX,DXA,DXB,A0,A1,A2 t A3,A4,C0)

CALL WMAT1 (A,C,NN,I, NCI ,NC2, AO, A1 , A2, A3,A4 ,C0)

- INTERFACE 12R

NC1=NTYWAL+1

NC2=NTY1R-1

DX= ( Y 1 R-YVALL ) /FLOAT ( NY 1 R

)

CALL AINF (1 I MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL VMAT1 (A,C, NN, I, NCI ,NC2, AO, A1 , A2,A3, A4,C0)

- INTERFACE 7L7U

DX=(Y7-Y1R)/NY7

NC1=NTY1R+1

NC2=NTY7-1

MA=4

MB= 4 CALL AINF (1,MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WMAT1 (A,C,NN,I,NC1 ,NC2,A0,A1 ,A2,A3,A4.C0)

A-60
C - INTERFACE 89 C

NC1=NTY7«-1 NC2=NC0L-1

DX=(Y8-Y7)/FLOAT(NY8)

MA=2

MB=1 CALL AINF ( 1 ,MA,MB,DX,DXA,DXB, A0,A1 t A2 t A3fA4,C0)

CALL VMAT1 ( A , C , NN, I , NC 1 ,NC2 , AO, A 1 , A2 , A3 , A4 , CO) C

C - 23 C

I=NTZFRP NCUNTY1U1

NC2=NTYVAL-1

MA=5

MB=2

DX = ( TVALL-T 1 L ) /FLOAT ( NYVALL

)

DIA = ( ZWALL-ZFRP ) /FLOAT ( NZWALL

)

DXB= ( ZFRP-ZPOLY ) /FLOAT ( NZFRP

)

CALL AINF ( 1 ,MA,MB,DX,DXA,DXB, AO, A1 , A2,A3,A4,CO)

CALL VMAT1 ( A, C, NN, I, NC 1, NC2 , AO, A1 f A2 , A3 , A4 ,C0) C

C - INTERFACE 54

C

J=NY5

NR 1=2 NR2=NZP0LY-1

DX=ZPOLY/FLOAT(NZPOLY-1

)

DXA=Y5/FLOAT(NY5-1)

DXB= ( Y4-Y5 ) /FLOAT ( NY4

)

MA=1

MB=3

CALL AINF (3,MA,HB,DX f DXA f DXB,AO f A1,A2 f A3,A4,CO)

CALL WMAT2 (A f C,NN, NR1 ,NR2,J, AO, A1 , A2, A3, A4.C0) C

C - INTERFACE 64 C

NR1=NZPOLY+1

NR2=NTZFRP-1

DX =( ZFRP-ZPOLY ) /FLOAT ( NZFR P

)

HA=2 CALL AINF (3,MA,MB t DX,DXA,DXB, A0,A1 , A2, A3, A4,C0) CALL VMAT2 (A,C,NN,NR1 ,NR2, J.AO.A1 ,A2,A3,A4,C0) C

C - INTERFACE 41L C

NR1=2

NR2=NZPOLY-1

J=NTY4

DX=ZPOLY/(NZ POLY-1

)

DXA=(Y4-Y5)/NY4

DXB= (Y1L-Y4) /FLOAT ( NY 1L)

MA=3

MB=1 CALL AINF (3tMA t MB,DX f DXA,DXB t A0,A1 , A2, A3« A4,C0) CALL VMAT2 (A,C, NN, NR1 ,NR2, J. A0,A1 , A2.A3, A4.C0) C

C - INTERFACE 42L

C

NR 1 =NZPOLY+1

NR2=NTZFRP-1

DX=( ZFRP-ZPOLY ) /FLOAT ( NZFRP

)

MA=3

MB=2 CALL AINF (S.MA^.DX^DXA.DXB.AO.AI.AE.AS.AA.CO)

CALL VMAT2 (A,C,NN,NR1 ,NR2, J.A0.A1 ,A2,A3, A4.C0) C

A-61
ooo

ooo

ooo

ooo

ooo

C - INTERFACE 1L1C

c NR1=2

NR2*NZP0LY-1 JsNTYIL DXs(ZPOLY)/FLOAT(NZPOLY-1)

DXAs(Y1L-Y4)/NY1L

DiB=(YWALL-Y 1L)/NYVALL

MA= 1 MBs 1 CALL AINF (3,MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WMAT2 (A, C,NN, NR1 , NR2 , J, AO, A1 , A2 , A3, A4 ,C0)

- INTERFACE 2L2C

NR1 sNZPOLY+1

NR2sNTZFRP-1

DXs(ZFRP-ZPOLY)/NZFRP

MA=2

MB=2

CALL AINF (3,MA,MB,DX,DXA,DXB,AO,A1,A2,A3,A4,CO)

CALL WMAT2 (A, C,NN, NR1 , NR2 , J. AO, A1 ,A2, A3, A4.C0)

- INTERFACE 1C1R

NR1s2

NR2sNZP0LY-1 JsNTYWAL DXs(ZPOLY)/(NZPOLY-1)

DXA s ( YWALL-Y 1 L ) /NYWALL DXBs ( Y 1 R-YVALL ) /NY 1

R

MAsI

MBs 1 CALL AINF (3,MA,MB,DX,DXA,DXB, A0,A1 , A2,A3,A4,C0)

CALL WMAT2 (A, C,NN, NR1 ,NR2 t J, AO, A1 , A2, A3, A4 ,C0)

- INTERFACE 2C2R

NR 1 sNZPOLY+1

NR2=NTZFRP-1 JsNTYWAL

DXs( ZFRP-ZPOLY ) /NZFRP

MAs2

MBs2 CALL AINF (3,MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WMAT2 (A,C,NN,NR1 ,NR2,J,A0,A1 ,A2,A3,A4,C0)

- INTERFACE 1R7

NR1s2

NR2sNZPOLY-1 JsNTYIR DXsZPOLY/(NZ POLY-1)

DXA=(Y1R-YWALL)/NY1R

DXB=(Y7-Y1R)/NY7

MAsI MBs4 CALL AINF (3,MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WMAT2 (A,C,NN,NR1 ,NR2,J.A0,A1,A2,A3,A4,C0)

- INTERFACE 2R7

NR 1 sNZPOLY+1

NR2sNTZFRP-1

DX=( ZFRP-ZPOLY ) /NZFRP

MAs2

MBs4 CALL AINF (3,MA,MB,DX,DXA,DXB,AO,A1,A2,A3,A4,CO)

CALL WMAT2 (A,C,NN,NR1 ,NR2, J.A0.A1 ,A2,A3,A4,C0)

A-62
ooo

ooo

ooooooo

ooo

oo - INTERFACE 78

NR 1*2 NR2*NZPOLY-1

J=NTY7

DX=ZPOLY/(NZPOLY-1)

DIA=(Y7-Y1R)/NY7

DXB=(Y8-Y7)/NY8

MA*4

MB= 1 CALL AINF (3,MA,HB,DX,DXA,DIB,A0,A1,A2,A3,A4,C0)

CALL WMAT2 (A,C, NN, NR1 ,NR2, J, A0,A1 , A2,A3, A4, CO)

- INTERFACE 79

NR1*NZP0LYf1

NR2=NTZFRP-1

DX=( ZFRP-ZPOLY )/NZFRP

MB=2

CALL AINF (3,MA,MB,DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL UMAT2 (A,C, NN, NR1 ,NR2 f J, AO, A1 , A2, A3, A4.C0)

4. CORNERS

- INTERIOR CORNERS, ALL CONDUCTION

- 654L4U

I=NZPOLT

J*NY5

DX1=Y5/(NY5-1)

DX2=(Y4-J5)/NY4

DI3=ZP0LY/(NZP0LY-1)

DX4= ( ZFRP-ZPOLY )/NZFRP

MA=3

MB=2

MC=1

MD=3

CALL AIC1 (MA,MB,MC,MD,DX1,DX2,DX3,DX4,A0,A1,A2,A3,A4,C0)

CALL VPT (A,C,NN,I,J,A0,A1,A2,A3,A4,C0)

- 4L401L2L

IcNZPOLY

J=NTY4

DX1=(Y4-Y5)/NY4

DX2=(Y1L-Y4)/NY1L

MA=2

MB=3

MC=3

MD=1

CALL AIC1 (MA ,MB,MC,MD,DX1 ,DX2,DX3,DX4, AO, A1 , A2, A3, A4,C0) CALL VPT (A,C,NN,I,J,A0,A1,A2,A3,A4,C0)

- 2C2L1L1C

J=NTY1L

DX1=(Y1L-Y4)/NY1L

DX 2= ( YVALL-Y 1 L ) /NYVALL

HA =2 MB=2 MCkI

KD=1

CALL AIC1 (HA,HB,HC,MD,DX1,DX2,DI3,DX4,A0,A1,A2,A3,A4,C0)

CALL VPT U.C.NN.I.J.AO.AI'AP^.AA.CO)

A- 63
non

no

o

ooo

oo 2R2C1C1R

JsHTYVAL

DX 1 = ( YVALL-Y 1L ) /HYWALL

DX2=(Y1R-YVALL)/HY1R

CALL AIC1 (MA,MB,MC,MD f DX1,DI2 f DX3,DX4,A0 I A1,A2,A3 ) AA,C0) CALL VPT (A,C f HH,I, J, AO, A1 , A2, A3, A4.C0)

- 72B1B

J=HTY1B

DX1=(Y1R-YVALL)/NY1R

DX2=(Y7-Y1R)/HY7

HA =4 HB=2

MC=1

MD=4

CALL AIC 1 (MA,HB I HC,KD,DX1,DX2,DX3,DX4 I A0,A1 I A2 I A3,A4 I C0) CALL VPT (A,C,NH,I, J,A0 ( A1,A2,A3,A4,C0)

- 978

J=HTT7

DX1=(Y7-Y1R)/HY7

DX2=(Y8-Y7)/HY8

MA=2

KB=4

MC=4

HD=1 CALL AIC1 (HA , HB , HC , HD , DX 1

, DX2 , DX3 » DX4 , AO , A 1

, A2 , A3 , A4 , CO

)

CALL VPT (A f C,HN,I t J,A0,A1,A2,A3,A4,C0)

- REENTRANT CORKERS

I=NTZFRP JsNTXIL

N1 = 1 H2=4

KR1=7

KH2=6

DX1=(Y1L-Y4)/NY1L

DX2= ( YVALL-Y 1L ) /NYWALL

DX3= ( ZFRP-ZPOLY ) /KZFRP

DX 4= ( ZVALL-ZFR P ) /HZWALL

HA =5 HB=6

HC=2 HD=2 CALL AIC2 (N1,N2,HH1,NB2,HA,HB,HC,HD,

• DX1 f DX2,DX3,DX4,A0,A1,A2,A3,A4,C0)

CALL VPT (A,C,HH,I f J,A0,A1,A2,A3,A4,C0)

C JsHTYVAL

HI =2 2=4

HH 1 =3 KH2=4

DX1=( YVALL-Y 1L)/HYVALL

DX2=(Y1B-YVALL)/KY1R

HA =6 MB=5

MC=2

MD=2

CALL AIC2 (H1 f H2,KH1,HH2 t HA f HB,HC f HD,

• DX1 f DX2,DX3,DX4,A0,A1,A2,A3,A4 f C0) CALL VPT (A,C,HH l I l J > A0 l A1 t A2 l A3,A4 t C0)

A-64
c

c CONVECTION ON ONE SIDE

J=1

IbKZPOLI

Ns 1 NHs8

DXsI5/(NI5-1

)

DXA =ZPOLY/ ( NZ POLI- 1

)

DXBs ( ZFRP-ZPOLI ) /NZPKP HAsi MBs2 CALL AIC3 (W.NH.MA.HB,

• DX,DIA,DIB, A0,A1 , A2,A3,A4,C0)

CALL VPT (A,C,NX,I,J,A0,A1,A2,A3,A4,C0)

Id

JsNT5

N=3 KHsI DXsZPOLI/(NZPOLI-1)

DXB=T5/(NT5-1)

DIAs(I4-I5)/m

KBsl

MA=3

CALL AIC3 (M.NH.HA.HB,

• DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL VPT (A,C,NN, I, J.A0.A1 ,A2,A3,A4,C0)

IsMTZPEP

Ns4 NH=7 DXs ( ZFHP-ZPOLT ) /NZFBP DXAsDXB

DXB= (Y4-T5 )/NT4 HA =2 KB=3 CALL AIC3 (N,NH,HA,MB f • DX, DXA, DIB, A0,A1,A2,A3,A4,C0)

CALL VPT (A f C,NN,I, J,A0,A1 ,A2,A3,A4,C0)

1=1 JsNTT4

M=3

VHsI DXsZPOLT/ ( NZPOLT-1

)

DXB=(T4-T5)/NY4

DXA=(T 1L-T4) /NT 1L MBs 3 MA = 1 CALL AIC3 (N,NH,MA,MB,

• DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL VPT (A,C, NN,I, J, A0,A1 ,A2,A3, A4,C0)

IsNTZFHP

1=4

IH=7 DI= ( ZFRP-ZPOLI ) /NZPRP DXAsDXB DXBr(T1L-T4)/NT1L

HA =3 MBs2

CALL AIC3 (N.NB.HA.MB,

• DX , DXA , DXB ,A0,A1,A2,A3iA4,C0)

CALL VPT (A,C,NN,I,J,A0,A1,A2,A3.A4,C0)

A- 65
1=1 JsHTYIL

N=3

NH= 1 DX=ZPOLY/(NZPOLY-1)

DXB=(Y1L-Y4)/NY1L

DXA = ( YWALL-Y 1 L ) /HIVALL

MA = 1 MB= 1 CALL AIC3 (N,NH,MA,HB,

• DX, DXA , DYE, A0,A1,A2,A3'A4,C0)

CALL WPT (A,C f NN,I,J,A0,A1,A2,A3,A4 f C0)

J=NTYWAL

DXB= ( YWALL-Y 1 L ) /NYWALL DIA=(Y1R-YWALL)/NY1R

CALL AIC3 (N,NH,MA,MB,

• DX,DIA,DXB,AO,A1 f A2,A3,A4,CO)

CALL WPT (A,C,NN,I,J,A0,A1,A2,A3,A4,C0)

J=NTY1R

DXB=(Y1R-YWALL)/NY1R

DXA=(Y7-Y1R)/NY7

HB=1

MA=4

CALL AIC3 (N, NH.MA,MB,

• DX,DXA f DXB, A0,A1,A2,A3,A4 f C0) CALL WPT (A,C,KN,I,J,A0,A1,A2,A3,A4 f C0)

I=NTZFRP

11=4 NH=3

DX=(ZFRP-ZPOLY)/NZFRP

DXA=DXB DXB=(Y7-Y1R)/NY7

HA =2 MB=4 CALL AIC3 (N,NH,HA,MB,

• DX,DXA,DXB,A0 f A1,A2,A3,A4,C0)

CALL WPT (A f C,NN,I,J f A0,A1 f A2,A3,A4,C0)

1=1 J=NTY7

N=3

NH=1 DX=ZP0LY/(NZP0LY-1

)

DXB=(Y7-Y1R)/NY7

DXA=(Y8-Y7)/NY8

MB=4 MA=1 CALL AIC3 (N,NH,MA,MB,

• DX,DXA,DXB, A0,A1 ,A2,A3,A4,C0) CALL WPT (A,C,NN,I,J,A0,A1 f A2,A3,A4,C0)

I=NTZFRP

H=4

HH=3

DX=( ZFRP-ZPOLY ) /NZFRP DXA=DXB

DXB=(Y8-Y7)/NY8

HA=4 KB=2

CALL AIC3 (N ( NH,MA,MB,

• DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C,NN,I,J,A0 f A1,A2,A3,A4,C0)

I=NZPOLY

JsNCOL

N=2

NH=2

DX=(Y8-Y7)/NY8

DXB=ZPOLY/(NZPOLY-1

)

DXA= ( ZFRP-ZPOLY ) /NZFRP

MB= 1 HA=2 CALL AIC3 (N.NH.MA.HB,

• DX,DXA,DXB,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C,NN,I,J,A0,A1,A2,A3,A4,C0)

A-66 C
C - EXTERIOR CORNERS c

I»1 J«1

N1«1 N2=3

NH1=8

NH2>1

DX1=Y5/(NY5-1)

DX 2= ZPOLY/ ( NZ POLY - 1

)

M=1 CALL AIC4 (N1,N2,NH1,NH2,M,

• DX1,DX2,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C,NN,I,J,A0,A1,A2,A3,A it,C0)

C

1=1 J=NCOL

N1=2

N2=3

NH1=2

NH2=1

DX1=(Y8-Y7)/NY8

CALL AIC4 (N1 ,N2,NH1 ,NH2,H,

• DX1,DX2,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C,NN,I,J,A0,A1,A2 f A3,A4,C0)

C I=NTZFRP

N1=2

N2=4

NH1=2

NH2=3

DX2= ( ZFRP-ZPOLY ) /NZFRP

M=2

CALL AIC4 ( N 1 ,N2,NK1 ,NH2,M,

• DX1,DX2,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C,NN,I,J f A0,A1,A2,A3,A4,C0)

C I=NTZWAL JrNTYWAL

N1=2

N2=4

NH1=4

NH2=5 DX 1 = ( YWALL-Y 1 L ) /NYWALL DX2= ( ZWALL-ZFRP ) /NZWALL

M=5 CALL AIC4 (N1 ,N2,NH1 ,NH2,M,

• DX1,DX2,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C,NN,I,J,A0 t A1,A2,A3,A4,C0)

C

J=NTY1L

N1 = 1 N2=4

NH1=6

NH2=5

CALL AIC4 (N1 ,N2,NH1 , NH2,M,

• DX1 ,DX2, AO, A1 ,A2,A3i A4,C0) CALL WPT (A,C,NN,I t J,A0,A1 f A2,A3,A4,C0)

C I=NTZFRP

J=1

N1 = 1 N2=4

NH1=8

NH2=7

DX1=Y5/(NY5-1)

DX2= ( ZFRP-ZPOLY ) /NZFRP

M=2

CALL AIC4 (N1,N2,NH1,NH2,M,

• DX1,DX2,A0,A1,A2,A3,A4,C0)

CALL WPT (A,C t NN, I, J, A0 ( A1 ,A2 t A3»A4,C0)

C

RETURN

C

END A-67
ooooo

ono

c SUBROUTINE BCONV (N,M,DY,DZ,MM,BO)

CCII*llllt«IOItt«ll«H«IIIIHIIIIHHIIIIII«lltl*<

C • SUBROUTINE TO FORM B VECTOR - LINES CHHIimiHIHmilllHHIIIIIHHIIIimtllftK

c

c

c

c

COMMON /HTCOEF/ HTC (8) ,TAMB(8) ,TAMBP(8) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

BO=HTC ( M

) «DY • ( ALFA»TAMBP ( M ) + ( 1

. -ALFA ) »TAMB ( M )

)

RETURN

END

SUBROUTINE B2 (N1 ,N2,NH1 ,NH2,MA,MB,MC,MD,DX1 ,DX2,DX3,DX4,B0)

• SUBROUTINE TO FORM B VECTOR - REENTRANT CORNERS

COMMON /HTCOEF/ HTC (8) ,TAMB( 8) ,TAMBP(8) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

C DDX 1 =DX

1

DDX2=DX4

IF (N1.EQ.2) DDX 1 =DX2

A0=0. 5*HTC (NH1 )*DDX1 AO 1 =0. 5*HTC ( NH2 ) »DDX2

B0=A0*(ALFA*TAMBP(NH1 )+( 1. -ALFA) # TAMB(NH1 ) ) +

• A0 1 • ( ALFA»TAMBP ( NH2 ) t-( 1 . -ALFA )*TAMB (NH2)

)

C RETURN

C END

CCHHIIttllllHHIlHHHIHIIIIIHIIHHIflHIIIIlfHIIIHHHHHl

c SUBROUTINE B3 ( N , NH , MA , MB , DX , DXA , DXB , BO

)

C

C* SUBROUTINE TO FORM B VECTOR - CORNER NODE, CONVECTION 1 SIDE ClltHHIIIIItllHimilHHHimHIfllllHItllllHIfltllfllllHH

C

COMMON /HTCOEF/ HTC(8) ,TAMB( 8) ,TAMBP(8) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

C

BO=0 . 5 •HTC ( NH ) • ( DXA+DXB )•( ALFA *TAMBP (NH )+( 1. -ALFA) #TAMB(NH)

)

C RETURN C

END

C

C»< C

C

c»

C«’ c

C

c

c

SUBROUTINE B4 (N1 ,N2,NH1 ,NH2,M,DX1,DX2,B0)

EXTERIOR CORNERS WITH HEAT TRANSFER - B VECTOR

COMMON /HTCOEF/ HTC (8) ,TAMB(8) ,TAMBP(8) COMMON /TIM/ TIMO.TIMF, DELTA, ALFA.IRFINE

A 0=0. 5*DX2*HTC ( NH 1

)

A01=0.5*DX1 #HTC(NH2)

B0=A0*(ALFA«TAMBP(NH1 )+( 1.-ALFA )»TAMB(NH1 )

)

• AO 1 • ( ALFA»TAMBP ( NH2 ) + ( 1

. -ALFA )«TAMB ( NH2 )

)

RETURN

END A- 68
cCHIHIIIHIHIIIIIIIIHHHmHHIIIUHHHKHHHIHIIHIIHMI

c SUBROUTINE WB2 (B,NN,NR1 ,HR2,J.B0)

CCHHIIIHIIHIIHHIlKIIIHKIIIIIIHMHIIIHimiHHIHHHmi

C» SUBROUTINE TO WRITE A ROW OF BODES INTO THE B VECTOR •

C

DIMENSION B( WN) , N( 5) C

DO 100 I = NR1,NR2

CALL GETNOD (I.J.N)

100 B(N(1))sB(N( 1))+BC RETURN C

ENT

CCHHHtHllllimHHIKHHItHItllllimHIIIIHIHItHIIHIHHI

C

SUBROUTINE WB1 <B, NN, I, NC 1 ,NC2 ,B0)

CCHHHHKIHIIIIHHIlHIHimHHHIIMmHmHHHIHHHHM

C« SUBROUTINE TO WRITE A COLUMN OF NODES IN THE B VECTOR * CHIIHIIHIIIIHIHIHIIHIHitHHItlllllimiHIHIHHIlHHHIl

c

DIMENSION B(NN),N(5)

C

DO 100 J=NC1,NC2

CALL GETNOD (I.J.N)

100 B(N(1))sB(N( 1))*BC RETURN C

END

CCHHiHimHWHimmuHimiiHnumumHHHtHH

C

SUBROUTINE WB (B.IN.I. J.BO) CCMIHiiiiiiiihihiihhhhiiiihikkiiiihihiiiiiiiiiih

C» SUBROUTINE TO WRITE A POINT IN THE B VECTOR • CHHIIIHIIHIIIIIIIIIHHIHIMIHIHHHIHHIIHIIIIIHII

c

DIMENSION B(NN),N(5)

C

CALL GETNOD (I.J.N)

B(N(1))sB(N(1))«.B0 RETURN C

END

A-69
ooo

ooo

ooo

ooooooo

ooo

ooooo

cCHHHIIIHIHIHIIIHIHItHIIHHIlHIIHHHHHIHIHHHHIlH

c SUBROUTINE FORHB ( ICON, B, BTEMP, NN) CCHHtlltlMHIHIHHHIIliHIMIIIIIHtlHHIHIHHIHtMIHMHI

C« SUBROUTINE TO FORM CONVECTION PART OF RHS VECTOR • CHIHHIHIIIIIHIIIHIHIHItlHHIHIIIIIIIHHIIIimUIMHIlH

c

COMMON /NUMZ/ NZPOLY, NZFRP, NZVALL

COMMON /NUMY/ NY5,NY4,NY1L,NYWALL,NY1R,NY7,NY8

COMMON /GRIDBY/ NTY4 , NTY1L, NTYWAL, NTY 1R, NTY7 , HCOL COMMON /GRIDBZ/ NTZFRP.NTZWAL

COMMON /ZCOOR/ ZPOLY, ZFRP, ZWALL COMMON /YCOOR/ T5,Y4,Y1L,YWALL,Y1R,Y7,Y8

C

DIMENSION B(NN) ,BTEMP(NN)

ICONr-1 - SET B=0 , REFORM RHS INTO BTEMP, REWRITE BTEMP INTO B ICON: 0 - FORM RHS INTO BTEMP, WRITE BTEMP INTO B ICON: 1 - SET B=0, WRITE BTEMP INTO B

IF (ICON, EQ. 0) GO TO 100

ZERO B VECTOR

DO 10 1:1 , NN

10 B(I)=0. IF (ICON. EQ. 1) GO TO 200

FORM CONVECTION LOAD VECTOR BTEMP

2. CONVECTION BOUNDARIES

- REGION 5B

100 DY=Y5/FLOAT(NY5-1)

DZ=ZP0LY/FL0AT(NZP0LY-1

)

1:1 NC1:2

NC2:NT5-1 CALL BC0NV (3, 1 ,DY, DZ, 1 ,B0) CALL WB1 (BTEMP,NN,I,NC1,NC2,B0)

- LINE 4B

DY= (Y4-Y5) /FLOAT (NY4)

NC1:NY5+1

NC2:NTY4-1 CALL BCONV (3, 1 ,DI, DZ.3.BO) CALL WB1 (BTEMP, NN,I,NC1 ,NC2, B0)

- LINE 1LB

DY: (Y 1L-Y4 ) /FLOAT ( NY 1L

)

NC1:NTY4+1

NC2:NTY 1L-1 CALL BCONV (3, 1 ,DY,DZ, 1 ,BO) CALL WB1 ( BTEMP , NN ,

1

1 NC 1

, NC2 , B0

)

- LINE 1CB

DY: ( YVALL-Y 1 L ) /FLOAT ( NYWALL

)

NC1:NTY1L+1

NC2=NTYWAL-1

A 70
non

ooo

ooo

ooo

ooo

ooo

ooo

ooo

CALL BCONV (3, 1 ,DY,DZ, 1 ,B0) CALL VB1 ( BTEMP , NN , I , NC 1

, NC 2 , BO

)

- LINE 1 RB

DI • ( 1 1 R-YNALL ) /FLOAT ( M1 1 R

)

NC 1 =NTYWAL+1 NC2=NTY 1R-1 CALL BCONV (3, 1 ,DY,DZ, 1 ,B0) CALL WB1 ( BTEMP , NN , I , NC 1

, NC2 , BO

)

- LIME 7B

DY*(Y7-Y1R)/FLOAT(NY7)

NC1=NTY1R«-1 NC2=NTY7-1

CALL BCONV (3, 1 ,DY,DZ,4,B0) CALL WB1 ( BTEHP , NN , I , NC 1 , NC2 , BO

)

- LINE 8B

DY= (Y8-Y7) /FLOAT (NY8)

NC 1=NTY7+1 NC2=NC0L-1 CALL BCONV (3, 1 ,DT,DZ, 1 ,B0) CALL WB1 ( BTEMP , NN , I , NC 1

, NC2 , BO

)

- LINE 5L

DT=T5/FLOAT(NY5-1)

DZ=ZPOLY/FLOAT(NZPOLY-1

)

Jrl

NR 1=2 NR2=NZ POLY-1 CALL BCONV ( 1

, 8,DZ,DY, 1 ,B0) CALL VB2 (BTEMP, NN.NR1 ,NR2, J. BO)

- LINE 6L

DZ= ( ZFRP-ZPOLY ) /PLOAT ( NZFRP

)

NR1=NZP0LY+1

NR2=NTZFRP-1

CALL BCONV (1,8,DZ,DY,2,B0)

CALL WB2 (BTEMP, NN.NR1 ,NR2, J,BO)

- LINE 6T

DY=Y5/FLOAT(NY5-1)

DZ=( ZFRP-ZPOLY) /FLOAT (NZFRP) IsNTZFRP

NCI =2 NC2=NY5-

1

CALL BCONV (»,7,DY,DZ,2,B0)

CALL WB1 ( BTEMP ,NN,I,NC1,NC2,B0)

- LINE AT

DY e(Y4-Y5) /FLOAT (NY4)

NC1=NY5f1

NC2=NTY4-1 CALL BCONV ( A,7,DY,DZ,3,BO)

CALL WB1 ( BTEMP , NN , I , NC 1 , NC2 , BO

)

- LINE 2LT

DY= ( Y 1L-Y4 ) /FLOAT ( MY 1L

)

A- 71
o

o

'ooo

oo

o

o

o

o

o

o

o

ooo

ooo

non

NC 1 =NTY4+1 NC2=KTY 1L-1 CALL BCONV (4.7.DY, DZ,2,B0) CALL WB1 ( BTEMP , NN , I , NC 1 , NC2 , BO

)

- LINE 3L

DY= ( YWALL-Y 1L ) /FLOAT ( NYWALL

)

DZ= ( ZWALL-ZFRP ) /FLOAT ( NZWALL

)

NR1=NTZFRP+1

NR2=NTZWAL-1 JsNTYIL

CALL BCONV ( 1,6,DZ,DY,5,B0) CALL VB2 ( BTEMP , NN , NR 1

, NR2 , J , BO

)

- LINE 3T

IsNTZWAL NCUNTY1L+1

NC2=NTYWAL-1

CALL BCONV (4,5,DY,DZ,5,B0)

CALL WB1 (BTEMP, NN,I, NCI f NC2, BO)

- LINE 3R

NR1 =NTZFRP+1 NR2=NTZWAL-1 J=NTYWAL

CALL BCONV (2,»«,DZ,DY,5,B0) CALL WB2 (BTEMP,NN, NR1 ,NR2, J,B0)

- LINE 2RT

DY= ( Y 1 R-YWALL ) /FLOAT ( NY 1 R

)

DZ= ( ZFRP-ZPOLY ) /FLOAT (NZFRP

)

I=NTZFRP

NC1=NTYWAL+1

NC2=NTY 1R-1 CALL BCONV (H,3,DY,DZ,2,B0) CALL NB1 (BTEMP, NN,I,NC1,NC2 f B0)

- LINE 7T

DY=(Y7-Y1R)/FLOAT(NY7)

NC 1 =NTY 1R+1 NC2=NTY7-1

CALL BCONV (fl,3,DY,DZ,U,B0) CALL WB1 (BTEMP, NN,I, NCI, NC2, BO)

- LINE 9T

DY= (Y8-Y7) /FLOAT (NY8)

NC1=NTY7+1

NC2=NCOL-1 CALL BCONV (*,3,DY,DZ,2,B0)

CALL WB1 (BTEMP,NN,I,NC1 ,NC2,B0)

- LINE 9R

J=NCOL

NR1=NZPOLY+1

NR2=NTZFRP-1 CALL BCONV (2,2,DZ,DY,2,B0) CALL VB2 (BTEMP, NN.NR1 , NR2, J, BO)

- LINE 8R

A-72
ooo

o

oooo

DZ=ZPOLY /FLOAT (IZ POLY- 1

)

IF 1«2 IR2=IZ POLY-1 CALL BCOSV (2,2,DZ,DY, 1 ,B0) CALL VB2 ( BTEKP , II , IB 1 , IR2 , J , BO

)

- IEEKTRAIT COMERS

ItITZFBP J=ITY1L

11:1 12: A IH1:7

HH2:6

DX 1:(Y 1L-YA )/IY 1L DX2= ( YV ALL-I 1L ) /IYVALL DI3=( ZPBP-ZPOLY )/IZPFP

DX A : ( ZVALL-ZFB P ) /IZVALL MA:5

KB=6

HC=2

HD:2

CALL B2 ( 1 1 , 12 , !H 1 , IH2 , MA , KB , HC , KD

,

• DX1,DX2,DX3,DXA,B0)

CALL VB (BTEKP, II, I, J, BO)

JrlTYWAL

11:2

I2:A IE1:3

IH2:A

DX 1 :( TVALL-I 1L ) /IYV ALL DI2=(Y 1B-YVALL ) /IT IE HA:6 KB=5 MC=2

KD:2

CALL B2 ( 1 1, 12 , IB 1 , HE2 , HA , KB , HC , KD t • DX 1 , DX2 , DX 3 , DXA , BO

)

CALL VB (BTEKP, II , I, J, BO)

- COITECTIOI 01 OIE SIDE

J: 1 I:IZPOLT

1=1 KB=8

DX=I5/(IYS-1)

DIA=ZP0LI/(KZP0LI-1)

DIB: ( ZPBP-ZPOLY )/IZPRP

MA=1 KB=2 CALL B3 (I, IB,KA,KB,

• DX,DXA,DXB,BO)

CALL VB (BTEKP, W, I, J, BO)

1:1 J=IY5

1=3

IH=1 DX=ZPOLI/(IZPOLY-1)

DXB=Y5/(IY5-1)

DXA*(IA-Y5)/IYA

KB=1

A-73
HA =3 CALL B3 (I, IB, HA, MB,

« DX,DXA,DE,BO)

CALL ME (BTEHF

,

55 , 1 , J, BO

)

IsVTZFBP

1=4

IH=7

DX=( ZFHP-ZPOLT )/IZEBP

DIA=BE

DXB= (14-15) /IT4 HA =2 MB=3

CALL B3 (I, IE, HA, KB,

• DX,DIA,DE,BO)

CALL MB (STEM? ,15,1, J, BO)

1=1 j=m4

1=3

15=1 CI=ZPOLI/(EZPOLT-1

)

DXB=(Y4-T5)/IT4

DIA=(Y1L-T4)/H71L

MB=3

MA = 1 CALL B3 (I, SB, HA, MB,

• EX,DIA,DE,BO)

CALL HB (27EM?,S5,I, J.BO)

IsITZTBP

1=4 IB=7 DX=(ZFBP-ZPOLT)/IZFSP

DXA=DX3

DXB=(I1L-T4)/IY1L

HA=3

ME= 2 CALL B3 (I, XH,HA,HB,

• DI, CIA, EE, 30) CALL WE (BTEHP,I5,I, J.BO)

1=1 J=ITY1L

1=3 15=1 EX=ZPOLT/(IZPOLY-1

)

DE=(Y1L-Y4)/IY1L

DIA= ( Y¥ALL-T 1L ) /IT* ALL

HA = 1 H3=1 CALL B3 (I, IE, HA, MB,

• DX,EIA,EE,BO)

CALL MB (BIBS’, II, I, J.BO)

j=ittval

DX3= ( IM ALL-T 1L ) /ITMALL DIA=(T 18-TVALL ) /IT 18 CALL B3 (I, IB, HA, KB,

• DX , DXA , CE , BO

)

CALL MB (STEM?, 15, I, J.BO)

J=ITT1B

DX5=(Y 18-TVALL) /IT 18 DXA=(T7-T 1B)/IT7

HB=1

A- 7 4
HA >4 CALL B3 (N.NH.MA.MB,

• DX,DXA,DXB,BO)

CALL VB ( BTEMP , NN , I , J , BO

)

C

I*NTZFRP

M»4 MH*3

DX=( ZFRP-ZPOLY ) /NZFRP DXAcDXB

DXB=(Y7-Y1R)/NY7

HA =2 MB=4

CALL B3 ( N , NH , MA , MB

,

DX , DXA , DXB , BO

)

CALL VB ( BTEMP , NN ,

I

( J . BO

)

C

1=1 J=NTT7

N=3

NH=1 DX=ZPOLY/(NZ POLY-1)

DXB=(Y7-Y1R)/NY7

DXA=(Y8-Y7)/NY8

MB=4

HA = 1 CALL B3 ( N, NH, MA , MB,

• DX, DXA, DXB, BO) CALL VB ( BTEMP , NN , I , J , BO

)

C

I=NTZFRP

N=4 KH=3

DX=( ZFRP-ZPOLY J/NZFRP DXA=DXB

DXB=(Y8-Y7)/NY8

MA=4

MB=2

CALL B3 (N , NH, MA , MB

,

• DX, DXA, DXB, BO) CALL VB ( BTEMP , NN , I , J , BO

)

C IrNZPOLY JsNCOL

N=2 NH=2

DX=(Y8-Y7)/NY8

DXB=ZPOLY/(NZPOLY-1

)

DXA= ( ZFRP-ZPOLY) /NZFRP

MB= 1 MA=2

CALL B3 (N, NH, MA , MB,

• DX, DXA, DXB, BO) CALL VB ( BTEMP , NN , I , J , BO

)

C

C - EXTERIOR CORNERS C

1=1 J=1

N1 = 1 N2=3

NH1 =8 NH2=1 DX1=Y5/(NY5-1)

DX2=ZP0LY/(NZP0LY-1)

M=1

A-75
o

o

o

CALL BA (N1,N2,NH1,NH2,M,

• DX1 ,DX2,B0) CALL VB ( BTEHP , NN , I , J , BO

)

C

1=1 JsNCOL

N1=2

N2=3

NH1=2

NH2=1 DX1=(Y8-Y7)/NY8

CALL BA (N1,N2,NH1,NH2,M,

• DX1 ,DX2,B0) CALL WB ( BTEMP , NN , I , J , BO

)

C IrNTZFRP

HI =2 N2=A

NH1=2

HH2=3

DX2=(ZFRP-ZPOLY)/NZFRP

M=2

CALL BA (N1 ,N2,HH1 , NH2,M,

• DX1 f DX2,B0) CALL WB (BTEMP,NN,I,J.BO)

C I=NTZWAL

J=NTYWAL

N1=2

N2=A

NH1=A

NH2=5 DX 1 = (YWALL-Y 1 L ) /WYWALL

DX2=( ZWALL-ZFRP ) /NZWALL

M=5

CALL BA (N1 ,N2,NH1 ,NH2 ,M, • DX1 , DX2,B0) CALL WB ( BTEHP , NN , I , J , BO

)

C J=NTY1L

N1 = 1 N2=A NH1=6

NH2=5

CALL BA (N1,N2,NH1,NH2,M,

• DX1.DX2.B0) CALL WB (BTEMP, NN,I,J, BO)

C IrNTZFRP

J=1

N1 = 1 N2=A NH1=8

NH2=7

DX1=Y5/(NY5-1)

DX2= ( ZFRP-ZPOLY ) /NZFRP

M=2 CALL BA (N1 ,N2,NH1 ,NH2,M,

• DX1.DX2.B0) CALL WB (BTEMP,NN, I, J,BO) c

200 CONTINUE

WRITE BTEMP INTO B VECTOR

DO 300 1=1, NN

300 B(I )=BTEMP(I) RETURN C

END

A-76
cCIHIIIIHIHIIHIIIIIHIIHIIHIIIIIIIHIIIMIIIIIIIIIIHIIIII

c SUBROUTINE QS0LV2 (ISTEP.KOUNT.T.NN.TIME)

CCIIHIIIIIIIIHIIIHHIIIIIIIHIIIIIIIIIHIIliHIIIIIIHIflHIIIH

C» SUBROUTINE TO EFFECT ENERGY BALANCES • CIIIHIHIIHIIIIIfllMHHIIIIIIIIHIIIIIIIItlllMIIIIIIHIIIIIH

c

COMMON /NDMZ/ NZPOLY , NZFRP. NZVALL

COMMON /NUMY/ NY5.NY4.NY1L.NYWALL.NY1R.NY7.NT8

COMMON /GRIDBY/ NTY4.NTY1L.NTYWAL.NTY1R.NTY7.NCOL

COMMON /GRIDBZ/ NTZFRP. NTZWAL COMMON /ZCOOR/ ZPOLY, ZFRP. ZWALL COMMON /YCOOR/ Y5.Y4.Y1L.YWALL.Y1R.Y7-Y8

COMMON /MATS/ RKAY (6 ) .RHO(6 ) ,CP(6

)

COMMON /HTCOEF/ HTC(8) ,TAMB(8),TAMBP(8) COMMON /HEAT/ Q( 21 ) ,IEB, IQPR. ICHECK( 5) C DIMENSION T(NN) t HT(13).N(5)

C CALL HEADER 1 WRITE (6.3001) TIME

C HT(1).Q(1)+Q(2)+Q(3)

HT(2)=Q(4)

HT(3)=Q(5)+Q(6)*Q(7)

HT(4)=Q(8)+Q(9)«-Q(10)

HT(5)=Q( 11)

HT(6)=Q( 12)

HT(7)=Q(13)+Q(14)+Q(15)

HT(8)=Q( 16

)

HT(9)=Q( 17)

HT( 10)sQ( 18) HT (11 )=Q( 19) HT(12)=Q(20)

HT(13)=Q(21)

IF ( I CHECK ( 1 ) . EQ. 0 ) GO TO 100 WRITE (6.2003) (HT(I ) ,1=1 . 13) C

C HEAT BALANCES FOR VARIOUS REGIONS

C

c

C REGION 5-6-4L-4U-1L-2L

C H1=HT(7)-HT(8)+HT(1)-HT(12)

C

C REGION 1C C

100 IF (ICHECK(2).EQ.O) GO TO 200

H2=HT(8)-HT( 11 )+HT(2)-HT(9)

C

C REGION 1 R-2R-7L-7U-8-9

C 200 IF (ICHECK(3).EQ.O) GO TO 300

H3=HT(9)+HT(10)-*HT(4HHT(3)

C

C REGION 2C C

300 IF (ICHECK(4). EQ. 0) GO TO 400

H4=-HT( 13)+HT( 12)+HT( 1 1 )-HT( 10) C

C REGION 3 C 400 IF (ICHECK(5).EQ.O) GO TO 500

H5=HT(6)+HT(5)+HT(13)

C

C TEST FOR PRINTING

C 500 IF (KOUNT.LT. IQPR) RETURN KOUNTsO

C

C PRINT ENERGY BALANCES

C

IF ( ICHECK ( 1 ) . GE. 0 ) GO TO 600

CALL HEADER 1 WRITE (6,2000) H1.H2.H3.H4.H5

CC

A- 7 7
oo

o

o

o

o

o

oooo

oooo

COMPUTE FRAME COEFFICIENT

- AVERAGE MC SPECIMEN SURFACE TEMP

DELTAZ=ZWALL-ZFRP

PER=15. 236 AW=

1

3 . 93 600 QT=PER *HT ( 7 ) +AW*HT ( 6 ) /DELTAZ

PRINT NET ENERGY INTO OR OUT OF MC

WRITE (6.2001) QT

AHTMC=AW»HT( 6 ) /DELTAZ FCMCA = 1

. -ABS ( AHTMC ) /ABS ( QT

)

FRAME COEFFICIENT USING 1 - D HEAT FLOW THROUGH WALL

JsNTYWAL

IsNTZWAL

CALL GETNOD (I.J.N)

TR=T(N(1))

J=NTY1L

CALL GETNOD (I.J.N)

TL=T(N( 1 )

)

HTW=AW • ( TR-TL ) *RKAY ( 5 ) / (YWALL-Y 1 L

)

FCHXS=1. -ABS(HTW )/ ABS(QT)

PRINT RELATIVE FRAME COEFICIENTS

WRITE (6.2002) HTW. FCMCA. FCHXS

CALL HEADER

1

RETURN

2000 FORMAT (//

• / 20X."E NERGY BALA N C E S"//

• / 15X."H1 (REGION 5-6-4U-4L- 1L-2L) =",E1H.7/

• / 15X. "H2 (REGION 1C ) =",EH*.7/

• / 15X."H3 (REGION 1R-2R-7L-7U-8-9) =",E14.7/

• / 15X."H4 (REGION 2C ) r".EH*.7/

• / 15X."H5 (REGION 3 ) =",EU.7/

2001 FORMAT (//

• / 20X. "NET MC ENERGY =",E14.7///)

2002 FORMAT (//

• / 20X. "FRAME COEFFICIENTS-REL"/

• / 15X."HTW =".EU.7/

• / 15X. "FCMCA =".E14.7/

• / 15X, "FCHXS =".E1H.7///)

2003 FORMAT (//

• / 20X, "T OTAL LINE Q ' S"//

• / 15X. "HT (1 ) =",E14.7/

• / 15X. "HT (2 ) =".E14.7/

• / 15X.-HT (3 ) =",E14.7/

• / 15X.-HT (* ) =",EU.7/

• / 15X. "HT (5 ) =".E1H.7/

• / 15X. "HT (6 ) =",E1H.7/

• / 15X. "HT (7 ) =".EU.7/

* / 15X. "HT (8 ) =".E14.7/

• / 1 5X . "HT (9 ) =".E14.7/

• / 15X."HT (10) =",E14.7/

• / 1 5X . "HT (11) r'.E14.7/

• / 15X. "HT (12) s",E14.7/

• / 15X, "HT (13) =".E14.7///)

3001 FORMAT (//

• / 20X. "E NERGY BALA N C E S A T"/

• / 20X, " TIME = " .E14.7///)

END

A-78
c SUBROUTINE QCONVH (NH f I t NC1 f NC2,DX,T,NN,F)

CCIIIIIIIIHIIIIIIIUIIIIHIIIIIIHIIIIHHHIIIIIIHHI

C« SUBROUTINE TO COMPUTE HORIZONTAL CONVECTIVE Q'S

CHIMI<HIHIIIIIIIIIHIIII<IIIIMIMIIIIIIIIM<HIIII

C

COMMON /HTCOEF/ HTC( 8) ,TAMB( 8) , TAMBP(8)

COMMON /TIM/ TIMO.TIMF, DELTA, ALFA .IRFINE C

DIMENSION T(NN),N(5)

C

TI = ALFA »TAMBP ( NH )( 1 . -ALFA )*TAMB ( NH

)

NDIV»NC2-NC1

IF (NDIV.GE. 1) GO TO 100 CALL HEADER 1 WRITE (6,3000) I, NCI ,NC2,NDIV CALL ADIOS (2)

C

100 J1=NCH1

J2=NC2-1

SUM=0. C

C FIRST NODE C

JrNCI CALL GETNOD (I,J,N)

SUM=SUM+0. 5*(TI-T(N( 1 ) )

)

C

C LAST NODE C

J=NC2

CALL GETNOD (I,J,N)

SUM=SUM+0. 5*(TI-T(N( 1 ) )

)

C

C INTERIOR NODES C

DO 200 J=J1,J2

CALL GETNOD (I,J,N)

200 SUM=SUM+(TI-T(N( 1 ) )

)

C FrHTC(NH)»DX*SUM

C

RETURN C

3000 FORMAT (//15X,« ERROR TERMINATION /

• / 20X, "FROM QCONVH - NDIV=0"/

• / 20X , "I , NC 1

, NC2 , NDIV r »,5I5) C

END

A-79
non

oo

o

o

o

o

o

oo

ooooo

ooo

SOBRODTINE QCONVV (NH.NR1 ,NR2, J,DX,T,NN,F)

SUBROUTINE TO COMPUTE HORIZONTAL CONVECTIVE Q'S

COMMON /HTCOEF/ HTC(8) ,TAMB(8) ,TAMBP(8)

COMMON /TIM/ TIMO.TIMF, DELTA, ALFA, IRFINE

DIMENSION T(NN),N(5)

TI =ALFA*TAMBP ( NH ) ( 1

. -ALFA ) *TAMB ( NH

)

NDIV=NR2-NR1

IF (NDIV.GE. 1) GO TO 100 CALL HEADER 1 WRITE (6,3000) NR1,NR2,J,NDIV

CALL ADIOS (2)

100 I1=NR1+1

I2=NR2-1

SUM=0.

FIRST NODE

I=NR1 CALL GETNOD (I,J,N)

SUM=SUM+0. 5*(TI-T(N( 1 ) )

)

LAST NODE

I=NR2

CALL GETNOD (I,J,N)

SUM=SUM+0. 5*(TI-T(N( 1 ) )

)

INTERIOR NODES

DO 200 1=11,12

CALL GETNOD (I,J,N)

200 SUM=SUM+(TI-T(N( 1 ) )

)

C F=HTC(NH)«DX*SUM

C

RETURN

C3000 FORMAT (//15X," ERROR TERMINATION "/

• / 20X, "FROM QCONVV - NDIV=0"/

• / 20X , " NR 1 , NR2 , J , NDIV = ",5I5) C

END

A-80
cCHIHIIIHIIHmilHtlllllllllllllllHIHIIIIIIIIIIII

C SUBROUTINE QCONDV (M,NR1 ,NR2, J,ND,DX,DY, C,T,NN,F)

CCIHIHIIimtllllllllHHIfltlllllllllllllHIfllHIHI

C« SUBROUTINE TO COMPUTE VERTICAL CONDUCTION FLUX CltlllHIHHHIHIIHHIIIIIIIIIIIHIIIHIIHIHIIIIII

c

COMMON /MATS/ RJCAY (6) ,RH0(6 ) ,CP(6

)

COMMON /DBUG/ IDBUG(IO)

C DIMENSION T(NN),N(5) C

NDIV=NR2-NR1

IF (NDIV.GE. 1 ) GO TO 100 CALL HEADER 1 WRITE (6,3000) NR1 ,NR2, J.NDIV

CALL ADI0S (2) 100 IlsNRI+1

I2=NR2-1

SUM=0. c

C FIRST POINT

C IsNRI CALL GETN0D (I,J,N)

SUM=SUM+0. 5*(T(N(1))-T(N(1 )+ND)

)

C

C LAST NODE C

I=NR2

CALL GETNOD (I,J,N)

SUM=SUM+0. 5»(T(N( 1 ) )-T(N( 1 )+ND)

)

C

C INTERIOR NODES

C

DO 200 1:11,12

CALL GETNOD (I,J.N)

200 SUM=SUM+T(N( 1 ) )-T(N( 1 )+ND)

C F=C#RKAY(M)*DI*SUM/(DX)

C

C

C WRITE DEBUG INFO

C

IF ( IDBUG ( 6 ) . EQ. 0 ) RETURN WRITE (6,2000) M,NR1 ,NR2, J.F

C

RETURN

C2000 FORMAT (//

• / 20X,"IN QCONDV - DEBUG INFO"/

• / 15X,"M,NR1,NR2,J,F = ",4I5,5X,EU.7//)

3000 FORMAT (//15X,» ERROR TERMINATION "/

• / 20X, "IN QCONDV - NDIVrO"/

• / 20X, "NR 1.NR2, J.NDIV * ",5I5) C

END

A-81
oono

ooo

ooo

ooo

o

o

ooooo

C SUBROUTINE QCONDH (M,I,NC1,NC2,ND,DX,DY,C,T,NN,F)

SUBROUTINE TO COMPUTE HORIZONTAL CONDUCTION FLUX

COMMON /MATS/ RKAY(6) ,RH0(6) ,CP(6)

COMMON /DBUG/ IDBUG(IO)

DIMENSION T(NN) ,N(5)

NDIV=NC2-NC1

IF (NDIV.GE. 1) GO TO 100 CALL HEADER 1 WRITE (6,3000) I,NC1,NC2,NDIV

CALL ADI0S (2) 100 J1=NC1+1

J2=NC2-1

SUM=0.

FIRST POINT

J=NC 1 CALL GETN0D (I,J,N)

SUM=SUM+0. 5*(T(N( 1 ) )-T(N( 1 )+ND)

)

LAST NODE

J=NC2

CALL GETNOD (I,J,N)

SUM=SUM+0. 5*(T(N( 1 ) )-T(N( 1 )+ND)

)

INTERIOR NODES

DO 200 J=J1,J2

CALL GETNOD (I,J,N)

200 SUM=SUM+T(N( 1 )

) —T ( N ( 1 )+ND)

F=C»RKAY(M)»DY*SUM/(DX)

WRITE DEBUG INFO

IF (IDBUG(8).EQ.0) RETURN WRITE (6,2000) M,I,NC1,NC2,F

2000 FORMAT (//

• / 20X, "IN QCONDH - DEBUG INFO"/

• / 1 5X , "M, I , NC 1 , NC2 , F = \4I5,5X,E14.7//)

C RETURN

C3000 FORMAT (//15X," ERROR TERMINATION "/

• / 20X, "IN QCONDH - NDIV=0"/

• / 20X , "I , NC 1 , NC2 , NDIV = ",5I5)

G END

A-82
SUBROUTINE PLOTH (IND,TIME,T,NN)

SUBROOTIHE TO PERFORM CONTOUR PLOTS INPUT

COMMON /PLOTS/ IPLT,NUMCON,NUMHIS,IOPT,NOOT,RANGE(4),NSTEPH,

• ICHAR (10), NSTEPC

COMMON /PHIST/ ITITH(5,72) ,ITITIH(5,36) , ITITYH ( 5 , 36 ) ,MH(5)

,

• NODESH ( 5 , 5 ) , N , XMI N , ZMA2

DIMENSION T(NN)

DIMENSION ITITLE(l44),YH(8l,5,5),Y(8l,5),X(8l),IMAG4(5151)

DATA INC, IT /I, 81/

EH=KH«-1 X(KH)*TIME

FILL PLOT ARRAY

DO 100 1=1 , NUMCON

M=MH(I

)

DO 100 J= 1 ,M NODEsNODESH ( I , J

)

100 YH(KH,J,I)=T( NODE

)

IF (IND.EQ.O) RETURN

PERFORM PLOTTING

CALL HEADER 1 WRITE (6,4000) TIME

DO 400 1=1 .NUMHIS MsMH(I)

FILL TITLE

L=0

DO 210 Jsl ,72 L=L+1 210 ITITLE ( L ) =ITITH ( I , J

)

DO 220 J=1,36

L=L+1 220 ITITLE ( L ) =ITITXH ( I , J

)

DO 230 J=1,36

LsL+1

230 ITITLE ( L ) sITITTH ( I , J

)

FILL THIS PLOT ARRAY

DO 300 Lcl.N

DO 300 J=1 ,M 300 Y(J,L)=YH(L,J,I)

A-83
SET RANGE

RANGE ( 1 )sXMIN RANGE (2 ) eXMAX

PERFORM PLOT

IF (I.NE. 1) CALL HEADER

1

WRITE (6,2000) I DO 395 J=1,M

395 WRITE (6,2001) J.NODESH(I, J)

400 CALL OSPLT (X,Y,IY,N,KH, INC, ITITLE, RANGE, ICHAR.IOPT,

• IMAG4 ,IER)

RETURN

2000 FORMAT (//25X,"C URVE OUTPUT"/

• / 15X, "CURVE NO. ",I5//

• / 10X, "LABEL", 5X, "NODE"/)

2001 FORMAT ( • 5(/7X,I5, 10X,I5)

)

4000 FORMAT (////15X,"T IME HISTORY PLOTS

• / 30X,"T IME =",E12.5////) END

A T"/

A-84
ooo

ooo

ooo

ooo

o

on

non

ooooo

C

SUBROUTINE PLOTIN (ITRAN)

SUBROUTINE TO PERFORM BASE PLOT INPUT

COMMON /PLOTS/ IPLT , NUMCON , NUMHIS , IOPT , NOUT , RANGE ( 4 ) , NSTEPH

,

• ICHAR (10), NSTEPC

READ PLOT PARAMETERS

READ (5,1000) NUMCON, NUMHIS, IOPT, NOUT, NSTEPC

READ (5,1001) THIN, THAI CALL HEADER 1 WRITE (6,2000) NUMCON, NUMHIS, IOPT, NOUT, NSTEPC, TMIN, THAI

IF (NOUT. EQ. 2. OR. NOUT. EQ. 6) GO TO 100 WRITE (6,3000) WRITE (6,3001) NOUT CALL ADIOS (2)

100 CALL UGETIO (3,NIN,NOUT)

IF ( NUMCON. LE. 5) GO TO 200 WRITE (6,3000) WRITE (6,3002) NUMCON

200 IF (NUMHIS. LE. 5) GO TO 300 WRITE (6,3000) WRITE (6,3003) NUMHIS

READ CONTOUR INFORMATION

300 IF (NUMCON. NE.O) CALL PLTCIN (NUMCON)

READ TIME HISTORY INPUT

IF (NUMHIS. NE.O) CALL PLTHIN (ITRAN)

SET CURVE LABELS

ICHAR ( 1

) = 1

H

SET RANGE

RANGE(3)=TMIN

RANGE(4)*TMAX

A-85
RETURN

1000

1001 2000

FORMAT

FORMAT

FORMAT

3000

3001 3002

3003

20X,

15X,

15X, 15X f 25X,

25X,

15X,

25X,

25X,

15X,

15X,

25X,

25X,

(1615)

(8E10.0)

(//

"PLOT

"NUMCON "NUMHIS

"IOPT

" 80 CHARACTERS

"126 CHARACTERS

"NOUT

"FILE IMSLS. OUT

"FILE FLANK. OUT

BASE INFORMATIO N"//

=",15/

=",15/

=",15/ IOPT

IOPT

0"/

1 "/

',15/

FORMAT

FORMAT

FORMAT

FORMAT

NOUT =2"/

NOUT = 5"/ "STEP INTERVAL FOR CONTOUR PLOTS =",15/ "RANGE FOR PLOTS"/

"TMIN r",ElH.5/

"TMAX r",E1i».5//>

(//" ERROR TERMINATION FROM PLOTIN"//)

(15X,"NOUT = ",I5," SHOULD BE 2 OR 6")

(15X, "NUMCON = ",I5," SHOULD BE < OR = 5")

(15X, "NUMHIS = ",I5," SHOULD BE < OR = 5")

END

A-86
CIIIHIIHIIIIIIIIIIHIItllllllllHIIIIIIIIIIIIIIIIIIIIIIIItIHII

C SUBROUTINE PLOTC (IND.TIME.T.NN)

CCHIHIIUIIIIIIIIIIHIIIIIIIIHIIIIIIHIIIIHIIUIIIIIIinilllH

C* SUBROUTINE TO PERFORM CONTOUR PLOTS INPUT CHIHIIIIIIHfllHIHIIIIIIIIIIIHIHIHHHItllHHIHHIIIIIH

c COMMON /DBUG/ IDBUG(IO) COMMON /PLOTS/ IPLT , NUMCON , NUMHIS , IOPT , NOUT , RANGE ( 4 )

, NSTEPH

,

• ICHAR( 10) |NSTEPC

COMMON /PCONT/ ITITC(5,72),ITITXC(5,36),ITITYC(5,36),NC(5),

• N0DESC(5,20)

C DIMENSION T(NN)

DIMENSION ITITLE(144),YC(20,5,5),Y(20,5),X(20),

• IMAG4(5151),TIMSAV(5)

C

DATA INC.IYC /I, 5/

C

100 KCsKC+1

IF (IND.EQ. 1.0R.KC.EQ.6) GO TO 300

TIMSAV ( KC ) =TIME

C

C FILL PLOT ARRAY

C

DO 200 Is 1 .NUMCON N=NC(I

)

DO 200 J=1,N

NODEsNODESC(I.J)

200 YC ( J , KC , I ) =T ( NODE

)

C RETURN C

C PERFORM PLOTTING

C

300 KC=KC-1

CALL HEADER

1

C WRITE (6,4000) KC , TIME

DO 600 1=1, NUMCON

IF (I.NE. 1) CALL HEADER

1

WRITE (6,2000) I DO 305 J=1,KC

305 WRITE (6,2001) J.TIMSAV(J)

C

C FILL TITLE

C L=0

DO 310 J=1,72

L=L+1

310 ITITLE ( L ) =ITITC ( I , J

)

DO 320 J=1,36

L=U1

320 ITITLE ( L ) sITITXC ( I , J

)

DO 330 J=1,36

L=U1

330 ITITLE ( L ) sITITYC ( I , J

)

A-87
C FILL THIS PLOT ARRAY

C

N=NC(I

)

DO 400 J=1 ,N X ( J ) =NODESC ( I , J

)

DO 400 L=1,EC

400 Y(J,EC)=YC(J,EC,I)

C

C PICE OOT MAX, MIN C XMIN=1.E50

DO 500 J=1,N

500 XMINsAMINI (XMIN,X( J )

)

C XMAX=-1.E50

DO 510 Jsl.N

510 XMAX=AMAX1(XMAX,X(J))

C RANGE ( 1

) sXMIN RANGE(2)=XMAX

IF (IDBDG(IO).EQ.O) GO TO 520 WRITE (6,5000) I, TIME WRITE (6,5001) X WRITE (6,5002) Y WRITE (6,5003) IYC,N,KC,INC WRITE (6,5004) RANGE WRITE (6,5005) IOPT 520 CONTINUE

C

C PERFORM PLOT

C600 CALL USPLT (X, Y,IYC,N, EC, INC, ITITLE, RANGE, ICHAR, IOPT,

• IMAG4.IER)

C EC=0

IF (IND.NE. 1) GO TO 100 C RETURN

C2000 FORMAT (//25X,"C URVE OUTPUT"/

• / 15X, "CURVE NO. ",I5//

• / 10X, "LABEL", 5X, "TIME"/)

2001 FORMAT ( • 5 ( /7X , 15 , 6X , E 1 2. 5 )

)

4000 FORMAT (15X.I5," NO. OF CONTOUR PLOTS AT"/

• / 30X,"T I M E =",E12. 5////)

5000 FORMAT (5X, "DEBUG ( 10) OUTPUT, PLOT NO.", 15," TIME =",E12.5//)

5001 FORMAT (5X,"X - MATRIX ",/10(/10( 1X,E1 1. 4) )

)

5002 FORMAT (5X,"Y - MATRIX ",/10(/10( IX, El 1. 4) )

)

5003 FORMAT (5X,"IYC,N,EC,INC =",1615)

5004 FORMAT (5X, "RANGE «,5(2X,E12.4))

5005 FORMAT (5X,"IOPT = «,I5) END

A-88
cCiiiiiiiiuiiHiiiiHiiimiiiiiiiniiiiiiiiiiHiiiiiimHiiiiHiii

C SUBROUTINE PLTCIN ( NUMCON

)

CCHIIIflttttflMHIlHIIHIIIIIIIIIHIIIIIIIIIIHIIIMMIIIIIIUIII

C" SUBROUTINE TO PERFORM CONTOUR PLOTS INPUT • CHtllHItHIIIIHIIHHIHHHIIIIIIIIIIHHIIIHIIHIIIIIMHIHI

c

COMMON /PCONT/ ITITC(5,72) ,ITITIC(5, B6 ) ,ITinC{5, 36 ) ,NC(5)

,

• NODESC (5.20)

CALL HEADER

1

WRITE (6(4000)

C

DO 100 1=1 ,NUMC0N WRITE (6,2000) I C

C READ TITLES

C

READ (5,1002) (ITITC(I, J) , J=1 ,72) READ (5,1002) (ITITXCU, J) , J=1 ,36) READ (5,1002) (ITITYC(I, J),J=1,36) WRITE (6,2001) (ITITC(I, J) , J=1 ,72) WRITE (6,2002) (ITITXC(I, J) , J=1 ,36) WRITE (6,2003) (ITITYC(I,J),J=1,36) C

C READ NUMBER OF NODE POINTS C

READ (5,1000) NC(I)

NCI=NC(I

)

IF (NCI. GT. 0. AND. NCI. LE. 20) GO TO 10 WRITE (6,3000) WRITE (6,3001) I, NCI CALL ADIOS (2) C

10 READ (5,1000) (NODESC(I, J) , J=1 ,NCI) WRITE (6,2004) I, NCI, (N0DESC(I, J),J«1,HCI)

C

100 CONTINUE

C RETURN

C

1000 FORMAT (1615)

1002 FORMAT (72A1)

2000 FORMAT (/////25X," CONTOUR GRAPH NO. ",15/

• / 15X,"T I T L E S"/)

2001 PORMAT (/ 15X, "GRAPH TITLE"

• / 5X.72A1)

2002 FORMAT (/20X, "X-AXIS TITLE"/ 1 OX, 72A 1

)

2003 FORMAT (/20X, "Y-AXIS TITLE"/10X,72A1

)

2004 FORMAT (/

• 15X, "GRAPH NO. ",I5,5X,"NO. OF POINTS ",I5

• / 15X, "POINTS =",20I5)

3000 FORMAT (15X," ERROR TERMINATION FROM PLTCIN ") 3001 FORMAT (/15X,"NO. OF POINTS IS WRONG, GRAPH NO. «",I5/

• / 30X, "NO OF POINTS = " , 15 , 5X-SH0ULD BE > 0 AND < OR * 20")

4000 FORMAT (/20X,"C ONTOUR PLOT INPUT ECH 0*//) C

RETURN

C END

At-89
EXAMPLE PROBLEM 1: FLANKING LOSS - STATIC CASE

FLANK LOSS POLYSTYRENE B-B 9-01 -82 TEMPB--STATIC CASE

2 -1 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 -1

-1 -1 -1 -1 -1

4 3 4 5 4 3 4 .254 .273 .400 .502 .578 0.610 .864

5 3 5 .456 .457 .957

2.0 0.0 7.45 7.45 0.0 2.33 2. 33 25.0 100. 65.0 65.0 100- 38. 38.

0. 022 0.430 0.180 0.180 .036

32. 00 1920. 1100. 1100. 15.10

921.1 963-0 2000. 2000. 1214.

1 38.0 2 38.2

0.0 3600. 900. 1.0 25.0 110 6 1 20. 80. TEST PLOT FOR CONTOURS - RJO 09/01/82

Z - AXIS NODES

Y - AXIS TEMPERATURE (C)

20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

17 18 19 20

TEST PLOT - TIME HISTORIES 08/27/82 BY PJB

X - AXIS TIME (SEC) Y - AXIS TEMPERATURE (C)

5

1 2 3 4 5

A-90
1

FLAK LOSS POLTSTTREHE B-B 9-01-82 TEMPB-STATIC CASE

CONTEOL DATA

NUMTBC

NUMIT

IDATA

IPR

IRPINE

IDGT

IEB

IQPR

ITRAN

IPLT

2

-1

0

0

0

DEBOG PARAMETERS

1 (WRITE A.C MATRICES. B VECTOR)

2 (WRITE DECOMPOSED C MATRIX )

3 (WRITE NODAL ASSOC. )

4 (WRITE ELEMENT MATRICES )

5 (WRITE INTERIOR DATA )

6 (WRITE C BEFORE DECOMP. )

7 (WRITE B VECTOR )

8 ( )

9 ( )

10 ( )

0

0

0

0

0

0

0

0

0

-1

A-91
ENERGY BALANCES

ICHECK(I)

I CHECK (2)

ICHECK(3)

ICHECK (4)

ICHECK (5)

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

GRID SPECIFICATION

Y NUMBERS FOLLOW

NY5 = 4

NY 4 = 3

NY1L r 4

NYWALL S 5

NY1R * 4

NY7 = 3

NY8 = 4

Y’S FOLLOW

Y5 = 0.2540E+00

Y4 S 0.2730E+00

Y1L = 0. 4000E+00

YWALL = 0. 5020E+00

Y1R = 0. 5780E+00

17 = 0. 6 1 00E+00

Y8 = 0. 8640E+00

Z NUMBERS FOLLOW

NZPOLY = 5

NZFRP * 3

NZWALL = 5

Z'S FOLLOW

ZPOLY = 0. 4560000E+00

ZFRP = 0. 4570000E+00

ZWALL — 0. 9570000E+00

i A- 9 2
FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

CUMULATIVE NODE NU

Y - DIRECTION NODE NUMBERS

NTY4

NTI1L

NTYUAL

NTY1R

NTY7

NCOL

Z - DIRECTION NODE NUMBERS

NTZFRP

NTZWAL

7

11

16

20

23

27

8

13

TOTAL NUMBER OF NODES (EQUATIONS) >

LOWER NUMBER OF NODE POINTS «

M B E R S

246

216

A-93
FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

CONVECTIVE COEFFICIENTS

HEAT TRANSFER COEFFICIENTS FOLLOW

0. 20000E+01 O.OOOOOE+OO 0.74500E+01 0.74500E+01

AMBIENT TEMPERATURES FOLLOW

0. OOOOOE+OO 0.23300E+01 0.23300E+01 0. OOOOOE+OO

0. 25000E+02 0. 10000E+03 0.65000E+02 0.65000E+02

1

0. 10000E+03 0. 38000E+02 0. 38000E+02 0. 10000E+03

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

MATERIAL PROPERTIES

I RKAY RHO CP

1

1 0. 2200E-01

2 0. 4300E+00

3 0. 1800E+00

4 0. 1800E+00

5 0.3600E-01

6 0. OOOOE+OO

0. 3200E+02

0. 1920E+04

0. 1100E+04

0. 1100E+04

0. 1510E+02

0. 0000E+00

0. 9211E+03

0. 9630E+03

0. 2000E+04

0. 2000E+04

0. 1214E+04

0. OOOOE+OO

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

1 TEMPERATURE BOUNDARY COND

NUMBER OF TEMPERATURE BCS = 2

1

0. 380E+02 2 0. 382E+02

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE
INITIAL TEMPERATURES

1 INITIAL TEMPERATURE 0. OOOOOE-f 00

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

INPUT PHASE COMPLETE

1

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

EQUATION DATA

DETERMINANT

ACCURACY TEST PARAMETER

NUMBER OF EQUATIONS

ROW SIZE IN CALLING ROUTINE

ERROR PARAMETER ( 129-SIN. ,34-ACC.

)

0.1 3749 133*76

0. 192466 1E+06

246

246

0

A- 95
FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

H

TEMPER A TOR E STATE

AT T I M E = O.OOOOOE+OO

1 0. OOOOOOOE+OO 51 0. OOOOOOOE+OO 101 0. OOOOOOOE+OO 151 0. OOOOOOOE+OO 201 0. OOOOOOOE+OO

2 0. OOOOOOOE+OO 52 0. OOOOOOOE+OO 102 0. OOOOOOOE+OO 152 0. OOOOOOOE+OO 202 0. OOOOOOOE+OO

3 0. OOOOOOOE+OO 53 0. OOOOOOOE+OO 103 0. OOOOOOOE+OO 153 0. OOOOOOOE+OO 203 0. OOOOOOOE+OO

4 0. OOOOOOOE+OO 54 0. OOOOOOOE+OO 104 0. OOOOOOOE+OO 154 0. OOOOOOOE+OO 204 0. OOOOOOOE+OO

5 0. OOOOOOOE+OO 55 0. OOOOOOOE+OO 105 0- OOOOOOOE+OO 155 0. OOOOOOOE+OO 205 0. OOOOOOOE+OO

6 0. OOOOOOOE+OO 56 0. OOOOOOOE+OO 106 0. OOOOOOOE+OO 156 0. OOOOOOOE+OO 206 0. OOOOOOOE+OO

7 0. OOOOOOOE+OO 57 0. OOOOOOOE+OO 107 0. OOOOOOOE+OO 157 0. OOOOOOOE+OO 207 0. OOOOOOOE+OO

8 0. OOOOOOOE+OO 58 0. OOOOOOOE+OO 108 0. OOOOOOOE+OO 158 0- OOOOOOOE+OO 208 0. OOOOOOOE+OO

9 0. OOOOOOOE+OO 59 0. OOOOOOOE+OO 109 0. OOOOOOOE+OO 159 0. OOOOOOOE+OO 209 0. OOOOOOOE+OO

10 0. OOOOOOOE+OO 60 0. OOOOOOOE+OO 110 0. OOOOOOOE+OO 160 0. OOOOOOOE+OO 210 0. OOOOOOOE+OO

11 0. OOOOOOOE+OO 61 0. OOOOOOOE+OO 111 0. OOOOOOOE+OO 161 0. OOOOOOOE+OO 211 0. OOOOOOOE+OO

12 0. OOOOOOOE+OO 62 0. OOOOOOOE+OO 112 0. OOOOOOOE+OO 162 0. OOOOOOOE+OO 212 0. OOOOOOOE+OO

13 0. OOOOOOOE+OO 63 0. OOOOOOOE+OO 113 0. OOOOOOOE+OO 163 0. OOOOOOOE+OO 213 0. OOOOOOOE+OO

14 0. OOOOOOOE+OO 64 0. OOOOOOOE+OO 114 0. OOOOOOOE+OO 164 0. OOOOOOOE+OO 214 0. OOOOOOOE+OO

15 0. OOOOOOOE+OO 65 0. OOOOOOOE+OO 115 0. OOOOOOOE+OO 165 0. OOOOOOOE+OO 215 0. OOOOOOOE+OO

16 0. OOOOOOOE+OO 66 0. OOOOOOOE+OO 116 0. OOOOOOOE+OO 166 0. OOOOOOOE+OO 216 0. OOOOOOOE+OO

17 0. OOOOOOOE+OO 67 0. OOOOOOOE+OO 117 0. OOOOOOOE+OO 167 0. OOOOOOOE+OO 217 0- OOOOOOOE+OO

18 0. OOOOOOOE+OO 68 0. OOOOOOOE+OO 118 0. OOOOOOOE+OO 168 0. OOOOOOOE+OO 218 0. OOOOOOOE+OO

19 0. OOOOOOOE+OO 69 0. OOOOOOOE+OO 119 0. OOOOOOOE+OO 169 0. OOOOOOOE+OO 219 0. OOOOOOOE+OO

20 0. OOOOOOOE+OO 70 0. OOOOOOOE+OO 120 0. OOOOOOOE+OO 170 0- OOOOOOOE+OO 220 0. OOOOOOOE+OO

21 0. OOOOOOOE+OO 71 0. OOOOOOOE+OO 121 0. OOOOOOOE+OO 171 0. OOOOOOOE+OO 221 0. OOOOOOOE+OO

22 0. OOOOOOOE+OO 72 0. OOOOOOOE+OO 122 0. OOOOOOOE+OO 172 0. OOOOOOOE+OO 222 0. OOOOOOOE+OO

23 0. OOOOOOOE+OO 73 0. OOOOOOOE+OO 123 0. OOOOOOOE+OO 173 0. OOOOOOOE+OO 223 0. OOOOOOOE+OO

24 0. OOOOOOOE+OO 74 0. OOOOOOOE+OO 124 0. OOOOOOOE+OO 174 0. OOOOOOOE+OO 224 0. OOOOOOOE+OO

25 0. OOOOOOOE+OO 75 0. OOOOOOOE+OO 125 0. OOOOOOOE+OO 175 0- OOOOOOOE+OO 225 0. OOOOOOOE+OO

26 0. OOOOOOOE+OO 76 0. OOOOOOOE+OO 126 0. OOOOOOOE+OO 176 0. OOOOOOOE+OO 226 0. OOOOOOOE+OO

27 0. OOOOOOOE+OO 77 0. OOOOOOOE+OO 127 0. OOOOOOOE+OO 177 0. OOOOOOOE+OO 227 0. OOOOOOOE+OO

28 0. OOOOOOOE+OO 78 0. OOOOOOOE+OO 128 0. OOOOOOOE+OO 178 0. OOOOOOOE+OO 228 0. OOOOOOOE+OO

29 0. OOOOOOOE+OO 79 0. OOOOOOOE+OO 129 0. OOOOOOOE+OO 179 0. OOOOOOOE+OO 229 0. OOOOOOOE+OO

30 0. OOOOOOOE+OO 80 0. OOOOOOOE+OO 130 0. OOOOOOOE+OO 180 0. OOOOOOOE+OO 230 0. OOOOOOOE+OO

31 0. OOOOOOOE+OO 81 0. OOOOOOOE+OO 131 0. OOOOOOOE+OO 181 0- OOOOOOOE+OO 231 0. OOOOOOOE+OO

32 0. OOOOOOOE+OO 82 0. OOOOOOOE+OO 132 0. OOOOOOOE+OO 182 0. OOOOOOOE+OO 232 0. OOOOOOOE+OO

33 0. OOOOOOOE+OO 83 0. OOOOOOOE+OO 133 0. OOOOOOOE+OO 183 0. OOOOOOOE+OO 233 0. OOOOOOOE+OO

34 0. OOOOOOOE+OO 84 0. OOOOOOOE+OO 134 0. OOOOOOOE+OO 184 0- OOOOOOOE+OO 234 0- OOOOOOOE+OO

35 0. OOOOOOOE+OO 85 0. OOOOOOOE+OO 135 0. OOOOOOOE+OO 185 0- OOOOOOOE+OO 235 0. OOOOOOOE+OO

36 0. OOOOOOOE+OO 86 0. OOOOOOOE+OO 136 0. OOOOOOOE+OO 186 0. OOOOOOOE+OO 236 0. OOOOOOOE+OO

37 0. OOOOOOOE+OO 87 0. OOOOOOOE+OO 137 0. OOOOOOOE+OO 187 0. OOOOOOOE+OO 237 0. OOOOOOOE+OO

38 0. OOOOOOOE+OO 68 0. OOOOOOOE+OO 138 0. OOOOOOOE+OO 188 0. OOOOOOOE+OO 238 0. OOOOOOOE+OO

39 0. OOOOOOOE+OO 89 0. OOOOOOOE+OO 139 0. OOOOOOOE+OO 189 0. OOOOOOOE+OO 239 0- OOOOOOOE+OO

40 0. OOOOOOOE+OO 90 0. OOOOOOOE+OO 140 0. OOOOOOOE+OO 190 0. OOOOOOOE+OO 240 0. OOOOOOOE+OO

41 0. OOOOOOOE+OO 91 0. OOOOOOOE+OO 141 0. OOOOOOOE+OO 191 0- OOOOOOOE+OO 241 0. OOOOOOOE+OO

42 0. OOOOOOOE+OO 92 0. OOOOOOOE+OO 142 0- OOOOOOOE+OO 192 0. OOOOOOOE+OO 242 0. OOOOOOOE+OO

43 0. OOOOOOOE+OO 93 0. OOOOOOOE+OO 143 0. OOOOOOOE+OO 193 0. OOOOOOOE+OO 243 0. OOOOOOOE+OO

44 0. OOOOOOOE+OO 94 0. OOOOOOOE+OO 144 0. OOOOOOOE+OO 194 0. OOOOOOOE+OO 244 0. OOOOOOOE+OO

45 0. OOOOOOOE+OO 95 0. OOOOOOOE+OO 145 0. OOOOOOOE+OO 195 0- OOOOOOOE+OO 245 0. OOOOOOOE+OO

46 0. OOOOOOOE+OO 96 0. OOOOOOOE+OO 146 0. OOOOOOOE+OO 196 0. OOOOOOOE+OO 246 0. OOOOOOOE+OO

47 0. OOOOOOOE+OO 97 0. OOOOOOOE+OO 147 0. OOOOOOOE+OO 197 0. OOOOOOOE+OO

48 0. OOOOOOOE+OO 98 0. OOOOOOOE+OO 148 0. OOOOOOOE+OO 198 0. OOOOOOOE+OO

49 0. OOOOOOOE+OO 99 0. OOOOOOOE+OO 149 0. OOOOOOOE+OO 199 0. OOOOOOOE+OO

50 0. OOOOOOOE+OO 100 0. OOOOOOOE+OO 150 0. OOOOOOOE+OO 200 0. OOOOOOOE+OO

l

A-96

!
run loss polystyrene b-b 9-01-62 tkmpb-static cask

TEMPERATURE STATE

AT TIME * 0. 10000E+01

1 0. 3800000E+02 51 0.3522 16 0E+02 101 0. 53 1 2222E+02 151 0.6328056E+02 201 0.4586673E+02

2 0. 382000 1E+02 52 0. 3524992E+02 102 0.5316809E+02 152 0. 6422050E+02 202 0.5002994E+02

3 0. 2653259E+02 53 0. 35283 10E+02 103 0. 5320329E+02 153 0.6447458E+02 203 0. 5428506 E+02

4 0.26 186 10E+02 54 0.35296 5 1E+02 104 0. 5322792E+02 154 0.6441223E+02 204 0.58685078+02

5 0. 2621045E+02 55 0. 3572284E+02 105 0.5394386E+02 155 0.6389954E+02 205 0.6330580E+02

6 0.2621477 E+02 56 0. 3547023E+02 106 0. 5432954E+02 156 0. 6385086 E+02 206 0.6422902E+02

7 0.261 99 13E+02 57 0. 34833 16E+02 107 0.5451755E+02 157 0.6388768E+02 207 0.6448042E+02

8 0.2573653E+02 58 0. 3465353E+02 108 0.5457301E+02 158 0.6401218E+02 208 0.6441872E+02

9 0. 2562080E+02 59 0. 34657 17E+02 109 0. 3790329E+02 159 0.6469182E+02 209 0.6391583E+02

10 0. 256 1020E+02 60 0.3466478E+02 110 0. 3789857E+02 160 0.6473489E+02 210 0.6388184E+02

11 0.2563477E+02 61 0.346 76 39E+02 111 0. 3789522E+02 161 0.6474092E+02 211 0.6391768E+02

12 0. 2565886E+02 62 0. 35286 27E+02 112 0. 3787305E+02 162 0-6473187E+02 212 0.6 40245 1E+02

13 0.2568828E+02 63 0. 360094 1E+02 113 0. 3787134E+02 163 0.3790366E+02 213 0.6469533E+02

14 0.2572438E+02 64 0. 368521 1E+02 114 0.3787 346 E+02 164 0. 3789894E+02 214 0.6473795E+02

15 0. 257721 9E+02 65 0. 3780980E+02 115 0. 3787935E+02 165 0. 3789558E+02 215 0.6474391E+02

16 0. 258442 1E+02 66 0. 38474 17E+02 116 0. 3808525E+02 166 0. 3787350E+02 216 0.6473492E+02

17 0. 2595937E+02 67 0.39l6658Ef02 117 0. 3842395E+02 167 0. 3787259E+02 217 0.4 1400948+02

18 0.2617430E+02 68 0. 3987474E+02 118 0. 3927539E+02 168 0.3787469E+02 218 0.45884068+02

19 0.266 048 1E+02 69 0. 4058566E+02 119 0.41 757 02E+02 169 0. 3787997E+02 219 0.5036676E+02

20 0. 2750252E+02 70 0.41 287 00E+02 120 0. 4586623E+02 170 0. 3808499E+02 220 0.54861478+02

21 0. 2755939E+02 71 0.4192239E+02 121 0.5002855E+02 171 0. 3842264E+0

2

221 0. 5937775E+02

22 0.2753908E+02 72 0. 425356 0E+02 122 0. 5428287E+02 172 0.3927151E+02 222 0.6392088E+02

23 0.2744057E+02 73 0.4312404E+02 123 0. 58681 93E+02 173 0.4 1746 57E+02 223 0.41 4093 9E+02

24 0.2607216E+02 74 0.4368668E+02 124 0. 6326938E+02 174 0. 45866 59E+02 224 0.45911078+02

25 0. 2592133E+02 75 0. 4372028E+02 125 0.6421651E+02 175 0.5002953E+02 225 0.5041 36 4E+02

26 0. 2590688E+02 76 0. 437457 1E+02 126 0.6447176E+02 176 0. 542844 1E+02 226 0. 5491788E+02

27 0.25 906 27E+02 77 0.4376297E+02 127 0.6440915E+02 177 0. 586841 1E+02 227 0. 59424 15E+02 28 0. 3561974E+02 78 0. 4426237E+02 128 0.6389107E+02 178 0. 632926 9E+02 228 0- 639321 3E+02

29 0. 3493483E+02 79 0.4456305E+02 129 0.6383527E+02 179 0.6 42246 7E+02 229 0. 4 1 4 1 1 29E+02

30 0.3215359E+02 80 0.4472310E+02 130 0. 6387257E+02 180 0.6447746E+02 230 0.4591534E+02

31 0. 3099481E+02 81 0.4477321E+02 131 0.6 4006 00E+02 181 0.6441542E+02 231 0.5041951E+02

32 0. 3098889E+02 82 0. 3674188E+02 132 0.6469005E+02 182 0.6390779E+02 232 0. 5492386E+02

33 0- 30986 37E+02 83 0. 367026 1E+02 133 0.6473334E+02 183 0. 6386638E+02 233 0.59428368+02

34 0.3098721E+02 64 0. 3668343E+02 134 0.6473940E+02 184 0. 639027 1E+02 234 0.6393298E+02

35 0. 31 13385E+02 85 0.3701643E+02 135 0.6473035E+02 185 0.6401834E+02 235 0. 4141 151E+02

36 0.3 1377 06 E+02 86 0. 3702560E+02 136 0. 3790347E+02 186 0.6469357E+02 236 0.459 1586 E+02

37 0.3170750E+02 87 0. 3703944E+02 137 0. 3789876E+02 187 0.6473643E+02 237 0. 5042020E+02

38 0.3211189E+02 88 0.37 0580 1&f02 138 0. 3789540E+02 188 0.6474242E+02 238 0.5492451E+02

39 0. 3240359E+02 89 0. 379791 1 Ei-02 139 0. 3787328E+02 189 0.6473340E+02 239 0. 5942880E+02

40 0.3271690E+02 90 0. 3910089E+02 140 0. 3787199E+02 190 0. 3790384E+02 240 0.6393307E+02

41 0.3304875E+02 91 0.4051501E+02 141 0.3787 4 10E+02 191 0. 3789912E+02 241 0. 4141 150E+02

42 0.3339658E+02 92 0.4230946E+02 142 0. 3787968E+02 192 0. 3789578E+02 242 0.4591582E+02

*3 0.337 58 35E+02 93 0. 4366922E+02 143 0. 38085 13E+02 193 0. 3787375E+02 243 0. 5042014E+02

44 0.3410638E+02 94 0.4512498Ef02 144 0. 3842332E+02 194 0-3787 3 15E+02 244 0.5492444E+02

45 0.3446364E+02 95 0. 466 1453E+02 145 0. 3927354E+02 195 0.3787524E+02 245 0.59'12876E+02

46 0. 34826 92E+02 96 0. 4807436E+02 146 0.4175206 E+02 196 0. 3788023E+02 246 0.6393306E+02

47 0. 351881 1E+02 97 0. 49434332+02 147 0. 45866 42E+02 197 0. 3806484E+02

48 0. 3520857E+02 98 0. 50537 16E+02 148 0.5002907E+02 198 0.3842191E+02

49 0. 35221 48E+02 99 0.5149931E+02 149 0.5428368E+02 199 0.3926932E+02

50 0. 35227 02E+02 100 0. 52350 11E+02 150 0.5868307E+02 200 0.4174054E+02

A-97
FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

ENERGY BA LANCES AT

TIME = 0. 1000000E+01

TOTAL LIN E Q • S

HT (1 ) =-0. 392480 1E+01

HT (2 ) =-0. 1461942E+00

HT (3 ) =-0.9480131E+00

HT (4 ) = 0. 1400995E+01

HT (5 ) = 0.4217818E+01

HT (6 ) =-0.4009733E+01

HT (7 ) s-0. 1973910E+00

HT (8 ) =-0. 538786 9E+00

HT (9 ) »-0.5843179E+00

HT (10) =-0.9703678E-01

HT (11) =-0. 1521804E+00

HT (12) =-0. 8679169E-01

HT (13) =-0. 1501728E+00

1

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

ENERGY BALANCES

HI (REGION 5-6-40-4L-1L-2L)

H2 (REGION 1C )

H3 (REGION 1R-2R-7L-70-8-9)

H4 (REGION 2C )

H5 (REGION 3 )

=-0. 34966 13E+01

= 0-5151719E-01

=-0. 2283725E+00

= 0-8237422E-02

= 0-5791229E-01

A- 98
NET MC ENERGY •-0.1147186E+03

FRAME COEFFICIENTS-REL

HTW r 0. 1107266E+03

FCMCA « 0. 2621585E-01

FCHXS • 0.3479850E-01

1

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

1

FLANK LOSS POLYSTYRENE B-B 9-01-82 TEMPB-STATIC CASE

NORMAL TERMINATION

01/16/85 13:53 : 17 USER CPU TIME 7:01.843 99.5 * SVC CPU TIME 2.082 .4 % ROLL CPU TIME 0.000 0.0 i PROCESSOR TIME 7:03.925 100.0 % 53.7 % WAIT TIME 6:04.075 46.2 t ROLL TIME 0.000 0.0 t ELAPSED TIME 13:08 100.0 % ROLLS 0 I/O 677

A-99
?

EXAMPLE PROBLEM 2: FLANKING LOSS - DYNAMIC CASE

FLANK LOSS POLYSTYRENE 09-01-82 TEMPB DYNAMIC CASE

2 -1 0 4 1 1 1 8 1 0

0 0 0 0 0 0 0 0 0 -1

-1 -1 -1 -1 -1

4 3 4 5 4 3 4

0.254 0.273 0.400 0.502 0.578 0.610 0.864

5 3 5

0.456 0.457 .957

2.0 0.0 7.45 7.45 0.0 2.33 2.33

25.0 100. 65.0 65.0 100. 38.0 38.0

0.022 0.430 0.180 0. 180 0.0360

32.00 1920. 1100. 1100. 15.10

921.

1

963.0 2000. 2000. 1214.

1 38.0 2 38.2

0.0 86400. 900. 1.0

25.0

1 1 0 6 1 20. 80.

TEST PLOT FOR CONTOURS - RJO 09/01/82

X

Y

- AXIS

- AXIS

20

1 2

17 18

NODES

TEMPERATURE (C)

3

19

4

20

10 11 12 13 14

TEST PLOT - TIME HISTORIES 08/27/82 BY RJO

X - AXIS

Y - AXIS

TIME (SEC)

TEMPERATURE (C)

5

1 2 3 4 5

0.0

100 .

15 16

A-100
i 1 4A (REV. 9-7*1

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO. a. Gov't *>'4 I. *«»»*’* Accw&lon tToT m

:

4. TITLE AND SUBTITLE

Thermal Flanking Loss Calculations for the National Bureau

of Standards Calibrated Hot Box

5. Publication Date

•. Performing Organisation Code

7. AUTHOR(S)

Ronald J. Onega and Patrick J. Burns

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, DC 20234

10. Project/Task/Work Unit No.

11. Contract/Grant No.

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

Department of Energy

National Bureau of Standards

13. Type of Report & Period Covered

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

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

16. ABSTRACT (A 200-word or leas factual aummary of moat aighificant Information. If document includea a aignlficant bibliography or

literature survey, mention it here.) A two-dimensional, finite-difference model was developed to calculate the flanking

loss for the NBS Calibrated Hot Box. A new definition of flanking loss is presented, along with the thermal theory, a description of the computer code and some

results. This model applies to both steady-state and dynamic boundary conditions.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of the first key word unless a proper name

separated by semicolons)

calibrated hot box; dynamic; finite difference; flanking loss

18. AVAILABILITY [^Unlimited

I | For Official Distribution. Do Not Release to NTIS

| | Order From Sup. of Doc., U.S. Government Printing Office, Washington, DC

20402, SD Stock No. SN003-003-

| | Order From National Technical Information Service (NTIS), Springfield, VA. 22161

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF

PRINTED PAGES

22. Price

USCOMM-DC