NATIONAL INSTITUTE OF STANDARDS k

TECHNOLOGY

Research Information Center

Gaithersburg, MD 20890

A111D2 =102173

NISTIR 88-3868

A Discrete Thermal

Analysis Method (DTAM)

for Building Energy

Simulation with DTAM1

Users Manual

James Axley

Department of Architecture

Cornell University

Ithaca, New York 14853

October 1988

Prepared for National Institute of Standards and Technology

U.S. Department of Energy

-QO

100

4)56

# 88-3868

1988

NISTIR 88-3868

too

, U- (a ho. % 8-?g(c

ry Lst a A Discrete Thermal

Analysis Method (DTAM)

for Building Energy

Simulation with GTAM1

Users Manual

James Axley

Department of Architecture

Cornell University

Ithaca, New York 14853

October 1988

National Bureau of Standards became the National Institute of Standards and Technology on August 23, 1988, when the Omnibus Trade and

Competitiveness Act was signed. NIST retains

all NBS functions. Its new programs will encourage

improved use of technology by U.S. industry.

Prepared for National Institute of Standards and Technology

U.S. Department of Energy

U.S. DEPARTMENT OF COMMERCE

C. William Verity, Secretary

NATIONAL INSTITUTE OF STANDARDS

AND TECHNOLOGY

Ernest Ambler, Director

ABSTRACT

This document includes a report that describes the theoretical basis of the program

DTAM1 and a users manual for the program.

DTAM1 is a general purpose building energy simulation program that was developed to demonstrate an approach to building energy simulation based upon discrete analysis

techniques, including, but not limited to, the Finite Element Method, used in other fields of physical simulation. It is the product of a first phase of development of Discrete Thermal

Element Analysis Techniques for Building Energy Simulation that are expected to provide a means to unify existing building energy simulation theory.

DTAM1 provides a library of discrete thermal elements, that may be assembled to model

thermal systems idealized to have constant material and heat transfer properties (i.e., linear

idealizations), including:

- ID two-node thermal resistance elements - single-node lumped capacitance elements - two-node fluid flow loop element - ID two-to-four node isoparametric conduction Finite Elements - 2D four-node isoparametric conduction Finite Elements (planar and axisymmetric)

Equations defining a variable node mean radiant temperature element are also presented in the report.

Steady state and transient analysis capabilities are included. Temperature, heat flow rate, and convective boundary conditions may be modeled and system temperature variables may

be constrained to be equal so that mixed assemblages of ID and 2D elements may be

employed.

KEY WORDS: building energy simulation, building dynamics, computer simulation

techniques, discrete analysis techniques, discrete thermal elements, dynamic simulation, finite element analysis, DTAM1

CONTENTS

ABSTRACT 1 ACKNOWLEDGEMENTS 5

R E P ORT

1. Introduction 6

2. Theoretical Basis 6

2.1 Thermal Model 6

2.2 System Degrees of Freedom 7

2.3 General Form of Element Equations 7

2.4 System Equations 9

2.5 Boundary Conditions 11

2.6 Solution of System Equations 12

2.6.1 Steady Excitation

2.6.2 Harmonic Excitation

2.6.3 General Excitation

2.7 Specific Element Equations 14

2.7.1 Resistance Elements

2.7.2 Lumped Capacitance Elements

2.7.3 Flow Loop Elements

2.7.4 Transient Conduction Finite Elements

2.7.5 Variable-Node, Mean Radiant Temperature Element

3. Summary and Conclusions 2 4

4. References 2 6

DTAM1 USERS MANUAL

I. General Instructions 2 8

II. Preparation of Input File 3 0

A. Problem Initiation . . 3 2

B. Constraint Data Group 3 2

C. Coordinate Data Group 3 3

D. Element Data Groups 3 3

3
1) Resistance Elements

2) Lumped Capacitance Elements

3 ) Fluid Flow Loop Element

4) Isoparametric ID Two To Four-node Conduction Finite

Elements

5) Isoparametric 2D Four-node Conduction Finite Elements

E. Boundary Conditions Data Giqud 3 7

F. Convection Boundary Conditions Data Group 3 7

G. Solution Control 38

1) Steady State Excitation/Response

2) Harmonic Excitation/Response

3)

General Excitation/Response - Predictor-Corrector

Integration

H. Output Report Control 41

Appendix A - Free-Field Input Conventions of the CAL-SAP Software

Development System

Appendix B - Example Problems

B1 Response of a Semi-Infinite Slab to a Step Change of Heat Flux

B2 Residential Building Example

B3 Analysis of the Huron Building Thermal Bridge

Appendix C - ElflfliaiD—Simglm

Appendix D - FORTRAN 77 Source Code

4
ACKNOWLEDGEMENTS

Most of the development of the program DTAM1 was completed during the summer of 1985 at the National Bureau of Standards, Gaithersburg, Maryland. During that time the

author was employed as a visiting researcher with the Thermal Analysis Group of the Building

Physics Division. Dr. Kusuda, the Chief of the Building Physics Division, generously offered

support for the development of the program and provided direction in a most adept and

timely manner. Dr. Richard Grot and George Walton of the Thermal Analysis Group critically

reviewed the progress of development of the program on a day-to-day basis. They provided

invaluable insight that has led to a generalization and extension of the methods presented in this document currently under development. Daniel Clark of the Building Equipment Division

introduced the author to the HVACSIM+ program and the theoretical basis of the HVAC

component models it contains. This effort was also supported by the U.S. Department of Energy through an Interagency Agreement with Lawrence Berkeley Laboratory.

5
REPORT

1. Introduction

DTAM1 is a general purpose building energy simulation program that was developed to demonstrate an approach to building energy simulation based upon discrete analysis

techniques, including, but not limited to, the Finite Element Method, used in other fields of physical system simulation. It is the product of a first phase of development of Discrete

Thermal Element Analysis Techniques for Building Energy Simulation that are expected to provide a means to unify existing building energy simulation theory. Portions of the

program are based on the Finite Element conduction analysis program, *HEAT* .written by

Robert Taylor, Department of Civil Engineering, U. C. , Berkeley [1] . DTAM1 was

developed using the CAL-SAP library of subroutines developed by Ed Wilson, Department

of Civil Engineering, U.C., Berkeley [2] .

DTAM1 provides a library of discrete thermal elements, that may be assembled to model

thermal systems idealized to have constant material and heat transfer properties (i.e., linear

idealizations), including:

- ID two-node thermal resistance elements - single-node lumped capacitance elements - two-node fluid flow loop elements - ID two-to-four node isoparametric conduction Finite Elements - 2D four-node isoparametric conduction Finite Elements (planar and axisymmetric)

Equations defining a variable node mean radiant temperature element are also presented in the report.

Static and dynamic analysis capabilities are included. Temperature, heat flow rate, and

convective boundary conditions may be modeled and system temperature variables may be

constrained to be equal so that mixed assemblages of ID and 2D elements may be

employed.

2. Theoretical Basis

2.1 Thermal Model

A building thermal system may be modeled as an assemblage of discrete thermal elements

connected at discrete system nodes. Typically, these nodes will be associated with discrete

points in the physical system, although this is not, in general, necessary. It will be shown

that system equations, governing the behavior of the building thermal system, may be

assembled from element equations that govern the behavior of individual elements.

A hypothetical example is illustrated below. Here, a simple house is idealized to model

dynamic thermal exchanges between the house and its environment using resistance

6
elements, lumped capacitance elements, and ID and 2D conduction finite elements.

Idealization of a thermal system as an assemblage of discrete thermal elements is facilitated by the use of diagrams such as the one illustrated below. The process of mathematical

assembly of discrete element equations is analogous to the graphical assembly of element

icons used in such diagrams and provides the analyst with an intuitive aide to the process of thermal idealization.

System Node

Lumped

Capacitance

ID 2-Node _

Finite Element

ID 2-Node Finite Element

MRT Sub-Network

Simple Resistance Elems.

q j

Temp. DOF

Constrained to be Equal

2D 4-Node Finite Element

Fia. 1 Hypothetical Example of a Thermal Model for a Residence

2.2 System Degrees of Freedom

With each system node we associate a temperature degree of freedom (DOF), Tj , and a

direct heat flow rate DOF, Q j (i.e., heat flow rate directly from a source into the system

node). These variables define the state of the system, at any time, and will be referred to as

the system degrees of freedom . The system DOFs may be represented as vectors for symbolic mathematical manipulation:

{T)-{Ti,T2,T3 Tn } T

; the system temperature vector (2.1)

{Q} = { Qi , 02 , 03 Qn } T

; the system direct heat flow rate vector (2.2)

for a system with n system nodes. The system temperature vector will be treated as the

principal dependent variable.

2.3 General Form of Element Equations

One or more nodes may be associated with each element. With each element node, i, we

associate an element temperature DOF, Te

j, and net heat flow rate DOF, q

e

net .j (i.e., heat

flow rate from the element node into the element):

7
These element DOFs may also be represented as vectors for symbolic mathematical

manipulation:

{T6

} = { T®i , T®2 > T®3 , , T®m } T

; the element temperature vector (2.3)

{q

e

n et} = f q

enet-1

’ ^061-2 ^net-m

}T

: the element net heat flow rate vector ( 2 - 4 )

for an element with m nodes.

In general, we attempt to describe the overall energy balance of an element by equations

of the form:

{q

e

net } - [ke

](Te

} + [c e]{dTe

/dt} - {q

e

} (2.5)

where:

[k e

] = matrix of element conductance transfer coefficients = the element conductance matrix = [k e

(t, T

e

)] , in general

[c e

] = matrix of element thermal capacity coefficients; the element capacitance

matrix = [c®(t,T®)] , in general

{q®} = vector of element-generated heat flow rates; the element internal generation

rate vector

For elements employed in the program DTAM1, the element conductance matrices, [k e

],

will be positive semidefinite and the element capacitance matrices, [c®], will be positive

definite.

Physically, the terms on the right hand side of equation (2.5) will often represent heat

flow rate into the element, at each node, by conduction, convection and/or radiation, heat

flow into storage, and heat flow generated internally within the element, respectively.

8
Other forms of element equations could have been considered, but equations of this form

embrace a large variety of thermal models and lead to systems of equations that have been

extensively studied and for which many numerical solution strategies have been published.

(Other forms based upon conduction transfer functions, transmission matrix formulations,

and equations associated with fluid flow heat transfer phenomena have been formulated and

are presently under study.)

2.4 System Equations

Requiring conservation of energy at each system node we demand:

{ (direct heat flow rate) = X (net heat f |ow rate int0 element ) Jsystem

connected nocie elements

( 2 . 6 )

or, for an arbitrary system node, n, with connected elements "a", "b", "c", ...

Qn - 1

a

n .,-k + 1

b

net-l +

^net-i + -- < 2' 7>

where k, I, i, ... are the element nodes of elements a, b, c respectively, that are associated with the system node n.

Fig. 3 Conservation of Energy at System Node

If individual element equations, of the form of equation (2.5), are substituted directly into the equilibrium relation of equation (2.7) a system of equations that relate system direct heat

flow rate DOFs, (Q) , to the collection of element temperature DOFs, {T3}, {Tb}, {TC}, ... would result. As each element temperature DOF is associated with a specific system

temperature DOF it is possible to reform these equations into the so-called system

equations , that have the form:

[ K ]{ T } + [ C ]{dT /dt) = { E } (2.8)

9
where;

[ K ] = matrix of system conductance transfer coefficients; the system conductance

matrix

[ C ] = matrix of system thermal capacity coefficients; the system capacitance matrix

{ E } = vector of source/sink flow rates; the system excitation vector

For elements employed in the program DTAM1, the system conductance matrix, [K], will be

positive semidefinite and the system capacitance matrix, [C], will be positive definite.

The process of forming the system equations from the element equations is known as the

element assembly process . It may be represented formally by first establishing the one-to- one correspondence between each element’s DOFs and the system DOFs through simple

Boolean Transformations of the form:

(T*) = [ B

e

] {T) (2.9a)

or, equivalently, as the correspondence is one-to-one:

{T} = [ B®]T {T«} (2.9b)

where:

[ B

e

] = the transformation matrix for element "e" ; a Boolean transformation matrix

consisting of zeros and ones as the element DOFs are either equal to a system

DOF (i.e., 1) or not (i.e., 0).

The same correspondence exists between each element’s net heat flow rate DOFs and the

contribution of that element to the equilibrium at each of the system nodes , represented,

here, by {Qe

} , as:

(Qe

) = [ B

e]T {Qe

net} <2.10)

Rewriting equation (2.7) in terms of the element contributions, (Qe

) , for all system nodes:

{0} = {Qa

} + {Qb} + {Qc

} + ... (2.11a)

or

{Q} = I{Q6

} (2.11b)

a, b, c, ...

and substituting the transformation expressions (2.9a), (2.10), and (2.5) we obtain

equation (2.8), above, with:

10
[ K ] = [B ( 2 . 12) a]T[ka

][B a

] + [B b]T[kb

][B b

] + [Bc]T[kc

][B c

] + ...

[ C ] = [B a]T[ca

][Ba

] + [B b]T[cb

][B b

] + [B^T^KB0

] + ... (2.13)

{ E } . {Q} + { [Ba]T{qa

} + [B^qb) + [B^qCj + ... (2 . 14)

The formal representation of the assembly process, equations (2.12) , (2.13) , and (2.14),

is useful for theoretical consideration. Practically, however, the system matrices may be

assembled more directly and efficiently using the so-called "LM Algorithm" [3] . This

algorithm is used to assemble element equations to form the system equations in the

program DTAM1.

2.5 Boundary Conditions

Some nodal temperatures, {T

p

}, may be prescribed (e.g., outside air temperature), while

the rest, (T

{ ) may be considered variable or free. At those nodes where nodal temperatures

are variable, direct heat flow rate (or in the case of some of the elements, convective flow

rates) may be prescribed. Equation (2.8) , then , may be partitioned as:

i Kfp'

(

)

i o o CL 1

Q. *1 Kpp l Tp| 1 o T>

CL CL 1 O

dll)

dt /

(Note: for lumped mass idealizations [Cfp] = [Cpf] = [0] .)

(2.15)

The first equation of equations (2.15) provides the governing equation of the free

response of the system:

[KffHTf} + [C„]{dTf/dt} = (E f ) - [K{p](Tp

} - [Cfp]{dT

p

/dt} (2.16)

The right hand side of this equation defines an effective excitation, that includes the effect

of prescribed temperatures on system response.

Equation (2.16) may be solved by a variety of methods and the solution for (Tf ) may then

be substituted into the second equation above:

{Ep} = [K

p

f]{Tf ) + [KppHTp} + [CpfHdTf/dt} + [CppKdTp/dt) (2.17)

to determine the excitation quantities, {Ep

}, needed to maintain the prescribed temperatures,

{Tp

}

.

A prescribed nodal temperature, Tj , may also be imposed numerically by scaling the

11
diagonal terms corresponding to the prescribed temperature DOF of either the system

conductance matrix - term Kjj - and/or the system capacitance matrix - term Cjj - by a large

value and setting the corresponding excitation term, Ej , to the product of the prescribed

temperature and the appropriate combination of Kjj and Cjj (depending upon the solution

type employed, steady, harmonic, or predictor-corrector integration) . This method is simple,

effective, and easily implemented. It is, therefore, employed in DTAM1.

2.6 Solution of System Equations

The solution of equation (2.8) for (T(t)} defines the system response to a given

excitation, (E(t)} Response analysis may be classified by the nature of the excitation, the

nature of the system (i.e., linear or nonlinear; constant property or nonconstant property),

and the type and characteristics of the analytical or numerical method used to obtain the

solution. Here we shall limit consideration to linear systems with constant properties,

subjected to either steady excitation, harmonic excitation, or any general excitation that may

be approximated by a piece-wise linear function.

2.6.1 Steady Excitation

For linear systems with constant properties driven by a steady excitation the response of the system will, eventually, come to a steady state (i.e., {dT/dt} = 0 ) given by the solution of:

[K]{T} = {E} (2.18)

2.6.2 Harmonic Excitation

For linear systems with constant properties driven by a steady harmonic excitation of the

form:

{E} = Re( {E*} e'wt ) ;

i = fT (2.1 9a)

where;

{E*} = the complex excitation vector (i.e., the excitation at each DOF j, E*j, is represented in terms of amplitude or modulus, E’j, and phase angle or argument

, <f>, as;

E*j = E'j e'4> = E’j (cos(<)>) + i sin(<j>) ) (2.19b)

to = circular frequency of excitation [=] radians/time

the response of the system will, eventually, come to a steady harmonic response:

{T} = Re( (T* } e'wt ) (2.20)

12
given by the solution of:

[K*]{T*J = {E*} (2.21a)

where:

[K‘] = [K + itoC] (2.21b)

m = complex temperature response = defined in terms of amplitude, T'j , and phase angle <|> , as:

T*i = T'j ei(

l> = Tj (cos(4>) + i sin(<t>)

)

(2.21c)

2.6.3 General Excitation

A finite difference scheme for the approximate integration of the semi-discrete equation

(2.8) may be developed by dividing time domain into discrete steps:

Wl = t n + 8t ; n = 0,1 ,2,3 ... (2.22)

t Q = initial time

where:

5t = integration time step (often constant but may be variable)

demanding the satisfaction of equation (2.8) at each of these steps:

[K]{T} n+ i + [C]{dT/dt} n+ i = (E} n+ i (2.23)

where:

OWl - (T(t n+i)}

(dT/dt} n+ i = {dT(tn+ i)/dt}

(E} n+1 - (E(t n+1 )}

and substituting into this equation, (2.23), the consistent difference approximation

represented by:

fOn+ 1 - {T} n + (1-a)5t{dT/dt} n + a8t(dT/dt} n+1 (2.24)

where:

0 < a < 1 a = 0 corresponds to the Forward Difference scheme

a = 1/2 corresponds to the Crank-Nicholson scheme

a = 2/3 corresponds to the Galerkin scheme

a = 1 corresponds to the Backward Difference scheme

a general implicit finite difference scheme is formulated:

13
[ a5t[K] + [C] ]{dT/dt} n+1 - {E} n+ i

- [K]{ {T}

n + (1-o)5t{dT/dt} n } (2.25a)

or, equivalently:

[ [K] + (1/a5t)[C] ]{T} n+ i - {E} n+1 + (1/a5t)[C]{T}

n + ( 1 -cx)5t[C]{dT/dt}

n (2 . 25b)

Computationally it is strategic to implement this general finite difference scheme, equation

(2.25), as a three step predictor-corrector algorithm:

rnn+1 = {Tj n + (i-o)5t(dTdt}

n

; predictor step (2.26a)

[ a5t[K] + [C] ]{dT/dt} n+ i = {E} n+1 - [K]{T}

n

; (i.e., equation (2.25a) ) (2.26b)

{T

}

n+i={T}n+ i +a5t {dT/dt} n+1 ; corrector step (2.26c)

It should be noted that this algorithm is self-starting. Given initial conditions, (T(t 0 )} , equation (2.23) may be solved to obtain an estimate of the initial rate of change of nodal

temperatures, {dT(t 0)/dt} , and the first predictor step, equation (2.26a), may then be

computed.

This predictor-corrector scheme has been analyzed by Taylor [1] and Huebner [5] and a more general predictor-multicorrector scheme that includes this implicit scheme has been

analyzed by Hughes [6]. For a> 1/2 this scheme leads to an unconditionally stable

(approximate) solution; for a >3/4 (approximately) leads to an unconditionally stable nonoscillatory solution; beyond this Taylor makes some recommendations regarding selection of a and step size, St, to limit error while minimizing computational effort. In the program

DTAM1 the default value of a is set to 0.75, and may be reset by the user, and an estimate

of the time step needed to limit error is reported (for given initial conditions) using a method

developed by Taylor [1].

2.7 Specific Element Equations

This section presents the element equations for:

-ID two-node thermal resistance elements, - single-node lumped capacitance elements, - two-node fluid flow loop elements,

-ID two-to-four node isoparametric conduction Finite Elements, - 2D four-node isoparametric conduction Finite Elements (planar and axisymmetric), and - variable-node mean radiant temperature element.

2.7.1 Resistance Elements

14
Resistance elements may be developed directly from fundamental heat transfer theory. A

2-node resistance element results that is defined by an element conductance matrix of the

general form:

[ k®] = (A/R)

where:

A = the area available for heat transfer

R = the element resistance

1 -1

-1 1 (2.27)

Fia. 4 2-Node Resistance Element

2.7.2 Lumped Capacitance Element

Lumped capacitance elements may be developed directly from fundamental heat transfer

theory. A 1-node capacitance element results that is defined by an element capacitance

matrix of the general form:

[ce

] = (me) [1] (2.28)

where:

m = mass of element

c = specific heat capacity

j

2 m , c

Fia. 5 Lumped Capacitance Element

2.7.3 Flow Loop Elements

Energy transferred by convective flow may be modeled with flow-related elements. Flowrelated element equations are best formulated in terms of nodal temperatures and the net

rate of enthalpy change, he net _j , (i.e., enthalpy transported from the node, i, into the

element). For the simple case of incompressible fluid flow from a single node to another,

ignoring fluid capacitance effects, the element equation takes the form:

15
(2.29)

where:

w-|_2 = the mass flow rate from node 1 to node 2

Cp = the specific heat capacity at constant pressure of the fluid

h

T

e &

net-1

w,

h

net-2

O t!

P'

'1-2

Fia. 6 2-Node Simple Flow Element

In principal, element equations relating rate of enthalpy change to nodal temperatures may

assembled, along with nonflow related elements, to form system equations governing

coupled flow and nonflow heat transfer. Flow related element equations are, however,

characteristically nonsymmetric, as in the case presented above. The program DTAM1 is organized to take advantage of the property of symmetry of nonflow elements and,

therefore, even this simple flow related element may not be added to its library of elements.

A subassemblage of two simple flow elements, creating a flow loop, however, results in a symmetric system of equations of the form:

(2.30)

where:

w = the mass flow rate from node 1 to node 2 and return, (here we have used

superscripts el and e2 indicate the two simple flow elements el and e2 of the subassembly)

C , w

Fig. 7 2-Node Flow Loop Element

The flow loop element equations are seen to be similar in form to the simple resistance

16
element equations and, therefore, flow loops may be modeled with resistance elements by

setting the quotient (A/R) equal to the product of (wCp), as is commonly done in thermal

network analysis. This approach to formulating flow loop equations, based upon a

consideration of rate of enthalpy change, leads to a entire class of new thermal elements

that is currently under investigation.

2.7.4 Transient Conduction Finite Elements

The solution of transient heat conduction problems using the finite element method has

been discussed by several authors [1, 3, 5, 7, 8, 9, & 10]. Over a given region, Q , bounded

by a surface, r , of a solid continuum, the spacial and temporal variation of temperature,

T(x,y,z,t) , due to conduction is governed by the so-called transient heat conduction

equation:

(k v—

ax

x dx dy (2.31a)

with temperature boundary conditions specified on part of the boundary, r-| :

T = Tp (x, y, z, t=0) ; on r-j , t > 0 (2.31b)

direct heat flow and/or convective heat flow specified on part of the boundary, r2 :

k xf^ n x + k

V T~ n v + k z n z = q(x,y,z,t) + h(x,y,z,t) (T^-T) ; onr2

,t > 0

dx y dy y *dz (2.31c)

and initial conditions:

T = T0(x, y, z, t=0) ; in Q

where;

(2.31d)

kx » k

y

• kz P

C

q*

TP

nx , n

y

, nz q

h

Too

= the principal conductivities for a thermally anisotropic continuum = the density = the specific heat capacity = the internal heat generation rate per unit volume

= the prescribed temperature boundary condition = the direction cosines of a normal to the surface at (x, y, z) = the prescribed heat flow rate boundary condition on r2 = the convective heat transfer coefficient on r2 = the temperature of the external environment that drives the

convective heat flow

In the finite element method a solution for the temperature field, T(x, y, z, t), is approximated by assumed functions that are defined over a small, but finite, part of the

region (i.e., the finite element). These functions have the general form:

17
T(x, y, z, t) = {N(x, y, z)}{Te

(t)} (2.32)

where the spacial {unctions, {N(x, y, z)}, - the so-call shape functions - are selected to satisfy continuity requirements across inter-element boundaries. The spacial variation of the

temperature field is thus approximated in a piecewise manner, using the shape functions, and

the temporal variation is captured by the variation of discrete values of temperature, {T®},

that correspond to the element temperature DOFs discussed above.

Using the approximation represented by equation (2.32) and the Galerkin variation of the

method of weighted residuals, the transient heat conduction equation is transformed to equations of the form of (2.8) that may be assembled, in the manner discussed above, from

element equations of the form of equation (2.5) where, now:

- the element conductance matrix coefficients for elements in Q. and on T-j are:

(2.33a)

- the element conductance matrix coefficients for elements on r2

are:

(2.33b)

- the element capacitance matrix coefficients for all elements are:

(2.33c)

the element internal generation rate vector components are:

(2.33d)

- the system direct heat flow rate vector is augmented by the components:

(2.33e)

18
where Qj® , the element contribution to the system direct heat flow rate at element node i, includes two terms - a direct surface gain term and a convective gain term, respectively.

A practically unlimited variety of specific conduction elements may be formulated using

these equations; no restrictions on the number of nodes per element have been made and

the compatibility requirements placed upon admissible shape functions is not very restrictive. Furthermore, even after decisions regarding the number of nodes per element and the

choice of shape functions has been made the analyst may choose to employ analytical or various numerical integration strategies to evaluate element matrices thereby generating

subclasses of element types.

Elements having distorted geometries and/or having nonuniform node spacings may be

conveniently formulated by mapping the global (i.e., actual or physical) coordinate geometry

(x,y) to an undistorted local geometry (r,s) and evaluating the integral expressions, above,

in the undistorted geometry space using relatively simple shape functions defined in the

undistorted geometry space. If the mapping from the global geometry to the local geometry

is defined using the same shape functions used to approximate the temperature field (i.e., in equation (2.32)) then the element formulation is said to be isoparametric. (An introduction

to isoparametric element formulation may be found in [8] pp. 193-230 and the practical

implementation of the approach is presented in [11].)

Two finite elements are provided in the program DTAM1 ; a variable 2 to 4-node

isoparametric one-dimensional conduction finite element and a 4-node isoparametric twodimensional finite element.

Two to Four-Node One-dimensional Isoparametric Elements

A variable two to four-node isoparametric conduction finite element may be formulated for one-dimensional heat transfer by using the following shape functions, {N} = {N-|, N2, N3, N4},

defined relative to a biunit element in local coordinate r as (see [3] page 199);

Include only if Include only if node

node 3 is present 3 and 4 are present

N, = ~r(1 - r) “0 - r 2

) 2 2

N2 = ^(1 - r) -l(1-r2

)

N3 = (1 - r 2

)

+ -^-(-9r 3 + r 2 + 9r -1) (2.34)

1

6

* —(9r 3 + r 2 - 9r - 1)

1

6

(27r 3 + 7r 2 - 27r - 7)

1

6

N4

= — ( - 2 7

3 - 9r 2 + 27r + 9)

1

6

where;

r = the local coordinate defined as indicated below:

19
Global Coordinates Local Coordinates

3

X)

r = - 1 + 1

r = - 1 0 + 1

r = - 1 -v3

+1 /3 + 1

2-

node

3-

node

4-

node

Fig. 8 Two to Fou r-Node Element Coordinate Systems

In the isoparametric formulation the element global coordinate, x, is related to the element

local coordinate, r, through a transformation defined by the shape functions given above:

x = {N}{x} ( 2 . 35 )

where:

{N} = {Ni , N2} for two-node elements

= {Ni , N2 , N3} for three-node elements = {Ni , N2 , N3 , N4} for four-node elements

{x} = {xi , X2} for two-node elements = {xi , X2 , X3} for three-node elements

= {xi , X2 , X3 , X4} for four-node elements

The integral equations ( 2 . 33) may then be rewritten in terms of the local coordinates as:

r

x2 dN, dNj

I dx dx

Cy = A pc /

X2 r +1 NjNjdx = Ape

J NjNj

A 1 "'-I

J dr

( 2 . 36a)

( 2 . 36 b)

( 2 . 36c)

where:

A = the area available for heat transfer

k = kx ; the element conductivity J = dx/dr = d/dr{N}{x} ; the Jacobian of the coordinate transformation.

20
It is convenient to evaluate these integrals numerically. Gauss quadrature may be used to exactly evaluate these integrals, but the resulting equations will tend to be "stiff" and, as a

result, will tend to demonstrate artificial oscillatory behavior [12,13]. Hughes has shown

that this problem can be mitigated through approximate evaluation of the capacitance matrix

using appropriate quadrature rules [6]. A diagonal or lumped capacitance matrix, which tends

to inhibit artificial oscillatory behavior at the expense of accuracy, will result if quadrature

points for the evaluation of the capacitance matrix, (2.36b), are selected to be equal to the

element nodal coordinates in the local coordinate system.

For the two-node element the element matrices include:

- the element conductance transfer matrix:

[k 6

) = (Ak/L)

1 -1

-1 1 12.37)

- the element thermal capacitance matrix:

[c®| = (ApcL) ( 1/2 r) r (1/2 - r) (2.38)

where:

r = 0 corresponds to a lumped capacitance model that results from using the

quadrature points (-1 , +1) = 1/12 corresponds to a higher order capacitance model that results from

using the quadrature points (-V 2/3 , +Y2/3 ) = 1/6 corresponds to a consistent capacitance model that results from using

Gauss quadrature (i.e., quadrature points -1/f3 , +1/VT ) to exactly evaluate

equation (2.36b)

and

- the element internally generated heat flow rate vector (e.g., a radiant floor) is:

(q

8

) = q‘L/2^

j (2.39)

It may be shown that the higher order capacitance matrix will minimize error in the

evaluation of temperature [6].

Ti A , L , k j,

qf P > c , q* qf

Fia. 9 One-Dimensional Two-Node Conduction Finite Element

When used with the implicit numerical integration scheme provided in the program

21
DTAM1, the higher order capacitance model, r = 1/12, results in the greatest accuracy for a given time step. The choice of r=0 results in equations identical to those associated with the

finite difference method solutions.

Four-Node Two-Dimensional Isoparametric Element

Element matrices for a four-node two-dimensional isoparametric element may be evaluated

in a similar manner as that presented above for the one-dimensional isoparametric element.

The program DTAM1 provides planar and axisymmetric elements that are evaluated using

Gauss integration of integral expressions corresponding to equations (2.33). As Gauss

integration results in the exact evaluation of these integrals these elements will result in

"stiff" equations, which tend to demonstrate artificial oscillatory behavior in regions where

the temperature field has rapidly changing and pronounced gradients.

2.7.5 Variable Node Mean Radiant Temperature Elements

A variable node nonlinear radiative element may be derived directly from Joseph Carroll's mean radiant temperature network method [14, 15, 16] that has a single node for each

surface considered and an additional node:

"The fictitious Mean Radiant Temperature node in the center of the network [that] acts as a clearinghouse for all radiative interchanges."

The surface nodes are linked to the MRT node by simple resistance elements whose

resistances are defined, nonlinearly, by Carroll's method. The element conductance transfer

matrix for this variable node radiative element will have the form:

V

Planar Element Axisymmetric Element

Fig. 10 Four-Node Isoparametric Finite Element

22
[0 = h b

KAjF’j) (-A -| F ’ 1 ) (-A2 F

' 2 ) (“A i F ' ) (A-jF’t) 0

(~A2 F

' 2 ) 0 (A 2 F' 2 )

(-A n F'n) 0 0

("AnF’n)

0

0

(An F n ) (2.40)

where:

The element DOF are numbered i=0, 1,2,3 ... n with DOF i=0 corresponding to the MRT node.

h b = 4a T3 a = Stefan-Boltzmann Constant

T = the average absolute temperature of surface

Aj = area available for radiative exchange from surface i

F'j = radiant interchange factor between surface i and the MRT node

= 1/ ( 1/Fj + (1*Ej) / £j)

Fj = 1/( 1— A j Fj/ X(AjFj) ) ; determined iteratively

£j = emissivity of surface i

Fia. 11 Variable Node MRT Element

23
3. Summary and Conclusions

An approach to building energy simulation has been presented that is based upon discrete

analysis techniques used in other fields of physical system simulation wherein system

behavior is governed by an overriding equilibrium principle (e.g., mass, momentum, or energy

conservation) that is used to define an assembly process that allows the formation of system equations from individual element equations. The element equations are developed to model the behavior of discrete regions (e.g., finite elements) or discrete processes within

the system. In the present context the governing equilibrium principle is the conservation of energy and element equations have been presented that model conductive, convective, and

radiative energy transport and that account for thermal storage processes in building

systems.

The discrete thermal analysis theory and its practical implementation in the program

DTAM1 has been presented to provide a first demonstration of this very general approach.

It has been shown that element equations may be developed from very basic considerations

(e.g., the ID thermal resistance element and the simple capacitance element), from Finite Element approximations to governing partial differential equations (e.g., the ID and 2D

isoparametric conduction Finite Elements), and from more heuristic, ad hoc physical

arguments (i.e., the variable-node MRT element). Although the discussion was limited to element formulations with constant thermal properties the extension to nonconstant and

nonlinear thermal properties is straightforward. Consequently, a practically unlimited variety

of elements could, conceivably be developed following the examples presented in this report.

In particular the nonsymmetric flow element discussed in section 2.7.3 should lead, quite

naturally, to the development of a variety of HVAC component elements.

One may also approach the formation of system equations based upon either conduction

response function theory or transmission matrix theory from an element assembly approach

and, thus, consider the development of response function elements and/or transmission

matrix elements. In these cases, however, the system equations that result will be algebraic,

rather than differential, the extension to nonlinear situations is not so straightforward, and

mixed assemblies of elements based upon the approach presented here with those based

upon these heat transfer theories will only be possible in special cases. Nevertheless, by

taking this approach to the formulation of these approaches one is able to bring the full range of building energy simulation possibilities under one unifying theory.

Steady state, steady periodic, and transient analysis options for thermal systems with

constant thermal properties (i.e., linear thermal systems) have been discussed. Additional

solution options may also be considered for linear thermal systems and a variety of

techniques exists for the analysis of thermal system with nonconstant and/or nonlinear

thermal properties; Table 1 outlines some of these possibilities.

From a theoretical point of view the discrete thermal analysis method serves to unify

existing building energy simulation theory and place it on a rigorous base. In so doing the

building energy simulation community will be in a better position to: - compare advantages and disadvantages of existing building energy simulation

24
methods, - make use of the enormous volume of research literature presently available in the

general area of discrete analysis techniques (e.g., the Finite Element and related

numerical methods), and - move on to the more demanding task of simulating coupled thermal, flow, and HVAC

system behavior for indoor air quality analysis.

From a computational point of view the discrete thermal analysis method, being based

upon distinct formal operations (e.g., formation of element matrices, assembly of element

matrices, solution of system equations, & post evaluation of element results) using a library

of thermal elements, leads to highly modular program structures and could be the basis of a

higher level programming language for building energy simulation program development.

From a practical point of view the discrete thermal analysis method allows an analyst to make use of all the "tools" that have emerged from building energy simulation research

community, to create detailed or simple idealization - as suits the purpose - that may readily

be modified as design evolves (i.e., through the addition, modification, or deletion of elements), and, for a given idealization, to consider simple to complex models of excitation

idealization as appropriate to the stage of design or level of analytical inquiry.

Excitation Model

Table 1 Solution Options for The Governing Discrete Thermal Analysis System Equations

t K ] {T } + [ C]{dT/dt } = {E}

25
4. References

[1] Taylor, Robert L., "HEAT*, A Finite Element Computer Program for Heat-Conduction

Analysis, Report 75-1, Prepared for: Civil Engineering Laboratory, Naval Construction

Battalion Center, Port Hueneme, California, May 1975

[2] Wilson, E. L. & Hoit, M. I., "A Computer Adaptive Language for the Development of Structural Analysis Programs", Computers & Structures Vol.19, No.3, pp 321-338, 1984

[3] Bathe, K.J., Finite Element Procedures in Engineering Analysis . Prentice-Hall, Inc. Englewood Cliffs, New Jersey, 1982

[4] Parlett, B., "Roundoff Error in the Solution of Finite Element Systems," Proceeding,

U.S. -Germany Symposium, Formulation and Computational Procedure in Finite Element

Analysis, M.I.T., August, 1976

[5] Huebner, K. H., The Finite Element for Engineers. Second Edition. John Wiley & Sons,

New York, 1982

[6] Hughes, T.J., "Analysis of Some Fully-Discrete Algorithms for the One-Dimensional Heat

Equation," International Journal for Numerical Methods in Engineering, Vol. 21, 1985

[7] Zienkiewicz, O.C., The Finite Element Method. Third Edition. McGraw-Hill Book Co.,

London, 1977

[8] Zienkiewicz, O.C. & Morgan, K., Finite Elements and Approximation . John Wiley & Sons,

New York, 1983

[9] Cook, R. D., Concepts and Applications of Finite Element Analysis. Second Edition . John Wiley & Sons, New York

[10] Desai, C.S., Elementary Finite Element Method, Prentice-Hall, Englewood Cliffs, New

Jersey, 1979; pp. 108-120

[11] Akin, J.E., Application and Implementation of Finite Element Methods, Academic Press,

London, 1982

[12] Emery, A.F. & Sugihara, K. "A Comparison of Some of the Thermal Characteristics of Finite-Element and Finite Difference Calculations of Transient Problems," Numerical Heat

Transfer, Vol 2, pp.97-113, 1979

[13] Schreyer, Howard L., "Nonlinear Finite-Element Heat Conduction Analysis with Direct

Implicit Time Integration," Numerical Heat Transfer, Vol 4, pp.377-391, 1981

[14] Carroll, J.A., "An MRT Method of Computing Radiant Energy Exchange in Rooms,"

Proceedings of the 2nd Systems Simulation & Economic Analysis Conference, San

26
Diego, 1980, pp. 343-348

[15] Carroll, J.A., "A Comparison of Radiant Interchange Algorithms", Proceedings of 3rd

Systems Simulation and Economic Analysis Conference, Reno, 1981

[16] Carroll, J.A. & Clinton, J.R., "A Thermal Network of a Passive Solar House",

Proceedings of 5th National Passive Solar Conference, Amherst, 1980

[17] Hoit, M.I., "New Computer Programming Techniques for Structural Engineering," Ph.D.

Dissertation, SESM, Dept, of Civil Engineering, Univ.Ca., Berkeley, 1983

27
DTAM1 USERS MANUAL

I. General Instructions

Thermal analysis of a building system, using DTAM1, involves three basic steps:

Step 1 Idealization of the Building System and Excitation

Fi g. 1.1 Idealization of the Building System and Excitation

Idealization of the building thermal system involves definition of a coordinate system,

discretization of the system as an assemblage of appropriate thermal elements connected at system nodes, identification of boundary conditions, numbering of system nodes optimally

(i.e., to minimize the bandwidth of system equations), and estimation of the bandwidth of the system equations (i.e., the largest coupled DOF number difference plus 1). Selection of appropriate thermal elements involves some skill and considerable knowledge of the

characteristics of building thermal behavior and the numerical and behavioral characteristics

of specific thermal elements.

The excitation - specified temperatures, direct nodal heat flow rates, and external

environmental temperatures associated with convection - may be modeled to be steady, to vary harmonically, or defined in terms of arbitrary piecewise-linear time histories. For the

latter case initial conditions of nodal temperatures will have to also be specified. (For those

problems where these initial conditions are not known apriori the analyst may use a false

start-up period to attempt to determine initial conditions.)

Step 2 : Preparation of Input Command/Data File

The program DTAM1 is, in essence, a command processor. The program reads an ASCII

text file that contains commands and associated data, collected together in distinct data

groups, that define the thermal idealization. The command/data input file may be prepared

with any available ASCII text editing program and given a file name, <filename>, specified by

the user. The <filename> must, however, consist of 8 or less alphanumeric characters and

Flq. 1.2 Preparation of Input Command/Data File

28
can not include an extension (i.e., characters separated from the filename by a period, ).

Step 3 : Execution of DTAM1

ooioft

01100

00011 I <filename>.DBS

ooio^i

01100 00011 <filename >.TMP

<filename>

<filename>.OUT

Fio. 1.3 Execution of DTAM1

DTAM1 is then executed. It will request the name of the command/data input file and

proceed to form element and system arrays and compute the solution to the posed problem.

DTAM1 reads the ASCII input file and creates a single ASCII (i.e., printable) output file and

up to four binary files:

Files Read

<filename> an ASCII input file specified by the user that contains problem data,

Files Created

<filename>.OUT a printable ASCII output file that contains analysis results,

<filename>.DBS a binary file containing system arrays that may be read by CAL-80 program

segments for further processing of system arrays, <filename>.TMP a binary file containing the computed time histories of nodal temperatures,

<filename>.PRE a binary file used for temporary disk storage of system arrays used by the

predictor-corrector integration subroutines in DTAM1 , and

<filename>.ELM a binary file used for temporary disk storage of element arrays.

The results of an analysis, <filename>.OUT, may be conveniently reviewed using an ASCII

editor and, from the editor, portions or all of the results may be printed out.

29
II. Preparation of Input File

Input data is defined in terms of data groups that are delimited by a data group separator

label, that may be thought to be a command, and an end flag. With the exception of the

problem initiation data group, data groups may be ordered as the user sees fit. Command/data input files follow the free-field input conventions of the CAL-SAP software

development system (see Appendix A).

The specific data groups needed to be included in a command/data input file will depend

upon the details of the thermal idealization, the nature of the excitation, and the type of analysis to be computed. For the dynamic analysis, using the predictor-corrector integration

option, an input file will have the following general form (where nn is a data value):

Column

Problem Initiation:

1DTAM1

<problem labeling information>

NODEmn B AND= nn SAVE

Constraint Data Group: CONSTRAIN

nn,nn,nn M=nn

END

Coordinate Data Group: COORDINATES

nn X=nn Y=nn Z=nn

nn X=nn Y=nn Z=nn GEN=nn,nn,nn

END

Element Data Groups:

Resistance Elements: RESIST

nn l=nn,nn R=nn A=nn

nn l=nn,nn R=nn A=nn GEN=nn

END

Capacitance Elements:

Isoparametric ID

2-Node Elements:

Isoparametric 2D

4-Node Elements:

CAPACI

nn l=nn M=nn C=nn

nn l=nn M=nn C=nn GEN=nn

END

ISOI D

nn l=nn,nn K=nn P=nn C=nn A=nn

nn l=nn,nn K=nn P=nn C=nn A=nn GEN=nn

END

IS02D4

nn l=nn,nn,nn,nn K=nn,nn,nn C=nn P=nn

nn l=nn,nn,nn,nn K=nn,nn,nn C=nn P=nn

END

Boundary Conditions Group: BOUNDARY COND

30
nn.nn.nn BC=nn

nn,nn,nn BC=nn

END

Convection Boundary

Conditions Group:

CONVECTION

nn l=nn,nn H=nn

nn l=nn,nn H=nn GEN=nn

Solution Control

- Predict-Correct. Integ.

Solution Control Data:

END

PREDICT

TIME=nn,nn,nn ALPHA=nn NEFN = nn

nn.nn.nn T=nn Q=nn TEX=nn

nn.nn.nn T=nn Q=nn TEX=nn

END

Initial Conditions Group: INITIAL COND

nn.nn.nn T=nn

nn.nn.nn T=nn

END

Excitation Function

Data:

EXCITATION

T=nn EFN=nn,nn,nn, ... T=nn EFN=nn,nn,nn, ...

Report Control Data: END

REPORT

WRITE PINT=nn

nn.nn.nn END

Details are given on the following pages for each data group.

31
A. Problem Initiation

An analysis problem must be initiated by the following lines of data:

DTAM1

<problem labeling information> : N0DE=I1 BAND = 82 [SAVE]

where:

11 = the number of nodes

12 = the (half) bandwidth; the largest DOF number difference between coupled

nodes plus 1

If the keyword SAVE is found at the beginning of the fourth line, system arrays will be

written to the binary file <filename>.DBS that may be read by CAL-80 program segments

for further processing. CAL-80 |16] is a high-level language developed to be used as an

educational/research tool that facilitates basic matrix algebraic operations and numerical

analysis methods (e.g., equation solving, eigenanalysis, fast fourier transformation, etc.). CAL-80 may be obtained from the National Information Service for Earthquake Engineering,

NISEE/Computer Applications, Davis Hall, University of California, Berkeley, California

94720, (415) 642-5113.

b. Const r aJi nl-Jaia - Gro up

System temperature DOFs may be constrained to be equal (e.g., see Fig. 1) by the

following lines of data:

CONSTRAIN

11,12,13 M=S4

END

where:

11, 12, 13 = first node, last node, node number increment of a group of slave

nodes constrained to be equal to the same master node; 13=1 default 14 = the master node to which the other, slave , nodes are constrained

Constraining system temperature DOF provides one means of using one-dimensional

elements in conjunction with two-dimensional elements. If not done with care, constraining

may have the effect of destroying the compactness of system arrays and thus result in unnecessary computational effort in solving the system equations.

32
C. Coordinate Data Group

Coordinate data is required for some, but not all, elements. It is defined by the following

lines of data:

COORDINATES

II X=R1 Y=R2 Z=R3 GEN=I2,I3,I4

END

where;

II = the node number

R1 = the x coordinate

R2 = the y coordinate

R3 = the z coordinate

12,13,14 = first node, last node, node number increment; 14=1 default

If the generation data is given, nodal coordinates are generated at equal intervals along a

straight line between the first and last nodes, 12 & 13, incrementing node numbers by the

generation increment, 14. Nodal coordinates not given will be assumed to be zero. Coordinate data is not needed for all elements; see specific element input description below

for details.

Either rectangular and cylindrical coordinate data may be defined but right hand

coordinate systems must be used:

X X (R)

D. Element Data Groups

Element data is entered by element type. Details for each element type is given below.

Each element type follows, however, the general form:

<element name>

II 1=12,13, ... GEN=I4 <element properties> [0=C1]

END

33
where:

II = the element number

12,13,... = the element node numbers

14 = generation increment (default = 1) Cl = Y if element nodal heat flow rate results are to reported (i.e., output) = N if element nodal heat flow rate results are not to be reported; N default.

Element data must be supplied in numerical order. Omitted element data is automatically

generated by incrementing the preceding node numbers by the current generation

increment. Properties of generated elements will be set equal to that of the current element.

*** The 0=Y option is presently not implemented. ***

1) Resistance Elements

Conventional 2-node resistance element equations may be added to the system equations

with the following lines of data:

RESIST

II 1=12,13 GEN=I4 R=R1 A=R2

END

where:

II

12, 13 14

R1

R2

- the element number

= the element node numbers

= generation increment (default = 1) = the element resistance = the area available for heat transfer

2} Lumped Capacitance Elements

Conventional 1-node capacitance element equations (i.e., lumped capacitance) may be

added to the system equations with the following lines of data:

CAPAC!

SI 3=12 GEN=I3 M=R1 C=R2

END

where;

S 1 = the element number

12 = the element node number

13 = generation increment (default = 1)

34
R1 = the thermal mass R2 = the specific heat capacity of the thermal mass

3) Fluid Flow Loop Elements

Elements equations that model heat transfer via a fluid flow loop between any two nodes

may be added to the system equations with the following lines of data:

FLOWLOOP

II 1=12,13 GEN=I4 W=R1 C=R2

END

where:

II = the element number

12, 13 = the element node numbers

14 = generation increment (default = 1) R1 = the mass flow rate from first node to second & second node to first R2 = the fluid specific heat capacity (at constant pressure)

4) Isoparametric ID Two to Four-node Conduction Finite Elements

Isoparametric ID two to four-node conduction finite element equations may be added to the system equations with the following lines of data:

ISOID

II 1=12, 13, [14, 15] GEN=I6 X=R1 ,R2,[R3,R4] K=R5 P=R6 C=R7 A=R8 [Q=R9]

[R=—

]

END

where:

II = the element number

12,13 = the element end node numbers (for node number order see sketch below)

14,15 «= the intermediate nodes (for node number order see sketch below):

if 13 < 0 a two-node linear element will be generated

if 13 > 0 & 14 < 0 a three-node quadratic element will be generated

if 13 > 0 & 14 > 0 a four-node cubic element will be generated

!6 = generation increment (default = 1)

R1,R2 = the element end-node coordinates (i.e., along the length of the element)

R3,R4 = the element intermediate node coordinates

R5 = the element conductivity

R6 = the element weight density

R7 = the element specific heat capacity

R8 = the area available for heat transfer

35
R9 = the internal heat generation rate per unit length of element

R=... = GAUSS or = RIO, R11, [R12.R13]

- if GAUSS Gauss quadrature will be used to evaluate the element

capacitance matrix

- RIO, R11, [R12.R13] are quadrature points to be used to evaluate

element capacitance matrix; the number of quadrature points used will be

equal to the number of element nodes (this assures exact integration when

Gauss quadrature points are used); the following quadrature points will result in a lumped capacitance modeling:

Element

2-Node Linear -1.00

RIO

1.00

Ell R 1 2 R12

3-Node Quadratic -1.00 0.00 1.00 default

4-Node Cubic -1.00 -1/3 +1/3 1.00 default

- for 2-node linear elements the quadrature points R = - V 2/3, + Y2/3 will minimize error; these quadrature points are the default quadrature points

for the 2-node linear element.

1 2

% % Two-Node Element

1 3 2 ^ Three-Node Element

1 3 4 2

% % % % Four-Node Element

Fia. 11.2 One-Dimensional Isoparametric Element Node Numbering

NOTE:

a) One may implement, in effect, a ID finite difference model, corresponding to a

linear difference approximation over the interval of the element length, by using

lumped capacitance two-node elements.

b) When using 3-node quadratic and 4-node cubic elements lumped capacitance

modeling may provide more reliable results than Gauss quadrature capacitance

modeling (see Example B1 discussed in Appendix B).

5) Isoparametric 2D Four-node Conduction Finite Elements

Isoparametric 2D 4-node conduction finite element equations may be added to the system

equations with the following lines of data:

IS02D4

II 1=12,13,14,15 GEN=I6 K=R1,R2,R3 C=R4 P=R5 [Q=R6] [D=R7] [REG=C1]

36
END

where;

II = the element number

12,13,14,15 = the element node numbers ordered counterclockwise around the

element

R1.R2.R3 = the components of the conductivity tensor, Kxx ,Kyy,K X y; for homogeneous, isotropic materials Kxx=Kyy=k , Kxy=0; for anisotropic materials

aligned to the chosen coordinate system kxx = kx , k

yy = ky

, kxy = 0

R4 = the specific heat capacity

R5 = the weight density

R6 = the internal heat generation per unit volume: R6=0.0 default

R7 = the thickness (planar regions) or subtended angle (axisymmetric regions) of the element: R7=1 .0 default

Cl = PLN for planar regions or AXI for axisymmetric regions: PLN default.

Triangular elements are possible: they are identified by repeating any two nodal point

numbers (e.g., I,J,K,K).

E. Boundary Conditions Data Group

Boundary condition flags are established with the following lines of data:

BOUNDARY CONDITIONS

11,12,13 BC=C1

END

where:

11,12,13 = the first node, last node, node number increment of a series of nodes with

identical boundary conditions

Cl = Q for heat flow rate prescribed boundaries: T for temperature prescribed

boundaries; Q default

The heat flow rate at the temperature - but not both - may be specified at each node to establish boundary conditions of prescribed temperature or heat flow rate. These conditions are specified with as many lines as needed.

If this data group is omitted all nodes will be assumed to be heat flow rate DOF.

F. Convection Boundary Conditions Data Group

Convection boundary conditions are established by the following lines of data:

37
CONVECTION BOUNDARY CONDITIONS

II 1=12,13 GEN=I4 H=R1 [D=R2] [REG=C1]

END

where:

II12,13 14R1

R2

Cl

= the convective surface segment number

= the surface segment end nodes

= the generation increment: 14=1 default = the convective transfer coefficient = the the thickness (planar regions) or subtended angle (axisymmetric regions)

of the element: R2=1.0 default = PLN for planar regions or AXI for axisymmetric regions: PLN default.

Convective boundary conditions may be specified for surface segments defined by pairs

of nodes that have been specified to be part of the Q boundary. For these nodes, then, it

is possible to account for both general direct heat flow excitation and convective excitation

(e.g., direct solar-to-surface gain and convective air-to-surface gain). These conditions are

specified with as many lines as needed.

(Convective boundary conditions would reasonably only be associated with

subassemblages of two-dimesional conduction elements. Furthermore, one would reasonably

only specify plane convective transport for subassemblages of plane elements and

axisymmetric convective transport for subassemblages of axisymmetric elements. No

attempt is made, however, to check the consistency of the model.)

G. Solution Control

Data groups needed to control the solution procedure and define the excitation are explained below.

1) Steady State Excitation/Response

The response of the system to steady excitation may be computed by including the

following lines of data:

STEADY

11,12,13 T=R1 Q=R2 TEX=R3

END

where:

11,12,13 * the first node number, last node number, and node number Increment of a

series of nodes with identical excitation conditions

R1 = the specified temperature: R1=0.0 default

R2 « the specified heat flow rate: R2=0.0 default

R3 = the specified convective fluid temperature: R3=0.0 default.

38
2) Harmonic Excitation/Response

The response of the system to steady harmonic excitation may be computed by including

the following lines of data:

HARMONIC

F-R1

11,12,13 T=R2,R3 Q=R4,R5 TEX=R6,R7

END

where:

R1 = the frequency of the harmonic excitation [=] cycles/time

11,12,13

= the first node number, last node number, and node number increment of a

series of nodes with identical excitation conditions

Harmonic Temperature Excitation:

R2 = Amplitude of Excitation: R2 = 0.0 default

R3 = Time Lag: R3 = 0.0 default Harmonic Direct (External) Flow Rate Excitation:

R4 = Amplitude of Excitation: R4 = 0.0 default

R5 = Time Lag: R5 = 0.0 default

Harmonic External (Convective Fluid) Temperature Excitation:

R6 = Amplitude of Excitation: R6 = 0.0 default

R7 = Time Lag: R7 = 0.0 default

note; Frequency = (circular frequency)/27t

Lag is defined in terms of the phase angle or argument of the complex

representations of both excitation components and response as: Lag = (-argument)^ 2n (frequency) )

(See discussion in section 2.6.2 of the REPORT, above, for clarification.)

3) General Excitation/Response • Predictor-Corrector Integration

The response of the system, including transients, to a general dynamic excitation that is defined in terms of piecewise linear functions may be computed by including the following

lines of data:

Solution Control Data :

PREDICTOR-CORRECTOR

TIME=R1,R2,R3 ALPHA=R4 NEFNsll

12,13,14

T=I5 Q=I6 TEX=I7

39
END

where:

R1.R2.R3

R4

II

12,13,14

151617

= beginning time, ending time, time step increment for dynamic analysis = integration parameter; 0.0 < R4 < 1.0 ; R4=0.75 default; instability may

result for R4 < 0.5 . = number of excitation functions that will be defined. = the first node number, last node number, and node number increment of a

series of nodes with identical excitation conditions = the excitation function number for specified temperature = the excitation function number for specified external flow rates = the excitation function number for specified external fluid temperature for convective boundary conditions

Initial Conditions Data :

INITIAL CONDITIONS

11,12,13 TrRI

END

INITIAL CONDITIONS

or RE AD = C1

END

where:

11,12,13 = the first node, last node, node number increment of a series of nodes with

identical initial temperature conditions

R1 = the initial temperature condition: R1=0.0 default

Cl = the filename of file to read initial conditions from.

Initial temperature conditions are specified by as many lines as needed. If this data group is omitted initial conditions of temperature at all nodes will be assumed to be zero.

If the READ=C1 option is used instead, initial conditions will be read from the last record

of <filename>.TMP in the form (TIME, (T(N), N=1,NNOD) ) where NNOD is the number of nodes in the system and the value TIME is ignored (see REPORT below). This option

facilitates restarting a problem to continue computation or as a means to get an estimate of

initial conditions that are not known apriori.

Excitation Function Data:

EXCITATION

T=R1 EFN=R2,R3,R4, ... END

where:

R1 =time

R1,R2,R3, ... = excitation function values for time=R1 for excitation functions 1, 2, 3,

40
... respectively

H. Output Report Control

The following lines of data control output reporting:

REPORT

[WRITE] PINTsll

12 , 13,14

END

where:

If the keyword WRITE is included the time and computed temperature solution for all nodes will be written to a binary file <filename>.TMP as a series of records of the form

(TIME, (T(N),N=1,NNOD)), where NNOD is the number of nodes in the system

II = the temperature solution "print" output interval (the solution will be written to a printable ASCII file <filename>.OUT; used by predictor-corrector integration

procedure only - otherwise ignored

12,13,14 = the first node number, last node number, and node number increment of a

series of nodes selected for printable output.

The WRITE option is useful, when used with the READ option used to specify initial conditions (above) to restart analysis.

41
.

.

'

.
Appendix A - Free-Fleld Input Conventions of the CAL-SAP Software

Development System

A "C" in column 1 of any line will cause the line to be echoed as a comment on the console.

A backslash V at the end of information on the first line will allow the next line to be

interpreted as a continuation of the first line: therefore, a 160 character record is possible.

A colon indicates the end of information on any line. Information entered to the right of the colon is ignored by the program. Therefore, it can be used to annotate an input file.

Data may be identified by an identifier of the form "<identifier=" (e.g. NUM=). Data not

identified by an identifier must be placed first in a line. If fewer data items are supplied than

required the missing data will be assumed to be zero or blank as appropriate.

Decimal points are not needed for real data. Scientific notation formats of the form E± nn

(e.g. 5.5 E-15) may be used. Simple arithmetic expressions may be used using the standard

operators +, *, *, and / . The order of evaluation is sequential from left to right - unlike

FORTRAN.

A -
5

.

.
Appendix B - Example Problems

Bl Response of a Seml-Inflnlte Slab to a Step Change of Surface Heat Flux

To demonstrate the use and characteristics of the one- and two-dimensional isoparametric

conduction finite elements the temperature response of a semi-infinite slab, T(x,t), to a step

change of surface flux, Q(0,t) = 10.0 Btu/sec ; for t > 0 , was considered. All material

properties were assigned unit values;

Four different models of this problem were considered, each representing a 6 inch deep

solid with a perfectly insulated boundary at the 6 inch depth. These models are illustrated below in Figure Bl. It may be seen that these models are similar in that an even spacial

discretization at 0.5" intervals is used in all models. Variants of each of the three models

employing one-dimensional isoparametric element assemblages based on different capacity

modeling strategies were also investigated bringing the total number of models studied to eight;

Model A1: 12 Two-node ID isoparametric elements with "consistent" (i.e. , exact

evaluation of the element capacitance matrices using two-point Gauss quadrature),

Model A2; 12 Two-node ID isoparametric elements with "lumped" element capacitance

matrices (i.e., using two-point quadrature with quadrature points of (-1.0, +1.0)),

Model A3: 12 Two-node ID isoparametric elements with "optimal" evaluation of the

element capacitance matrices (i.e., using two-point quadrature with quadrature points

of (- V273, + V2/3)),

Model Bl: 6 Three-node ID isoparametric elements with "consistent" (i.e., exact)

1 Carslaw, H.S. & Jaeger, J.C., Conduction of Heat in Solids . 2nd Ed. Oxford University Press,

Conductivity = k =1.0 Btu/sec-in-°F

Specific Heat = c = 1.0 Btu-in/sec

2

-lb-°F Mass Density = p =1.0 sec2 -lb/in4

and initial conditions were assumed to be zero, T(x,0) = 0.0 .

The exact solution to this problem is given by1

;

where;

k = k/pc

1969.

B -
evaluation of the element capacitance matrices (i.e., using three-point Gauss

quadrature),

Model B2: 6 Three-node ID isoparametric elements with "lumped" element capacitance

matrices (i.e., using two-point quadrature with quadrature points of (-1.0, 0.0, +1.0)),

Model Cl: 4 Four-node ID isoparametric elements with "consistent" (i.e., exact)

evaluation of the element capacitance matrices (i.e., using four-point Gauss

quadrature),

Model C2: 4 Four-node ID isoparametric elements with "lumped" element capacitance

matrices (i.e., using two-point quadrature with quadrature points of (-1.0, -1/3, +1/3, + 1 - 0 )),

Model D: 12 Four-node 2D isoparametric elements with "consistent" (i.e., exact

evaluation of the element capacitance matrices using four-point Gauss quadrature).

B - 2
Temperature distributions were computed for t=0.02 sec and t=0.10 sec. The results of Model A1 and Model D , the consistent ID and 2D isoparametric elements repectively, were

practically identical in this case of 1 D heat transfer. This is to be expected since identical

shape functions and numerical integration schemes were used for these cases. Results of analysis for the seven one-dimensional element assemblages are plotted below.

X (in)

Fia. B2 Temperature Profiles After a Long Time Interval

The results of all modelings are very close after long time intervals, converging to the exact

solution. It may be seen, however, that the higher order elements (i.e., the four-node and

three-node elements) provide marginally better approximations of the temperature profiles.

Although it appears that the use of "consistent" capacitance modeling for these higher order

elements provides slightly better approximations of the temperature profiles than the

"lumped" approximations the difference is slight and the "lumped" capacitance modeling

results in a diagonal system capacitance matrix, a fact that can be used to advantage in transient solution algorithms.

For short time intervals the two-node element assemblages, Models A1-A3, tend to predict

results that overshoot the exact solution, as may be seen in Fig. B3. Comparing the results of the three capacitance modeling alternatives reveals that the "consistent" capacitance

modeling provides the best estimation of surface temperature but overshoots at the second

node the most, the "optimal" capacitance modeling provides minimal overall error in the

estimation of temperature but reveals a slight overshoot at the second node, and the

"lumped" capacitance modeling suffers no overshoot at the second node but proves least

accurate in modeling the surface temperature. The lumped capacitance modeling results in a model that is mathematically equivalent to a finite difference model and as such provides a

useful comparison between finite-element and finite-difference analysis of this particular

problem.

B - 3
Fia. B3 Temperature Profiles After a Short Time Interval for Models A1-A3

Fta. B4 Temperature Profiles After a Short Time Interval for Model A and Models B & C

Finally, Fig. B4 provides a comparison of the "optimal" two-node element assemblage with

the three-node assemblages, Models B1 & B2, and the four-node assemblages, Models Cl &

C2. The higher order elements, again, provide better estimates of surface temperature than

the lower order elements with the "consistent" higher order elements out-performing the

"lumped" higher order elements. The "lumped" higher order elements. Models B2 and C2, on

the other hand, tend to exhibit less overshoot and again lead to diagonal system

capacitance matrices that can result in substantial savings in memory and solution time.

The command/data input file for Model B, using Gauss quadrature for capacitance modeling

is listed below. The output results file for this problem is also listed.

B - 4
Model B-Gauss Command/Data Input File

DTAMl

Semi Infinite Slab: MODEL B-Gauss

NODE=l 3 BAND=3

ISOID

1 1=1, 3, 2 X=0 . 0, 1 . 0, 0 . 5 K=1 . 0 P=1 . 0 C=1 . 0 A=1.0 R=GAUSS

6 1=11,13,12 GEN=2 X=0 . 0, 1.0, 0.5 K=1 . 0 P=1 . 0 C=1 .0 A=1.0 R=GAUSS:

END:

PREDICT:

TIME=0, 0.1, 0.005 ALPHA=0 . 75 NEFN=1

1 Q=1

END:

EXCITATION:

T=0 EFN=10 . 0

T=1 . 0 EFN=10 . 0

END:

REPORT:

PINT=4

lr 7,1

END:

Mpdel...B -.G3uss Output Results Fil e

| DTAMl: Building Energy Simulation for Linear Systems I Ver-6-86--

Jim Axley - Cornell & NBS

MTOT : 50000 ==== PROBLEM LABEL

SEMI INFINITE SLAB

===== PROBLEM CONTROL VARIABLES

Number of nodes 13

Maximum probable bandwidth ... 3

==== EQUALITY CONSTRAINT CONDITIONS

-- NOTE: Constraint condition data not found.

No DOFs will be constrained to be equal.

==== NODAL COORDINATES

-- NOTE: Coordinate data not found.

==== element data

== ISOPARAMETRIC ID 2 TO 4-NODE CONDUCTION ELEMENTS

Kx = Conductivity P = Density C = Capacity

A = Area Q = Int. Heat

B - 5
Elem Node Coords Kx P C A Q' Quad. Points

1 1 .00 1.0 1.0 1.0 1.0 .00 GAUSS

3 1.0

2 .50

2 3 .00 1.0 1.0 1.0 1.0 .00 GAUSS

5 1.0

A .50

3 5 .00 1.0 1.0 1.0 1.0 .00 GAUSS

7 1.0 :> 6 .50

A 7 .00 1.0 1 .

0

1.0 1.0 .00 GAUSS

9 1.0

8 .50

5 9 .00 1.0 1.0 1.0 1.0 .00 GAUSS

11 1.0

10 .50

6 11 .00 1.0 1.0 1.0 1.0 .00 GAUSS

13 1.0

12 .50 — NOTE : Needed bandwidth 3 Requested bandwidth .

.

3 — NOTE : Time to form equations

4 seconds

.

==== BOUNDARY . CONDITIONS AND EQUATION NUMBERS — NOTE : Boundary condition data not found.

All nodes assumed to be flux-prescribed nodes

Equation numbers = node numbers.

===== CONVECTIVE BOUNDARY CONDITIONS

-— NOTE: Convection boundary data not found.

==== SOLUTION

== DYNAMIC SOLUTION: PREDICTOR-CORRECTOR INTEGRATION

-- SOLUTION CONTROL INFORMATION

Start time .000

End time .100

Time step increment 0.500E-02

Integration parameter, alpha .750

Number of excitation functions ... 1 — EXCITATION FUNCTION NUMBERS

Node Temp Ext. Flux Ext. Temp.

1 10 — INITIAL CONDITIONS

•— NOTE: Initial condition data not found.

B - 6
Initial temperatures assumed to be zero. — TIME STEP

-- NOTE: Estimated time step to limit error to approx. 5.00% is:

Specified time step is: 0.500E-02 — RESPONSE

Time: 0.200E-01

Node Temp. Node

1 1.27 2

6 -0.101E-03 7

Time: 0.400E-01

Node Temp. Node

1 2.09 2

6 0 . 142E-02 7

Temp. Node Temp.

-0 . 7 60E-01 3 .117

0.48 6E-03

Temp. Node Temp.

0 . 1 67E-01 3 .102

-0 . 863E-03

Node Temp. Node

4 -0.51 6E-02 5

Node Temp. Node

4 0 . 905E-02 5

Time: 0.600E-01

Node Temp. Node Temp. Node

1 2.68 2 .178 3 6 0.126E-02 7 -0.641E-03

Temp. Node

0 . 678E-01 4

Temp. Node

0 . 188E-01 5

Time: 0.800E-01

Node Temp. Node

1 3.16 2 6 0 . 4 12E-03 7

Temp. Node Temp. Node

. 366 3 0 . 472E-01 4 0 . 101E-03

Temp. Node

0.221E-01 5

Time: .100

Node Temp. Node

1 3.56 2

6 -0.21 IE-03 7

Temp. Node Temp. Node

.564 3 0.47 6E-01 4 0. 611E-03

Temp. Node

0.216E-01 5

-- NOTE: Solution time: 9 seconds.

. 277E-02

Temp

.

0. 932E-02

Temp

.

-0 . 131E-03

Temp.

-0.44 7E-02

Temp

.

-0. 313E-02

Temp

.

0 . 853E-03

B - 7
B2 Residential Building Example

A section of a hypothetical town house and a corresponding thermal idealization is illustrated below.

Fig B5 Hypothetical Town House and Thermal Idealization

It may be seen that the model is primarily a R/C network model with a single linear one- dimensional two-node isoparametric finite element, LI, used to model the dynamics of the mass wall.

The command/input file to compute both steady and harmonic indoor air temperature

responses for specified conditions of outside air temperature and solar gain is listed below;

DTAMl

Townhouse Model:

NODE=T 9 BAND=5:

RESIST: Consistant Units of Btu, ft, hr

1 1=1,2 A=17 R= 1/1.5 R1 Air Film at Sunspace

2 1=2,3 A=17 R=0 .

9

R2 Sunspace Glass

3 1=3, 4 A=17 R= 1/1.5 R3 Air Film at Sunspace

4 1=4,5 A=8 R= 1 / 1 .

5

R4 Air Film at Mass Wall

5 1=6,7 A=8 R=1 / 1 .

5

R5 Air Film at Mass Wall

6 1=7, 8 A=8 R-l/1.5 R6 Air Film at N. Wall

7 1=8, 9 A=8 R=ll R7 N. Wall

8 1=9, 10 A=8 R-l/1.5 R8 Air Film at N. Wall

9 1=7,11 A=24 R—l/1 .

R9 Air Film at Ceiling

10 1=11, 12 A=2 4 R=1

RIO Ceiling

11 1=12, 13 A=24 R—l/1 .

6

Rll Air Film at Ceiling

12 1=13, 14 A=13 R= 1/1.3 R12 Air Film at Roof

13 1=13,17 A=13 R=l/1 .

3

R13 Air Film at Roof

14 1=14, 15 A=13 R=1

R14 Roof

15 1=17, 18 A=13 R=1

R15 Roof

16 1=15, 16 A=13 R=l/1 .

6

R16 Air Film at Roof

17 1=18,19 A=13 R=1 / 1 .

6

R17 Air Film at Roof

END

CAPACI : Consistant Units of Btu, ft, hr

B-8
1 1=4 M=150*12/3 C=0.16 :C1 Sunspace Quickmass (4" cone, floor)

2 1=7 M=150*24/3 C=0.16 :C1 Room Quickmass (4" cone, floor)

END:

IS01D: Consistant Units of Btu, ft, hr

1 1=5,6 X=0. 0,8.0/12.0 K=1 . 0 P=150 C=0.2 A=8 :L1 8" Concrete Wall

END:

BOUNDARY

1 BC=T

10 BC=T

16 BC=T

19 BC=T

END:

:Outside air temp to be specified

:Outside air temp to be specified

:Outside air temp to be specified

:Outside air temp to be specified

STEADY:

1 T=7 .

5

10 T=7 . 16 T=7 . 19 T=7 .

5

4 Q=200*13/3.1416

END

HARMONIC:

: Out side air temp - Steady Component

: Outside air temp - Steady Component

: Outside air temp - Steady Component

: Outside air temp - Steady Component

: Sunspace Solar Gain - Steady Component

F=1 .0/24

1 T=10 , 10

10 T=10, 10

16 T=10, 10

19 T=1 0 , 1

0

4 Q=200*13*0.5,6

END

REPORT

P INT=1

If 19,

: Outside air temp - Harmonic Component

: Outside air temp - Harmonic Component

: Outside air temp - Harmonic Component

: Outside air temp - Harmonic Component

: Sunspace Solar Gain

END:

B - 9
B3 Analysis of the Huron Building Thermal Bridge

In this example, assemblages of two-dimensional isoparametric conduction elements and

simple resistance elements are used to model the details of heat transfer in an exterior

masonry wall section, at the intersection of a concrete floor slab, Figure B6, of a recently

constructed federal office building in Huron, South Dakota, studied earlier by Grot et. al

2

.

4" Face Brick

6" Lt. Wt. Cone. Block

3" Fiber Glass Insul.

5/8" Gypsum Board

Expansion Joint

Reinforced Concrete

Slab

A A A ,A A AAA

Steel Beam

Fireproofing

! V/7?

v

Air Space

Fia. B6 Huron Building Wall Section

2 Grot, R.A., Childs, K.W., Fang, J.B., & Courville, G.E., "The Measurement and Quantification of Thermal Bridges in Four Office Buildings", ASHFtAE Trans. 1985, V. 91, Pt. 1, CH-85-11 No. 4

B - 10
Field thermographic measurements of this building revealed a thermal bridge (i.e., highly

conductive path through the wall construction) due, evidently, to the penetration of the

floor slab that acted to significantly increase building heat loss. Analysis of heat transfer in the vicinty of this thermal bridge is complicated by the complex two-dimensional nature of the geometry and the air space between the steel beam and the outer wall.

Two models of the wall section were considered. In the first model, Model 1, an assemblage of two-dimensional isoparametric elements was used to model the simpler

construction from the center of the concrete slab upward, Fig. B8; heat transfer was

assumed to be symmetric about the slab centerline. In the second model, Model 2, an assemblage of two-dimensional isoparametric elements were used to model, in detail, the

geometry and discontinuous material properties of the actual construction and a subassembly

of simple resistance elements was used to model convective/radiative heat transfer in the air space adjacent to the steel beam, Fig. B9. The response of both models to a steady state temperature drop of 30.5°F, from inside-air to outside-air, and a steady excitation of a harmonic variation of outside air of a 12.5°F amplitude and 24 hour period were computed.

The response of both models was similar. Representative results, for the temperature field through the construction, are shown in Fig. B7. It is seen that the steady state response

dominates the behavior of the construction and the thermal path provided by the slab is evident. Surface fluxes, calculated from the computed results closely agreed with measured

behavior and other analytical results.

0 1 2 X 3

Fig, B7 Computed Temperatu re .fjfild far Model .1

Model 1 employed 34 elements connected at 46 nodes and Model 2 employed 105

elements connected at 139 nodes. Computation to form and assemble the element

equations for Model 2 took 73 seconds and solution for both the steady state and steady

B- 11
harmonic excitation took an additional 48 seconds using a Macintosh microcomputer with a MC 68000 microprocessor and 51 2K of central memory.

NV/VN \ ' Zero Flux Boundary

S. \ \ > S s s s > \ \ \ ' 5 / / / 5 \ \ \ ' j \ \ s ' i / / / j N N N ' J

/ / i \ \ \ ' & / * y \ \ \ ' s / s * >

\ \ \ '

fc / / / a

^ \ s ' 3 / / / 3 \ \ \ ' 3 / / / 3 \ \ \ ' 3 N S N ' 5 \ N N ' J / / / 3 S \ N ' 3 / / / 3 V/V S W5

^ \ '3

/ / / 3

• / / s

* \ ' S / / 5 ^ s '5

s \ \ 'N

/ / / 3

\ \ s '3

s \ \ ' S

/ / / 3

V s \ ' 3 / / / 3

»s \ \ '3

v n s ' 5 / / / 3 s s x '3

/ / / 3

X X N / / / 3

x x x '3

/ / / 3 x x x ' 3 x x x

x x x ' J / / / S

x x x ' J / / / ^

X x x ' S

X X X ' J / / / 3 x x x ' > / / / i w/ S

x \ N ' > / / / 3

X X X ' i

* & =

r-r;p 7-q

Mll-ssl

-A- -2-

0.09 50.0 0.26 5/8" Gyp Board

0.02 6.0 0.23 3" Fiber Glass Board

0.42 116.0 0.21 6” Reinforced Lt.Wt. Concrete Block*

0.75 130.0 0.19 4" Face Brick

0.04 23.0 0.14 Expansion Joint**

0.75 140.0 0.22 Reinforced Concrete

(Units: consistent units of Btu, ft, hr, F)

* Reinforced light weight concrete block with half of cores filled assumed to be equivalent to cement

mortar

** Expansion joint assumed to be equivalent to mineral

fiber board

1

Zero Flux Boundary

X

Fig, B8 Huron Building Wall .Seslign MsdeLL

B - 12
—P—

0.09 50.0 0.26 5/8" Gyp Board

0.02 6.0 0.23 3” Fiber Glass

0.42 116.0 0.21 6" Lt.Wt. Cone.

0.75 130.0 0.19 4" Face Brick

0.04 23.0 0.14 Expansion Joint

0.75 140.0 0.22 Reinforced Cone

26.2 489.0 0.12 Steel Beam

0.21 80.0 0.32 Fireproofing

X

Subassembly of Simple Resistance

Elements to Model Air Space

Eia B9 Huron Building Wall Section Model 2

B - 13
'

-

.
Appendix C - Program Structure

The program DTAM1 is structured in a manner similar to a command processor.

Structurally it was developed as a collection of trees of subroutines linked to the main

program through roots

.

Fia. Cl Generic Trees Structure

Each tree was developed to be largely independent, dynamically defining arrays and

seeking only the input data that it alone requires.

Root of Tree: - dynamically DEFINE's arrays used by Tree

Tree: - gets only input data needed by tree

- executes command operation

Locally shared labeled COMMON & SUBROUTINES

Fig. C2 Generic Tree of Subroutines

This independence allows incremental processing of input data and reporting of results so

that a single error in the command/data input data file will not result in program termination.

It is also hoped that the hierarchical structure will help to facilitate future program

development efforts.

The CAL-SAP development software subroutines [2] are used by all trees for dynamic array

management and input data interpretation. The CAL-SAP dynamic array management is accomplished by storing array values and a directory to those values in a single blank

common vector IA(MTOT) as indicated below in Fig. C3. One may change the array capacity

of the program by simply resetting the variable MTOT and redimensioning IA(MTOT) in the

C-1
main program declaration.

Arrays stored low

inblark common.

Directory of arrays with pointers

to arrays stored high in blank common.

t*

IA(1)

[A] {B} 1C] (D) DIR D DIR C DIR B DIR A

MTOT- 4

IA(MTOT)

Fia. C3 CAL-SAP Dynamic Array Management

Specifically, the main program DTAM1 calls six trees of subroutines;

Fia. C4 Structure of DTAM1

RID is the root to the ID Tree that establishes node-number to equation-number

correspondance, accounting for any constraint input; it creates an ID array. RCOORD is the root to the Coordinate Tree that inteprets coordinate input data and

generates nodal coordinates; it forms the XYZ array. RFORM is the root to the FORM Tree that forms element arrays and assembles them to form system arrays; it forms the system conductance matrix, K(NEQN.MBAND),

the system capacitance matrix, C(NEQN,MBAND), and the internal generation

heat flow rate vector, EO(NEQN) - where NEQN=the number of equations and

MBAND=the system (half) bandwidth.

RBOUND is the root to the BOUND Tree that establishes boundary conditions by

modifying the ID array. RCONV is the root to the CONV Tree that establishes convective boundary conditions by

modifying the system conductance matrix and stores convective coefficients in array CH(NNOD), where NNOD = the number of nodes in the system.

RLSOLV is the root to the LSOLV Tree that interprets solution specification and

excitation data and implements solution by calling STEADY, HARMON, and/or

PREDIC to affect steady state, steady harmonic, or predictor-corrector solution

options.

C - 2
APPENDIX D

FORTRAN 77 SOURCE CODE

C

The FORTRAN77 source code for DTAM1 is listed below.

PROGRAM DTAM1 C LINEAR: INTERNAL VARIABLES

ERR*. FALSE.

C —PRO: DTAM1 - A BUILDING ENERGY SIMULATION PROGRAM BASED ON A C C DISCRETE THERMAL ANALYSIS METHOD FOR LINEAR IDEALIZATIONS C—2.0 WRITE BANNER & OPEN INPUT AND OUTPUT DATA FILES

C C

c DEVELOPED BY JAMES AXLEY . CORNELL UNIVERSITY WRITE (NTM, 2200) MTOT

C NBS BUILDING PHYSICS DIV. 1985-86 CALL IFILE

C USING; CALL POPENCOUT ’) C A) CAL-SAP LIBRARY OF SUBROUTINES DEVELOPED BY ED WILSON. WRITE (NOT. 2200) MTOT

C U.C. BERKELEY

C B) MacFort ran V2 . 0 4 Microsoft Fortran V2.1 COMPILERS 2200 FORMAT (// C ANSI FORTRAN 77 STANDARD .* C • INCLUDING; • ' / C 1. USE OF SAVE /contnonname/ TO RETAIN DEFINITION OF .*1 DTAM1 : Building Energy Simulation for Linear Systems

C LOCAL, SUBROUTINE VARIABLES AND SOME NAMELIST VARIABLES . I'/ C ‘EXCEPT; . * C 1. TYPE MacFort ran EXTENSION "TYPE" USED FOR . Ver-6-86— •/ C PROMPTS TO ALLOW IN-LINE RESPONSES .55X,’Jlm Ax ley - Cornell 4 NBS*/

C*« •••*••*•••••••••••*••• . 65X, 'MTOT: ’110) c c IMPLICIT REAL* 8 (A-H.O-Z) C--3.0 WHILE ERR*. FALSE.

C C C**** CAL-SAP: DATA 4 COMMON STORAGE **•*•*•*•*••*•**••*•****•••••••** C 3.1 GET PROB LABEL 4 CONTROL INFORMATION

C C

CAL-SAP: LOGICAL UNIT NUMBERS

NTM - 9 NTR - 9 NIN « 10 NOT = 11 ND1 « 12 ND2 - 13 ND3 - 14 ND4 * 15

CHARACTER FIN*1.EXT*1 COMMON MTOT. NP. IA (50000) COMMON /DBSYS/ NUMA . NEXT. IDIR, IP (3

)

COMMON /IOLIST / NTM, NTR, NIN, NOT, ND1 , ND2 , ND3 . ND4 COMMON /PARC/ FIN (12) . EXT (8 0) CC DICTIONARY OF VARIABLES

C

C VARIABLE

C MTOT

C NP

C IA (MTOT) C NUMA

C NEXT

C IDIR

C IP (3) C NTM

C NTR

C NIN

C NOT

C ND1 thru ND4

C FIN (12) C EXT (80) c **..»**»*t..t*t*«i

CC**** LINEAR: DATA 4 COMMON STORAGE

C

DESCRIPTION

SIZE OF BLANK COMMON VECTOR IA CURRENT DATA TYPE: 1 "INTEGER; 2=REAL ; 3*LOGICAL

BLANK COMMON VECTOR

NUMBER OF ARRAYS IN BLANK COMMON DATA BASE NEXT AVAILABLE STORAGE LOCATION IN BLANK COMMON START OF DIRECTORY IN BLANK COf**ON NUMBER OF LOGICALS IN EACH DATA TYPE LOGICAL UNIT NUMBER FOR TERMINAL /SCREEN LOGICAL UNIT NUMBER FOR TERM INAL /KEY BOARD LOGICAL UNIT NUMBER FOR INPUT DATA FILE LOGICAL UNIT NUMBER FOR OUTPUT DATA FILE LOGICAL UNIT NUMBERS FOR GENERAL USE INPUT DATA FILE NAME

FILE EXTENSION NAME FOR GENERAL USE

LOGICAL*l CHARACTERERR YESNO*l . SAVE *4

C DICTIONARY OF VARIA L E S

VARIABLE DESCRIPTION

ERR DO-WHILE TERMINATOR FLAG YESNO CONTINUATION FLAG NNOD NUMBER OF NODES IN SYSTEM NEON NUMBER OF (SYSTEM) EQUATIONS

MBAN (HALF) BANDWIDTH OF SYSTEM EQUATIONS POINTERS TO BLANK COMMON LOCATIONS

C MPXYZ

C MPC

C MPK

C MPEO

VECTOR

C MPT

C MPID

ARRAY

CC

NODE

C

C MPCH

FLUX

C MPNPR

C

XYZ (NNOD , 3) : COORDINATE ARRAY

C(NEQN.MBAN) : CAPACITY MATRIX

K (NEON. MBAN) : CONDUCTANCE MATRIX EO (NNOD) : INTERNALLY GENERATED EXCITATION

T (NEQN) ID (NNOD)

CH (NNOD)

NPR (NNOD)

NODAL TEMPERATURE VECTOR

NODAL EQUATION NUMBERS B. C. FLAG

EQTN NO. « ABS (ID) WEG. ID - TEMPERATURE PRESCRIBED

POS. ID - FLUX PRESCRIBED NODE COEFS TO COMPUTE EFFECT. CONVECT.

OUTPUT PRINT CONTROL ARRAY

C —1.0 INITIALIZE INTERNAL VARIABLES

C CAL-SAP: DATABASE

10 MTOT * 50000 NUMA - 0 NEXT * 1 IDIR = MTOT

IP(1) 4 IP (2) * 8 IPO) « 1

CALL FINDN (

' DTAM1 ’ , 5 . KEY

)

IF (KEY . EQ. 1 ) THEN WRITE (NTM. 2310)

WRITE (NOT, 2310)

GO TO 999 ENDIF

2310 FORMAT ( ' *•*• ERROR: Separator "DTAMl * not. found.’)

CALL FREETY WRITE (NTM. 2312) WRITE (NOT. 2312) 2312 FORMAT (/’ **= PROBLEM LABEL’/) CALL FREE

CALL FREEPT

WRITE (NTM. 2316) WRITE (NOT, 2316) 2316 FORMAT (/’ «“= PROBLEM CONTROL VARIABLES’) NNOD=0

MBAN=0 CALL FREE

CALL FREETY

CALL FREEI ('E’ .NNOD. 1) IF (NNOD. LE. 0) THEN

WRITE (NTM, 2314)

WRITE (NOT. 2314)

ERR= . TRUE

.

ENDIF

2314 FORMAT ( ’ •*•* ERROR: Number of nodes must be greater than 0.’)

CALL FREEI (’D’ . MBAN. 1) IF (MBAN .LE . 0) MBAN - 5 + INT (SQRT (REAL (NNDO) )

)

WRITE (NOT. 2318) NNOD, MBAN

2318 FORMAT (/

.* Nurcber of nodes *,I5./

.' Maximum probable bandwidth ...*,15)

CALL FREE

CALL FREEH (' ’.SAVE, 4.1)

IF (ERR) OO TO 999 C

C 3.2 ESTABLISH EQUATION NUMBERS

C CALL RID (MPID, NNGD, NEQN, ERR)

IF (ERR) OO TO 999

CC 3.3 GET 4 GENERATE NOOAL COORDINATES

C

CALL RCOORD (MPXYZ. NNOD, ERR)

IF (ERR) OO TO 999

CC 3.4 FORM ELEMENT ARRAYS AND ADD TO SYSTEM ARRAYS

C CALL TIME (ITIME1

)

CALL RFORM (MPID, MPXYZ . MPC , MPK. MPEO, NNOD. NEQN . MBAN . ERR) CALL TIME (ITIME2)

WRITE (NTM. 2340) (ITIME2-ITIME1

)

WRITE (NOT. 2340) (ITIME2-ITIME1

)

2340 FORMAT(/' — NOTE: Time to form equations: *,I5,' seconds.*) D -1
IF (ERR) GO TO 999 C

C 3.5 GET BOUNDARY CONDITION FLAGS

C CALL RBCOND (MPID. NNOD. ERR) IF (ERR) GO TO 999 CC 3.6 MODIFY K (NEQN. MBAN) FOR CONVECTIVE TRANSFER

C CALL RCONV (MPID. MPXYZ. HPK, MPCH . NNOD, NEQN , MBAN . ERR) IF (ERR) GO TO 999 C

C 3.7 SAVE SYSTEM ARRAY DATABASE

C

IF (SAVE. EQ. ’SAVE’ ) CALL SAVE

C

C 3.8 SOLVE SYSTEM OF EQUATIONS C CALL RLSOLV (MPID, MPK. MPC, MPEO. MPCH, NNOD. NEQN, MBAN, ERR)

CC—4.0 CLOSE INPUT AND OUTPUT DATA FILES

C

999 CALL FCLOSE (NOT) CALL FCLOSE (NIN) C

C —5.0 SOLICIT ANOTHER PROBLEM

C TYPE 2500 READ (NTM, 2510) YESNO

IF ( (YESNO.EQ. *Y') .OR. (YESNO.EQ. *y’ ) ) GO TO 10 2500 FORMAT (/' •• Do you want to consider another problem? (Y/N) ’) 2510 FORMAT (LAI)

STOP END

C *M*t****t**..M«t .t«t** **t*tt*t*«**««*l.**M*t*«t*M**Mtt*.**.*

C TREE: ID - TREE OF SUBROUTINES TO ESTABLISH EQUATION NUMBERS

C RID SUBROUTINE RID (MPID. NNOD, NEQN. ERR) C—SUB: RID - ROOT TO ESTABLISH THE NODE NUMBER-EQUATION NUMBER

C CORRESPONDANCE

C PRESENTLY IMPOSES EQUALITY CONSTRAINTS AND NUMBERS EQUATIONS

C *•• IN NODE NUMBER ORDER

C *** TO BE USED SUBSEQUENTLY TO; C •** 1. OPTIMIZE EQUATION NUMBERING TO MINIMIZE BANDWIDTH

C 2. ALLOW USER REDEFINITION OF NUMBERING

C

C CAL-SAP: DATA 6 COMMON STORAGE

C COMMON MTOT.NP, IA (1) COMMON /DBSYS / NUMA . NEXT . IDIR, IP (3

)

COMMON /IOLIST/ NTM, NTR, NIN , NOT. ND1 , ND2 , ND3 . ND4 C

C COORD: DATA t COMMON STORAGE

C LOGICAL* 1 ERR

SAVE /DBSYS/. /IOLIST/

C

C —0.0 WRITE HEADER

C WRITE (NTM, 2000) WRITE (NOT. 2000) 2000 FORMAT (/' = EQUALITY CONSTRAINT CONDITIONS'/)

CC—1.0 DEFINE ID ARRAY & INITIALIZE

C CALL DELETE (' ID ') CALL DEFINI ('ID MPID . NNOD. 1

)

CALL ZEROI (IA (MPID) , NNOD. 1) C

C 2.0 FIND SEPARATOR <CGNSTR>AIN; NUMBER IN NODE ORDER IF NOT FOUND

C CALL FINDN ( 'CONSTR* , 6, KEY) IF (KEY . EQ. 1) THEN WRITE (NTM, 21 00) WRITE (NOT. 2100)

C IMPLICIT REAL* 8 (A-H.O-Z) C

C CAL-SAP: DATA 6 COMMON STORAGE

C COMMON /IOLIST/ NTM. NTR. NIN , NOT . ND1 . ND2 , ND3 . ND4 C

C ID: DATA i COMMON STORAGE

C INTEGER MSTR, MNODE. MDOF CHARACTER ENDFLAG* 3 , BC*1 DIMENSION ID (NNOD), IJK(3) LOGICAL* 1 ERR

SAVE /IOLIST/

C

C—1.0 READ CONSTRAINT DATA AND GENERATE CONSTRAINT FLAGS

C

10 CALL FREE

CALL FREETY

C CHECK FOR "END- CALL FREEH (

'

' .ENDFLAG. 3.1) IF (ENDFLAG. EQ. 'END' ) GO TO 20 C GET MASTER NODE

CALL FREEI ('M' .MSTR, 1) IF (MSTR. LE.O. OR. MSTR. GT. NNOD) THEN WRITE (NTM. 2100) MSTR

WRITE (NOT. 2100) MSTR

ERR - .TRUE. ENDIF

2100 FORMAT (’ ••*• ERROR: Master node* ,15.* Is out of range.’) C GET SLAVE NODE GENERATION DATA

CALL FREEI (' ’.IJK(1).3)

IF (I JK (2) .EQ. 0) IJK(2)-IJK(1)

IF (I JK (3) . EQ. 0) I JK (3) =1 WRITE(NOT. 2110) MSTR.IJK

2110 FORMAT (’ Master Node - ',13, • Slave Nodes = ' , 13. ' to ',13

* Step ’,13) DO 12 N-IJK(l) ,IJK(2).IJK(3)

IF (N . LE . 0 . OR. N . GT . NNOD) THEN WRITE (NTM. 2120) N WRITE (NOT, 2120) N ERR «=. TRUE. GO TO 10 ENDIF

12 ID (N ) = -MSTR GO TO 10 2120 FORMAT (' *•*• ERROR: (Generated) Node '.15.' is out of range.’)

CC--2.0 NUMBER UNCONSTRAINED DOF (ID=0) C

20 IF (ERR) RETURN

NEQN *= 0 DO 22 N*=l,NNOD

IF (ID (N) . EQ. 0) THEN NEQN = NEQN 1

ID (N) «= NEQN ENDIF

22 CONTINUE

CC--3.0 NUMBER CONSTRAINED DOF (ID=NEG) TO MASTER NODE (WODE) DOF

(MDOF) C DO 30 N=1 , NNOD

IF (ID (N) .LT. 0) THEN MNODEeABS (ID (N)

)

MDOF - ID(WODE)

IF (MDOF.LT. O) THEN

WRITE (NTM. 2300) N. MNODE

WRITE (NOT. 2300) N, MNODE

ERR « .TRUE. NEQN - NNOD

DO 20 N-MPID. MPID^NNOO

20 IA (N) - M-HPID+1

RETURN

ENDIF

2100 FORMAT ( * — NOTE: Constraint condition data not found.*/

.* No DOFs will be constrained to be etjual.') CALL FREETY

CC—3.0 CALL ID TO DO THE WORK

C CALL ID (IA (MPID) .NNOD. NEQN. ERR)

ELSE

ID (N) - »OF

ENDIF

ENDIF

30 CONTINUE

2300 FORMAT ( ' **•* ERROR: Slave node ',15.' Is constrained to slave node ',15,/

' Multilevel constraints are not permitted.')

RETURN

END

RETURN

END

C ID SUBROUTINE ID (ID. NNOO, NEQN . ERR) C-SUB : ID - READS AND GENERATES EQUALITY CONSTRAINT DATA

C AND NUMBERS DOFS VIA NODE ID ARRAY

C*********************************************************************

C TREE: COORD - TREE OF SUBROUTINES TO FORM XYZ COORDINATE ARRAY C ttt«t** *M*.****tM*t«*MMM»*ttt**t*M*t*tMt***t*Mt*.tt.ttt*.

C RCOORD SUBROUTINE RCOORD (MPXYZ .NNOD. ERR) C—SUB : RCOORD - ROOT SUBROUTINE TO COORD

C COORD - READS COORDINATE INFORMATION AND FORMS X C Y ARRAYS

C D -2
c dictionary of variables

C VARIABLE DESCRIPTION

C VARIABLES PASSED TO SUBROUTINE

C ERR DO-WHILE TERMINATOR FLAG

C NNOD NUMBER OF NODES IN SYSTEM

C MPXYZ POINTER TO XYZ (NNOD) IN BLANK COMMON

C XYZ (NNOD) COORDINATE ARRAY (ORDERED BY NODE NUMBER)

C-C- CAL-SAP: DATA 4 COMMON STORAGE

C COMMON MTOT.NP. IA (1) COMMON /DBSYS / NUMA. NEXT. IDIR. IP (3

)

COMMON /IOLIST/ NTM, NTR, NIN , NOT, ND1 . ND2 , ND3 , ND4

CC COORD: DATA 4 COMMON STORAGE

C LOGICAL*l ERR

SAVE /DBSYS/, /IOLIST/

CC--0.0 WRITE HEADER

C WRITE (NTM, 2000) WRITE (NOT, 2000) 2000 FORMAT (/' «= NODAL COORDINATES '

)

CC—1.0 FIND SEPARATOR <COOR>DINATES

C CALL FINDN ( 'COORD '

, 5, KEY) IF (KEY . EQ. 1) THEN WRITE (NTM. 2100)

WRITE (NOT. 2100)

RETURN

N2 - NN (2) N3 - NN (3) IF (N1 . EQ . N2) GO TO 300

IF (N2 . EQ . 0 ) GO TO 300

IF (N2 . EQ.N1) GO TO 300

IF (N3 . EQ. 0) N3 = 1 IF (N1 .GT.NNOD) N1 = NNOD

IF (N2.GT.NNOD) N2 » NNOD

NDIF « (N2-N1) /N3 DIF « DBLE(NDIF)

XDIF - ( XYZ (N2 , 1 ) - XYZ(Nl.l) ) /DIF YDIF « ( XYZ (N2, 2) - XYZ(N1,2) ) /DIF

ZDIF = ( XYZ (N2 , 3) - XYZ(N1,3) ) /DIF

C GENERATE ADDITIONAL COORDINATES

XX « XYZ (Nl. 1) YY - XYZ (Nl , 2) ZZ XYZ (Nl , 3) Nl - Nl N3

N2 - N2 - N3 DO 200 J-N1.N2.N3

XX - XX XDIF

XYZ ( J, 1) - XX

YY « YY YDIF

XYZ ( J, 2) « YY

ZZ - ZZ ZDIF XYZ ( J, 3) - ZZ 200 CONTINUE

300 CONTINUE

C PRINT FINAL COORDINATES

400 DO 500 I-l.NNOD

X - XYZ (1,1)

Y - XYZ (I, 2)

Z « XYZ (1.3) WRITE (NOT. 2000) I.X.Y.Z

500 CONTINUE

C

900 RETURN ENDIF CALL FREETY

CC —2 . 0 DEFINE ARRAY

C CALL DELETE (’XYZ ’) CALL DEFINE (’XYZ MPXYZ, NNOD. 3) C

C—3 . 0 INITIALIZE ARRAYS

C CALL ZEROR (IA (MPXYZ) .NNOO, 3) CC— 4.0 CALL COORD TO DO THE WORK

C CALL COORD (IA (MPXYZ) .NNOD, ERR)

RETURN

2100 FORMAT (/* — NOTE: Coordinate data not found.’)

END

CCOORD

SUBROUTINE COORD (XYZ. NNOD. ERR) C

C SUBROUTINE TO READ AND GENERATE COORDINATES

C DEVELOPED BY MARC HOIT - U.C. BERKELEY

CC IMPLICIT REAL* 8 (A-H.O-Z)

C CHARACTER ENDFLAG *3 DIMENSION XYZ (NNOD. 3) .NN (3) COMMON /IOLIST/ NTM, NTR , NIN . NOT, ND1 , ND2 . ND3 . ND4 SAVE /IOLIST/

C WRITE (NOT, 1000) C INITIALIZATION X - 0.0

Y - 0.0

Z = 0.0 NODE - 0 C READ LINE OF COORDINATE INFORMATION DO 300 I-l.NNOD CALL FREE

CALL FREETY

CALL FREEH (’ ’ , ENDFLAG, 3 , 1)

c FORMAT SPECIFICATIONS

C 1000 FORMAT (/

1 SX, ’ NODE X-COORD. Y -COORD . Z -COORD. ’

)

C2000 FORMAT (7X, 13. 9X.G10.3. 4X, G1 0 . 3. 4X. G10 . 3) C

END

TREE: FORM - TREE OF SUBROUTINES TO FORM AND ASSEMBLE ELEMENT EQTNS

RFORM SUBROUTINE RFORM (MP ID. MPXYZ . MPC, MPK, MPEO, NNOD , NEQN. MBAN. ERR) --SUB: RFORM - ROOT SUBROUTINE TO ALL ELEMENT SUBOUTINES THAT FORM 4 ASSEMBLES ELEMENT ARRAYS DICTIONARY OF VARIABLES

VARIABLE

VARIABLES PASSED

ERR NNOD

NEON

MBAN MPEO

MPK

MPC MPXYZ MPID

EO (NNOD) C (NEQN. MBAN) K(NEQN. MBAN)

ID (NNOD)

XYZ (NNOD. 3)

DESCRIPTION

TO SUBROUTINE DO-WHILE TERMINATOR FLAG

NUMBER OF NOOES IN SYSTEM NUMBER OF (SYSTEM) EQUATIONS

(HALF) BANDWIDTH OF SYSTEM EQUATIONS POINTER TO EO (NNOD) IN BLANK COMMON POINTER TO K (NEQN. MBAN) IN BLANK COMMON POINTER TO C (NEQN. MBAN) IN BLANK COMMON POINTER TO XYZ (NNOD. 3) IN BLANK COMMON POINTER TO ID (NNOD) IN BLANK COMMON INITIAL EXCITATION VECTOR SYSTEM CAPACITY MATRIX (COMPACT FORM) SYSTEM CONDUCTANCE TRANSFER MATRIX (COMPACT FORM) EQUATION NUMBER/B.C. FLAG ARRAY; ABS (ID (N) ) = EQUATION NUMBER OF NODE N SIGN (ID (N) ) - -1 - TEMPERATURE PRESCRIBED NODE SIGN (ID (N)) - +1 - FLUX PRESCRIBED NODE COORDINATE ARRAY (ORDERED BY NODE NUMBER)

CAL-SAP: DATA 4 COMMON STORAGE

COMMON MTOT.NP, IA(1) COMMON /DBSYS/ NUMA. NEXT. IDIR. IP (3) COMMON /IOLIST/ NTM. NTR. NIN . NOT. ND1 . ND2 . ND3 . ND4 C

C RFORM: DATA 4 COf*)ON STORAGE

C LOGICAL*! ERR. NOELEM

IF (ENDFLAG. EQ. ’END’) OO TO 400 CALL FREEK’ ’.NODE.l)

IF (NODE.EQ. 0) GO TO 400 CALL FREER ( »X\ X, 1) CALL FREER (’Y’.Y.l) CALL FREER (

' Z * , Z. 1) C STORE COORDINATES

XYZ (NODE, 1) - X

XYZ (NODE. 2) « Y XYZ (NODE. 3) «= Z C READ AND SET GENERATION PARAMETERS

NN (1) - 0 NN (3) -0

CALL FREE I

( 'N ’ , NN , 3) Nl - NN (1

)

IF(Nl.EQ.O) GO TO 300

SAVE /DBSYS/, /IOLIST/

CC—0.0 WRITE HEADER

C WRITE (NTM. 2000) WRITE (NOT. 2000) 2000 FORMAT (/’ -= ELEMENT DATA*)

CC—1.0 DEFINE ARRAYS

C CALL DELETE (’C *) CALL DELETE (’K ') CALL DELETE CEO ') CALL DEFINE (

’ EO ’. MPEO, NNOD. 1) CALL DEFINECK ’. MPK. NEQN . MBAN

)

CALL DEFINE (’C ’, MPC, NEQN . MBAN

)

D - 3
cC—2.0 INITIALIZE ARRAYS

C CALL ZEROR (IA (KPEO) .NNOO. 1) CALL ZEROR (IA (MPK) , NEON , MEAN) CALL ZEROR (IA (KPC) , NEQN , MBAN

)

CC—3.0 FIND ELEMENTS

C NOELEM -.TRUE. MAXBAN - 0 CC—3.1 FIND SEPARATOR "RESIST"

C CALL FINDN ( 'RESIST' , 6, KEY) IF (KEY . EQ. 0) THEN WQRLEM- . FALSE

.

CALL FREETY

CALL RES I ST (IA (HP ID) . IA (MPK) . NNOD, NEQN . MBAN, MAXBAN. ’RES' , ERR) ENDIF C

C—3 . 2 FIND SEPARATOR "rLOWLO"

C CALL FINDN ( 'FLOWLO* , 6, KEY) IF (KEY . EQ. 0) THEN

JKDELEN= . FALSE

.

CALL FREETY

CALL RESIST (IA(MPID) . IA (MPK) , NNOD. NEQN . MBAN . MAXBAN . ’FLO’ .ERR) END IF

CC—3.3 FIND SEPARATOR - CAPACI

" C CALL FINDN (’CAPACI ’ .6. KEY) IF(KET.EQ.O) THEN GELEM= . FALSE

.

CALL FREETY

CALL CAPACI (IA(MPID) .IA(MPC) . NNOD. NEQN . MBAN. ERR) ENDIF

CC—3.4 FIND SEPARATOR "ISOID"

C CALL FINDN (* ISOID’ .5, KEY) IF (KEY . EQ. 0) THEN N€ELEM= . FALSE

.

CALL FREETY

CALL ISOLD (IA(MP ID) .IA(MPXYZ) .IA(MPC) , IA (MPK) , IA (MPEO) ,

. NNOD. NEQN. MBAN, MAXBAN, ERR) ENDIF

CC—3.5 FIND SEPARATOR "IS02D4"

C CALL FINDN ( ’IS02D4 '

, 6. KEY) IF (KEY . EQ. 0) THEN SCELEM* . FALSE

.

CALL FREETY

CALL IS02D4 (IA (MPID) . IA (MPXYZ) . IA (MPC) , IA (MPK) , IA (MPEO) .

. NNOO. NEQN. MBAN. MAXBAN. ERR) ENDIF

CC—3.6 FIND SEPARATOR "VNMRT"

C CALL FINDN ('VNMRT' .5, KEY) IF (KEY .EQ. 0) THEN ROELEM=. FALSE.

CALL FREETY

CALL VNMRT (IA (MPID) , IA(MPK) . NNOD. NEQN , MBAN . MAXBAN . ERR; ENDIF

C

C—4 . 0 IF BO ELEMENTS ... C XFCNOCLEM) THEN WRITE (NTH. 2400)

WRITE (BOT, 2400)

ERR-.TRDE. ENDIF

C

C—5 . 0 REPORT BANDWIDTH NEEDED

C WRITE (NTM. 2500) MAXBAN . MBAN WRITE (NOT, 2500) MAXBAN. MBAN

RETURN

2400 FORMAT (' ERROR: No element data found.)

2500 FORMAT </’ — NOTE: Needed bandwidth ',15./ Requested bandwidth ..’,15)

END

RESISTSUBROUTINE RESIST (ID, K. NNOD. NEQN . MBAN . MAXBAN . TYPE. ERR) C-SUB : RESIST - FORMS 4 ASSEMBLES RESISTANCE OR FLOWLOOP ELEMENTS

C IMPLICIT REAL* 8 (A-H.O-Z)

C

C

C CAL -SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM. NTR, NIN . NOT . ND1 . ND2 , ND3 . ND4

CC RESIST: DATA 4 COWON STORAGE

C REAL *8 K (NEQN, MBAN) INTEGER ID (NNOD) , LMNEN (2) , LMOLD (2) LOGICAL* 1 ERR CHARACTER ENDFLAG* 3 . TYPE*3

SAVE KNEW, NOLD, LMNEM, LMOLD. R. A, /IOLIST/

NDOF - 2

CC—1.0 WRITE ELEMENT HEADER

C

IF (TYPE. EQ. 'RES’ ) THEN WRITE (NOT, 2100) ELSEIF (TYPE . EQ. ' FLO' ) THEN WRITE (NOT. 2110) ENDIF

2100 FORMAT

(

. /’ — RESISTANCE ELEMENTS’//

.' Elern I-Node J-Node Resist. Area') 2110 FORMAT

(

./’ -= FLOW-LOOP ELEMENTS’//

.* Elern I-Node J-Node Flow Rate Capac.')

C

C —2.0 GET FIRST ELEMENT, FORM 4 ADD ELEMENT TO SYS

C INCR « 0 R - 0.0 A - 0.0 CALL FREE

CALL FREETY CALL FREEI (* ’ .NOLD. 1) CALL FREEI (

' I

' , LMOLD (1 ) , 2) IF (TYPE. EQ. ’RES’ ) THEN CALL FREER CR’ ,R. 1) CALL FREER (’A* , A, 1) ELSEIF (TYPE. EQ. ’FLO' ) THEN CALL FREER (’W’ ,R, 1) CALL FREER (’C’ .A. 1) ENDIF

CALL RESIS0 (ID. NOLD. LMOLD . K . R . A , NNOD . NEQN . MBAN . MAXBAN. TYPE. ERR)

CC—3.0 GET NEXT LINE OF DATA

C

30 CALL FREE

CALL FREETY

C CHECK FOR "END" CALL FREEH (’ '. ENDFLAG. 3 . 1

)

IF (ENDFLAG. EQ. ’END’ ) THEN RETURN ENDIF

C GET NEW ELEMENT INFORMATION CALL FREEI (* ’.NNEW.l) CALL FREEI ('X'.LMJEWU) ,2) CALL FREEI ('N’ . INCR. 1) IF (INCR . EQ. 0) INCR=1

IF (TYPE. EQ. 'RES’ ) THEN CALL FREER ('R' ,R. 1) CALL FREER (’A' , A. 1) ELSEIF (TYPE. EQ. ’FLO’) THEN CALL FREER <’W’ ,R, 1) CALL FREER ('C' , A, 1) ENDIF

C CHECK NUMERICAL ORDER

IF (NNEW . LE . NOLD) THEN WRITE (NTM. 2300) NNEW

WRITE (MOT. 2300) NNEW

ERR- . TRUE

.

RETURN

ENDIF C OTNERATE MISSING ELEMENTS

IF (NNEW . GT . NOLD+ 1 ) THEN

DO 34 N—NOLD+ 1 , NNEW- 1 ,

DO 32 I-l.NDOF

32 LMOLD (I) - LMOLD (I) INCR

34 CALL RESIS0 (ID. N. LMOLD. K.R. A, NNOD. NEQN. MBAN. MAXBAN, TYPE. ERR) ENDIF

C DO NEW ELEMENT NOLD - NNEW

DO 36 I -1, NDOF

36 LMOLD (I) - LWEW(I)

CALL RESIS0 (ID. NOLD. LMOLD. K. R. A. NNOD. NEQN. MBAN. MAXBAN. TYPE. ERR)

D 4
GO TO 30 NDOF - 1 2300 FORMAT (' *•** ERROR: Element number ’.IS,* Is out of order.*)

END

RESISO SUBROUTINE

RESISO (ID, NELM. LM.K. R, A. NNOD, NEQN . KBAN . MAXBAN . TYPE . ERR) C-SUB : RESISO - REPORTS ELEM INFORMATION. CHECKS BANDWIDTH,

C FORMS ELEM ARRAYS 4 ADDS THEM TO SYS ARRAYS

C IMPLICIT REAL *8 (A-H.O-Z) C

CC CAL -SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM, NTR, NIN . NOT, ND1 , ND2 , ND3 , ND4

CC RESISO: DATA 4 COMMON STORAGE

C REAL* 8 K (NEQN . MBAN ) . S (2, 2) INTEGER ID (NNOD), LM(2) LOGICAL* 1 ERR CHARACTER TYPE*

3

SAVE /IOLIST/

C

C REPORT ELEMENT INFORMATION TO OUTPUT DATA FILE

C WRITE (NOT. 2000) NELM.LM(l) ,LM(2) ,R.A

CC NODE ERROR TRAP

C DO 200 N=1 .

NN*LM(N)

IF (NN.LE. O.OR.NN.GT.NNOD) THEN WRITE (NTM. 2010) NN

WRITE (NOT. 2010) NN

ERR= . TRUE

RETURN

200 ENDIF

C

C DETERMINE ELEMENT BANDWIDTH

C CALL MBAND (ID, LM, 2. NNOD. MBAN, MAXBAN, ERR)

C

C FORM ELEMENT ARRAYS

C IF (TYPE. EQ. ’FLO*) THEN

A = R*A

R = 1.0ENDIF

CALL RESIS (S, R. A) C

C ADD ELEMENT CONTRIBUTION TO SYSTEM ARRAYS

C CALL ADDCM (ID. LM. S , K, 2 , 2 . NNOD, NEQN , MBAN)

CC—1.0 WRITE ELEMENT HEADER

C WRITE (NOT, 2100) 2100 FORMAT

(

./* CAPACITANCE ELEMENTS’//

Elere I-Node Mass Capac.')

C

C—2.0 GET FIRST ELEMENT. FORM 4 ADD ELEMENT TO SYS

C INCR « 0 MASS « 0.0 CAP - 0.0 CALL FREE CALL FREETY CALLFREEK* *.NOLD,l) CALL FREEI (

* I ’ . LMOLD , 1

)

CALL FREER ( *M’ .MASS, 1) CALL FREER ( 'C* .CAP, 1)

CALL CAPAC0 (ID. NOLD. LMOLD. C, MASS. CAP. NNOD. NEQN, MBAN, ERR)

CC—3.0 GET NEXT LINE OF DATA

C

30 CALL FREE

CALL FREETY

C CHECK FOR "END- CALL FREEH (’

’ .ENDFLAG. 3. 1) IF (ENDFLAG.EQ. ’END’ ) THEN

RETURN

ENDIF

C GET NEW ELEMENT INFORMATION CALLFREEI(’ *,NNEW,1) CALL FREEI (

' I

' , LMNEW, 1

)

CALL FREEI CN’ . INCR. 1) IF (INCR . EQ. 0) INCR-1 CALL FREER CM’ .MASS, 1) CALL FREER CC' .CAP, 1) C CHECK NUMERICAL ORDER

IF (NNEW.LE.NOLD) THEN WRITE (NTM. 2300) NNEW

WRITE (NOT, 2300) NNEW

ERR= . TRUE

.

RETURN

ENDIF

C GENERATE MISSING ELEMENTS

IF (NNEW . GT . NOLD+1) THEN DO 34 N*NOLD+l , NNEW-1 ,

LMOLD « LMOLD INCR

34 CALL CAPAC0 (ID,N. LMOLD. C. MASS. CAP. NNOD, NEQN. MBAN. ERR) ENDIF

C DO NEW ELEMENT NOLD - NNEW LMOLD * LMNEW CALL CAPAC0 ( ID. NOLD, LMOLD. C , MASS . CAP , NNOD , NEQN. MBAN . ERR)

RETURN GO TO 30

2000 FORMAT (3 (5X, 15) .4X.2G10.3) 2010 FORMAT (’ ••** ERROR: (Generated) Node *,15.* Is out of range.’)

END

C RESIS SUBROUTINE RESIS (XK, R, A) C SUB: RESIS - 2-NODE CONVENTIONAL RESISTANCE ELEM

C IMPLICIT REAL *8 (A-H.O-Z) C DIMENSION XK (2,2) C FORM 2X2 CONDUCTANCE MATRIX XK (1,1) - A/R XK (1,2) - -A/R XK (2,1) - -A/R XK (2,2) - A/R RETURN

END

CAPACI SUBROUTINE CAPACI (ID. C. NNOD. NEQN . MBAN. ERR) C-SUB: CAPAC I - FORMS 4 ASSEMBLES CAPACITANCE ELEMENTS

C IMPLICIT REAL* 6 (A-H.O-Z) C

• c

C CAL-SAP: DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM. NTR. NIN . NOT. ND1 . ND2 . ND3 . ND4

C

C CAPACI: DATA 4 COMMON STORAGE

C REAL* 8 C (NEQN. MBAN) . MASS. CAP INTEGER ID (NNOD). L8tJEW, LMOLD

LOG I CAL *1 ERR CHARACTER ENDFLAG*

3

SAVE NNEW, NOLD, LHJEW, LMOLD. MASS . CAP . /IOLIST/

2300 FORMAT ( ’ •••• ERROR: Element number '.15. * Is out of order.')

END

C -CAPAC0 SUBROUTINE CAPAC0 (ID . NELM. LM, C. MASS . CAP . NNOD. NEQN , MBAN. ERR) C-SUB :CAPAC0 - REPORTS ELEM INFORMATION.

C FORMS ELEM ARRAYS 4 ADDS THEM TO SYS ARRAYS

C IMPLICIT REAL* 8 (A-H.O-Z) C

C

C CAL-SAP: DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM. NTR . NIN , NOT. ND1 , ND2 . ND3 . ND4 C

C FORM: DATA 4 COMMON STORAGE

C REAL *8 C (NEQN, MBAN) , HASS, CAP, HC INTEGER ID (NNOD). LM .NELM

LOG I CAL*1 ERR

SAVE /IOLIST/

C

C REPORT ELEMENT INFORMATION TO OUTPUT DATA FILE

C WRITE (NOT. 2000) NELM, LM, MASS, CAP C

C NODE ERROR TRAP

C NN-LM

IF (NN.LE. O.OR.NN.GT.NNOD) THEN WRITE (NTM, 2010) NN

WRITE (NOT, 2010) NN

ERR*. TRUE.

RETURNENDIF D -5
FORM CLEMENT ARRAY

C PR (2) - 0.0D0

c

c

c

c

c

HC - MASS • CAP

ADD ELX*»JT CONTRIBUTION TO SYSTEM ARRAYS

CALL ADOCM(ID.LM.MC.C.1,1.NNOD,NEQN.MBAN)

RETURN

PR (3) - 1.0D0

ELSEIF (NDOF.EQ.4) THEN

PR (1 ) - -l.ODO

PR (2) - -1.0DO/3.0DO PRO) - 1.0D0/3.0D0

2000 FORMAT (2 (5X, 15) . 4X, 2G10.3) 2010 FORMAT ( * •••• ERROR: (Generated) Node '.15.' Is out of range.')

END

CC—1.0 WRITE ELEMENT HEADER

C WRITE {KOT, 2100) 2100 FORMAT

(

./• «= ISOPARAMETRIC ID 2 TO 4 -NODE CONDUCTION ELEMENTS'//

Kx « Conductivity P * Density C = Capacity'/ A * Area Q = Int . Heat’//

.* Ele® *©de Coords. Kx P C A Q’

.* Q»j*d . Points')

CC—2.0 GET FIRST ELEMENT, FORM £ ADD ELEMENT TO SYS

C

INCR * 0 DO 10 1*1.

10 PX(I) « 0.0 PK « 0.0

PP = 0.0 PC = 0.0

PA * 0.0 PQ » 0.0

CALL FREE

CALL F1JEETY CALL FREEI (

' '.MOID.l) DO 15 1*1.4

15 LMOLD (I) « 0 CALL FREEI ('I' ,LMOLD(l) ,4) WDOr * 4 IF(LMOLD(4) .LZ.0) NDOF-3 SF(LMOLD(3) .LE.0) WDOF-2 CALL FREER (

* X ' , PX ( 1 ) . 4

)

CALL FREER ( ' K * ,PK, 1) CALL FREER (*P* , PP, 1) CALL FREER ('C. PC. 1) CALL FREER ( 'A* , PA, 1) CALL FREEACQ'.FQ.l)

QUADFLAG *12345’ CALL FREEH CR\ QUADFLAG. 5.1) GAUSS - .FALSE. IF (QUADFLAG. EQ. 'GAUSS') GAUSS - .TRUE. IF ( .WOT.GAUSS) THEN C SET DEFAULT QUAD POINTS

IF ‘WDOF.EQ.2) TfiEN

PR ( 4 ) - l.ODO

ENDIF

C GET QUAD POINTS FOR EVALUATION OF [C

)

CALL FREER <’R' .PR (1) .NDOF)

DO 20 1*1. NDOF

IF (ABS (PR(I) ) .GT. 1 .0D0) THEN

WRITE (NTH. 2200)

WRITE (NOT. 2200)

20 ENDIF

ENDIF

2200 FORMAT (

' WARNING: Quadrature points should normally

be In the rar^ge of -1.0 to **-1.0')

CALL ISO1D0 (ID. NOLD. LMOLD.C, K. EO. PX. PK. PA. PP. PC. PQ. PR, 4NDOF . GAUSS , NNOO . NEON . KBAN . MAX BAN . ERR

)

CC--3.0 GET NEXT LINE OF DATA

C

30 CALL FREE CALL FREETY

C CHECK FOR "END- CALL FREEH (

'

' . ENDFLAG. 3.1)

IF (ENDFLAG . EQ. ' END ' ) THEN RETURN

ENDIF

C GET NEW ELEMENT INFORMATION CALLFREEIC '.NNEW.l) DO 35 1=1. 4 35 LMNEW(I) = 0 CALL FREEI (

' I

’ . LMNEW (1 ) , 4) NDOF = 4 IF (LMNEW (4) .LE. 0) NDOF=3

IF (LMNEW (3) .LE.0) NDOF =2 CALL FREEI ('N' . INCR. 1) IF(INCR.EQ.O) INCR=1 CALL FREER ('X' . PX(1) , 4) CALL FREER ('K' . PK. 1) CALL FREER ('P' .PP. 1) CALL FREER CC' .PC. 1) CALL FREER ('A' ,PA,1) CALL FREER ('Q' .PQ. 1) QUADFLAG - '12345' CALL FREEH ('R' .QUADFLAG. 5. 1) GAUSS = .FALSE.

IF (QUADFLAG. EQ. ’GAUSS') GAUSS = .TRUE.

IF (.NOT. GAUSS) THEN

C SET DEFAULT QUAD POINTS

IF (NDOF. EQ. 2) THEN

PR (1) * - SQRT (2 . 0D0/3 . 0D0)

PR (2) - - PR(1)

ELSEIF (NDOF. EQ. 3) THEN

PR (1) - -1 . ODO

PR (2) - 0.0D0 PRO) - l.ODO

ELSEIF (NDOF.EQ.4) THEN

PR(1) - -1 . ODO

PR (2) - -1.OD0/3.OD0 PRO) - 1 . 0D0/3 . ODO

PR (4) - l.ODO

ENDIF C GET QUAD POINTS FOR EVALUATION OF (C) CALL FREER ('R' .PR (1) .NDOF)

ISOID SUBROUTINE ISOID (ID. XYZ.C, K. EO. NNOD, NEQN, MBAN. MAXBAN , ERR) C-SUB: ISOID - FORMS £ ASSEMBLES ISOPARAM.1D 2-4-NCOE FINITE ELEMENTS

C IMPLICIT REAL* 8 (A-H.O-Z) C

C

C CAL-SAP: DATA £ COMMON STORAGE

C COMMOM /IOLIST/ NTM . NTR . NIN . NOT . ND1 , ND2 , ND3 . ND4 C

C ISOID: DATA £ COMMON STORAGE

C

C DICTIONARY OF VARIABLES

C

C VARIABLE

C PK

C PA

C PP

C PC

C PQ

C PX (4) C PR (4) C* ••**•***••••'

REAL* 8 K (NEQN. MBAN) . C (NEQN , MBAN ) ,XYZ (NNOD, 3) , EO (NNOD)

.

4 PK. PA. PP . PC, PQ,PX(4) ,PR(4) INTEGER ID (NNOD). L«4EW(4), LMOLD ( 4 ) LOGICAL* 1 ERR. GAUSS CHARACTER ENDFLAG* 3. QUADFLAG*

5

SAVE RSEW. NOLD. LWiEW, LMOLD, PK, PA.PP. PC. FQ. PX, PR, GAUSS. /IOLIST/

DESCRIPTION

ELEMENT CONDUCTIVITY AREA AVAILABLE FOR HEAT TRANSFER ELEMENT DENSITY ELEMENT CAPACITY ELEMENT INTERNAL HEAT GENERATION RATE/LENGTH ELEMENT NODE LOCAL COORDINATES

QUAD POINTS -1.0 TO 1.0

PR (1) «= - SORT (2. 0D0./3. 0D0) DO 40 1=1. NDOF

PR (2) * - PRO ) IF (ABS (PR (I) ) .GT. l.ODO) THEN

EL-SEIF (NDOF . EQ. 3) THEN WRITE (NTM, 2200)

PRO) * -1-000 WRITE (NOT. 2200)

40 ENDIF D 6
ENDIF

C

C CHECK NUMERICAL ORDER IS0241

IF (NNEW. LE. NOLD) THEN

IfRITE (NTH, 2300) NNEW

WRITE (NOT. 2300) NNEW

ERR= . TRUE

.

RETURN

ENDIF

C GENERATE MISSING ELEMENTS

IF (NNEW .GT . NOLD+1 ) THEN DO 50 N-NOLD+l.NNEW-1,1

DO 45 I-l.NDOF

45 LMOLD(I) - LMOLD (I) INCR

50 CALL I SOI DO (ID. N , LMOLD, C . K . EO, PX, PK. PA, PP , PC, PQ, PR,

+ NDOF. GAUSS. NNOD. NEON. MBAN. MAXBAN. ERR) ENDIF

C DO NEW ELEMENT NOLD - NNEW

DO 55 I-l.NDOF

55 LMOLD (I ) - LMNEW(I) CALL ISO1D0 (ID. NOLD , LMOLD. C, K, EO. PX . PK, PA, PP , PC. PQ. PR, NDOF, GAUSS, NNOD.NEQN. MBAN, MAXBAN , ERR)

SUBROUTINE IS02 41 (KE. CE. EE. PK, PA. PP. PC. PQ. PX , NELNOO, RPT, GAUSS

)

C--SUB : IS0241 - 2 TO 4 NODE ISOPARAMETRIC CONDUCTION FINITE ELEMENT

NODE NUMBERING

C DICTIONARY OF VARIABLES

C VARIABLE

C DUMMY VARIABLES

C NELNOO

C NIP

C KE (NELNOD. NELNOD)

C CE (NELNOD. NELNOD)

C EE (NELNOD) PK

PA

PPPCPQ

PX (NELNOD) RPT ( 4

)

DESCRIPTION

NUMBER OF ELEMENT NODES (2 TO 4] NUMBER OF INTEGRATION POINTS ELEMENT CONDUCTANCE MATRIX ELEMENT CAPACITANCE MATRIX ELEMENT INTERNAL GENERATION RATE VECTOR CONDUCTIVITY AREA AVAILABLE FOR HEAT TRANSFER WEIGHT DENSITY

HEAT CAPACITY

INTERNAL GENERATION RATE PER UNIT VOL.

(GLOBAL) COORDINATES OF ELEM. NODAL POINTS

QUAD. POINTS USED TO EVALUATE [C] C GAUSS GO TO 30

2300 FORMAT (' BRROR : Element number QUADRATURE METHOD TO EVALUATE (C) Is out of order.') C

END

C ISO1D0 C SUBROUTINE ISO1D0 (ID . NELM. LM. C, K. EO. PX, PK. PA.PP. PC. PQ. PR. NDOF. GAUSS. NNOD.NEQN. MBAN. MAXBAN. ERR) C-SUB : ISO1D0 - REPORTS ELEM INFORMATION. CHECKS BANDWIDTH.

C FORMS ELEM ARRAYS 4 ADDS THEM TO SYS ARRAYS

C IMPLICIT REAL* 8 (A-H.O-Z) C

cC CAL-SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM, NTR, NIN , NOT, ND1 , ND2 . ND3 , ND4

C

C ISO1D0 : DATA 4 COMMON STORAGE

C REAL* 8 C (NEON. MBAN) . K (NEQN . MBAN ) .EO(NNOD) . CE (4 . 4 ) . KE (4 . 4 )

.

EE (4) . PX ( 4 ) .PR (4) INTEGER ID (NNOD). LM(4) LOGICAL* 1 ERR. GAUSS SAVE /IOLIST/

C

C REPORT ELEMENT INFORMATION TO OUTPUT DATA FILE

C

IF (GAUSS) THEN WRITE (NOT. 2000) NELM, LM (1 ) ,PX(1) , PK, PP . PC. PA. PQ

WRITE (NOT. 2005) (LM(I).PX(I) ,1-2. NDOF) ELSE

WRITE (NOT. 2010) NELM, LM (1 ) ,PX(1) . PK. PP . PC , PA. PQ. (PR(I) . I-l.NDOF)

WRITE (NOT. 2005) (LM (I ) . PX (I ) . 1=2 . NDOF) ENDIF

2000 FORMAT (215. 6G8.2E1. ’GAUSS’) 2005 FORMAT (5X.I5.G6.2)

2010 FORMAT (215, 6G8 . 2E1 . 4F6. 3) C

C NODE ERROR TRAP

C DO 200 N-l.NDOF

NN-LM (N) IF (NN.LE. O.OR.NN.GT.NNOD) THEN WRITE (NTM. 2020) NN

WRITE (NOT, 2020) NN

ERR-. TRUE.

RETURN

200 ENDIF

2020 FORMAT (' ERROR: (Generated) Node '.IS,' Is out of range.')

C

C DETERMINE ELEMENT BANDWIDTH

C CALL MBAND (ID. LM.NDOF, NNOD, MBAN. MAXBAN, ERR) C

C FORM ELEMENT ARRAYS

C CALL IS0241 (KE, CE. EE, PK, PA.PP. PC, PQ.PX, NDOF, PR, GAUSS)

C

C ADD ELEMENT CONTRIBUTION TO SYSTEM ARRAYS

C CALL ADDCM (ID, LM, KE, K, NDOF. 4 . NNOD, NEQN . MBAN

)

CALL ADDCM ( ID, LM.CE.C, NDOF, 4. NNOD. NECK. MBAN)

DO 300 N-l.NDOF

300 EO (LM (N) ) = EE (N)

RETURN

END

- .TRUE. USE GAUSS QUADRATURE

- .FALSE. USE GIVEN QUAD POINTS

C INTERNAL VARIABLES

C QPT ( 4

)

GAUSS QUAD. POINTS

C QWT (4) GAUSS QUAD. WEIGHTS

C H ( 4

)

ELEMENT SHAPE FUNCTION VALUES

C DEL (4) ELEMENT LOCAL DERIVATIVE VALUES

C R LOCAL COORDINATE

C JACOB JACOBIAN

C IMPLICIT REAL* 8 (A-H.O-Z) REAL* 8 KE (4. 4) . CE (4. 4) . EE (4) . PX (4) , RPT(4) REAL*8 JACOB, QPT(4), QWT(4), H (4) . DEL ( 4

)

LOGICAL* 1 GAUSS

C-1.0 ZERO ELEMENT ARRAYS

DO 10 1=1.4

EE (I ) = 0.0D0

EX) 10 J=1 ,

KE(I.J) = 0.0D0

CE(I.J) = 0.0D0

10 CONTINUE

C-2.0 EVALUATE CONDUCTANCE MATRIX AND INTERNAL GENERATION RATE VECTOR

C-2.1 GET GAUSS QUAD POINTS AND WEIGHTS

C NOTE: USE NIP-NELNOD-1 FOR EXACT INTEGRATION OF [KE] 4 (EE) NIP - NELNOD - 1 CALL GAUSCO (NIP , QPT , QWT) DO 30 IP-1, NIP C-2.2 GET SHAPE FUNCTION VALUES AT QUAD POINT CALL SHP241 (QPT (IP) ,H. NELNOO) C-2.3 GET LOCAL DERIVATIVE VALUES AT QUAD POINT CALL DER241 (QPT (IP) .DEL, NELNOD) C-2.4 COMPUTE JACOBIAN AT THE POINT (NOTE: FOR ID JACOBIAN IS 1X1) CALL JACOB 1 (NELNOD . DEL . PX . JACOB

)

C-2.5 ACCUMULATE QUADRATURE CONTRIBUTION TO ELEMENT MATRICES DO 20 1-1 . NELNOO

C ELEMENT INTERNAL GENERATION RATE VECTOR

EE (I) - EE (I) QWT (IP) *PQ*H (I) ‘JACOB DO 20 J-l, NELNOO

C ELEMENT CONDUCTANCE MATRIX

KE(I.J) - KE(I.J) QWT (IP) *DEL (I) *PK*DEL ( J) / JACOB

20 CONTINUE

30 CONTINUE

C-3.0 EVALUATE CAPACITANCE MATRIX

C-3.1 GET QUAD POINTS AND WEIGHTS

C NOTE: USE NIP^ELNOD FOR EXACT INTEGRATION OF (CJ NIP - NELNOO

C USE GAUSS QUADRATURE IF SPECIFIED

C EVALUATE HEIGHT COEFFICIENTS IF QUAD POINTS ARE SPECIFIED ir (GAUSS) CALL GAUSCO (NIP, RPT, QWT)

IF (.NOT. GAUSS) CALL QUADCO (NIP . QWT . RPT) DO 50 IP-1, NIP C-3.2 GET SHAPE FUNCTION VALUES AT QUAD POINT CALL SHP241 (RPT (IP) ,H, NELNOO) C-3.3 GET LOCAL DERIVATIVE VALUES AT QUAD POINT CALL DER241 (RPT (IP) .DEL, NELNOD) C-3.4 COMPUTE JACOB I AN AT THE POINT (NOTE: FOR ID JACOBIAN IS 1X1) CALL JACOBI (NELNOD. DEL. PX. JACOB) C-3.5 ACCUMULATE QUADRATURE CONTRIBUTION TO ELEMENT MATRICES

DO 40 I -1, NELNOD

DO 40 J-l, NELNOD

CE (I, J) - CE(I.J) QWT(IP) *H (I)*PP*PC*H(J) ‘JACOB D -7
40 CONTINUE

50 CONTINUE

RETURN

END

CSHP241 SUBROUTINE SHP241 (R.H, NELNOD) C-SUB : SRP24L - DETERMINES THE VALUE OF 2-4 NODE ISOPARAMETRIC ID C SHAPE FUNCTIONS, H, AT LOCAL COORDINATE R

C

C

NODE NUMBERING 1 3 4 2 DICTIONARY OF VARIABLES

C VARIABLE DESCRIPTION

C R LOCAL COORDINATE -1 TO +1 C H(2 TO 4) VALUE OF SHAPE FUNCTIONS

C NELNOD NUMBER OF ELEMENT NODES l 2 TO 4

]

IMPLICIT REAL* 8 (A-H.O-Z) DIMENSION H (NELNOD)

H (1 ) - 0.5*(1.0 - R) 3(2) - 0.5* (1.0 R) IF (NELNOD. GE. 3) THEN

B (1 ) » H (1) - 0 . 5* (1 . 0 - R*R)

H<2) ~ H (2) - 0.5M1.0 - R*R)

fl(3> » (1.0 - R*R) ENDIF

IF (NELNOD. EQ. 4) THEN

3(1) « H (1) (-9 . 0*R*R*R R*R 9.0*R - 1.0)/16.0

3(2) « 3(2) > (9. 0*R*R*R + R*R - 9.0*R - 1.0)/16.0

H (3 ) - H (3) + (27 . 0*R*R*R + 7.0*R*R - 27.0*R - 7.0)/16.0

3(4) « (-27 . 0 *R*R*R - 9.0*R*R 27.0*R + 9.0J/16.0 ENDIF RETURN END

CDER241

C-SUB

C

c

c

cc~

c

SUBROUTINE DER2 4 1 (R. DEL, NELNOD)

:DER241 - DETERMINES THE VALUE OF LOCAL DERIVATIVES. DEL. OF

2-4 NODE ISOPARAMETRIC ID SHAPE FUNCTIONS AT LOCAL COORDINATE R

NODE NUMBERING 1—3 4 2 DICTIONARY OF VARIABLES

VARIABLE DESCRIPTION

R LOCAL COORDINATE -1 TO +1 DEL (2 TO 4) VALUE OF LOCAL DERIVATIVE OF SHAPE FUNCTIONS NELNOD NUMBER OF ELEMENT NODES [ 2 TO 4

)

IMPLICIT REAL* 8 (A-H.O-Z) DIMENSION DEL (NELNOD)

DEL (1 ) - -0.5 DEL (2) * 0.5 IF (NELNOD. GE. 3) THEN

DEL (1 ) - DEL (1) + R

REAL*0 JACOB DIMENSION DEL (NELNOD). XCOORD (NELNOD) JACOB - 0.0 DO 10 N«l. NELNOD

10 JACOB - JACOB DEL (N) • XCOORD (N) RETURN

END

CGAUSCO

SUBROUTINE GAUSCO (N , A. W) C—SUBrGAUSCO - GAUSS QUADRATURE ABSCISSAE AND WEIGHT COEFFICIENTS

CC DICTIONARY OF VARIABLES

VARIABLE

NA (N) M(N)

DESCRIPTION

NUMBER OF QUADRATURE POINTS ABSCISSAE

WEIGHTS

IMPLICIT REAL*0 (A-B.O-Z) DIMENSION A (N) . W(N) IF(N.EQ.l) THEN A (1) » 0.0D0

W(l) - 2.0DC ELSEIF (N .EQ. 2) THEN

A (1) = -1.0D0/SQRT(3. 0D0)

A (2) - "A (1)

W(l) - 1.0D0

W (2) 1 . 0D0 ELSEIF (N.EQ. 3) THEN A (I) - -SQRT (3. 0/5.0)

A (2) - 0.OD0

A (3) - -A(l)

W(l) - 5. 0D0/9.0D0

W(2) - 8. 0D0/9.0D0

N(3) . W(l) ELSEIF (N.EQ. 4) THEN

A (1) - -0.861136311S94053D0

A(2) « -0. 33998104358485 6D0

A (3) = -A (2)

A (4) = -A(l)

W(l) - 0.347854845137454D0

W(2) - 0 . 6521 451 54 862546D0

W(3) = W(2)

W (4) B W(l) ELSEPAUSE ' ***• ERROR: SUB: GAUSCO: Gauss coefs. not available.’ ENDIF

RETURN

END

DEL (2) « DEL (2) R

DEI. (3) « -2 0*R ENDIF

IF (NELNOD. EQ. 4) THEN DEL (1) - DEL (1) + (-27 . 0*R*R 4 2.0*R * 9.0)/16.0

DEL (2) *= DEL (2) (27.0*R*R + 2.0*R - 9.0)/16.C

D£L(3) * DEL (3) + <81.0*R*R + 14.0*R - 27.0) /16.0

DEL (4) « (-81 *R*R - 18*R 27.0)/16.0 ENDIF

RETURN

END

CQUADCO

SUBROUTINE QUADCO (N.W.R) C--SUB :QUADCO - EVALUATES WEIGHT COEFS. GIVEN QUAD POINTS

CC DICTIONARY OF VARIABLES

C

C VARIABLE DESCRIPTION

C INPUT

C N NUMBER OF QUADRATURE POINTS

C R (N) SPECIFIED QUAD POINTS

C OUTPUT

C W(N) WEIGHT COEFFICIENTS

C LOCAL

C C (4 , 4)

JACCBlSUBROUTINE JACOBI (NELNOD . DEL . XCOORD , JACOB) e-SUS: JACOBI - COMPUTES THE JACOBIAN, JACOB, FOR 2-4 NODE ISOPARA. C ELJSK. USING VALUES OF LOCAL DERIVATIVES OF SHAPE FUNCTION IN DEL

C6 »J « dx/dr - d/dr {HI H2 H3 H4)(X1 X2 X3 X4)T ; dEi/dr * DELI ee KOBE NUMBERING 1 3 4 2

C DICTIONARY OF VARIABLES

C

C VARIABLE DESCRIPTION

C* NELNOD NUMBER OF ELEMENT NODES [2 TO 4) C DEL (2 TO 4) VALUE OF LOCAL DERIVATIVE OF SHAPE FUNCTIONS C XCOORD (2 TO 4) (GLOBAL) COORDINATES OF NODES C JAC C ILLICIT REAL* 8 (A-a.O-Z)

QUADRATURE RULE EQUATION COEFFICIENTS

C IMPLICIT REAL* 8 (A-H.O-Z) DIMENSION W( 4) . R(4), C(4,4>

DO 10 I-l.N

W(I) - (1.0D0 - (-1 . 0D0) **I) /DBLE (I) DO 10 J*1,N

10 C(I, J) - R ( J) ** (1-1)

CALL CROUT (C,W,N,4,1,0)

RETURN

END

C IS02D4 SUBROUTINE IS02D4 (ID. XYZ, C, K. EO, NNOD , NEQN , MBAN , MAXBAN , ERR) C-SUB :IS02D4 - FORMS AND ASSEMBLES ISOPARAMETRIC 2D 4-NODE ELEMENTS

C IMPLICIT REAL* 0 (A-H.O-Z) D -8
IMPLICIT REAL* 8 (A-H.O-Z) C CAL-SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM. NTR. NIN , NOT. ND1 . ND2 . ND3 , ND4

C

C RESIST: DATA 4 COMMON STORAGE

C REAL* 8 K (NEQN, MBAN) . C (NEQN . MEAN ) .XYZ CNNOD, 3) .EO(NNOD) ,PK(3) INTEGER ID (NNOD). U«EW(4), LMOLD (4) LOGICAL* 1 ERR CHARACTER ENDFLAG* 3 , REG*

3

SAVE LMNEW, LMOLD, NNEW, NOLD. PK, PC, PP , PQ, PD, REG, /IOLIST/

NDOF * 4

CC—1.0 WRITE ELEMENT HEADER

C WRITE (NOT, 2100) 2100 FORMAT

(

./' «* ISOPARAMETRIC 2D 4 -NODE CONDUCTION ELEMENTS'//

. * c - Capacity P * Density 0 * Internal Heat Rate’/

. * D - Thickness for REG*PLN; Subtended angle for REG*AX I

’ /

/

. * Elem PODI J K L Rn Kyy Kxy C REG’) C

C —2 . 0 GET FIRST ELEMENT. FORM 4 ADD ELEMENT TO SYS

C PQ • 0.0 PD * 1.0 REG * * PLN

INCR « 1 CALL FREE CALL FREETY CALLFREEIC ’,NOLD,l) CALL FREEI (

’ I

’ . LMOLD ( 1 ) , 4

)

CALL FREER (’ K

’ , PK (1 ) , 3

)

CALL FREER ('C', PC. 1) CALL FREER ( 'P ’ .PP.l) CALL FREER CQ’. PQ.l) CALL FREER CD’, PD, 1) CALL FREEH ('G'. REG, 3,1)

CALL ISO2D0 (ID. XYZ . NOLD, LMOLD. K, C, EO. PK. PC. PP, PQ. PD. REG.

. NNOD, NEQN . MBAN , MAXBAN , ERR

)

CC —3 . 0 GET NEXT LINE OF DATA

C 30 CALL FREE CALL FREETY

C CHECK FOR "END- CALL FREEH ( ' ’ .ENDFLAG. 3. 1) IF(ENDFLAG.EQ. ’END’ ) THEN RETURN

ENDIF

C GET NEW ELEMENT INFORMATION CALL FREEI (’ ’ .NNEW.l) CALL FREEI (’I’ .LMNEW(l) , 4) CALL FREEI (’N’ . INCR, 1) IF (INCR .EQ. 0) INCR=1 CALL FREER ('K' ,PK (1) . 3) CALL FREER (*C’ .PC. 1) CALL FREER (’P’ .PP. 1) CALL FREER CQ’ .PQ. 1) CALL FREER CD’, PD. 1) CALL FREEH (’G’ .REG. 3, 1) C CHECK NUMERICAL ORDER

IF (NNEW. LE. NOLD) THEN WRITE (NTM, 2300) NNEW

WRITE (NOT. 2300) NNEW

ERR- . TRUE

RETURNENDIF

C GENERATE MISSING ELEMENTS

IF (NNEW .GT . NOLD+1 ) THEN DO 34 N-NOLD+1. NNEW- 1.1 DO 32 I-l.NDOF

32 LMOLD (I) * LMOLD (I) INCR

34 CALL ISO2D0 (ID, XYZ , N. LMOLD, K, C, EO, PK, PC. PP, PQ. PD, REG.

. NNOD. NEQN , MBAN . MAXBAN . ERR

)

ENDIF C DO NEW ELEMENT NOLD - NNEW

DO 36 1*1, NDOF

36 LMOLD (I) * LMNEW(I) CALL ISO2D0 (ID. XYZ . NOLD. LMOLD, K, C, EO. PK. PC, PP, PQ, PD. REG,

. NNOD , NEQN , MBAN , MAXBAN , ERR

)

GO TO 30

2300 FORMAT (’ **** ERROR: Element number ’.15.’ Is out of order.')

END

c ISO2D0 SUBROUTINE ISO2D0 (ID. XYZ . NELM, LM.K. C. EO. PK. PC, PP. PQ, PD. REG.

. NNOD . NEQN . MBAN . MAXBAN , ERR

)

C-SUB : ISO2D0 - REPORTS ELEM INFORMATION. CHECKS BANDWIDTH,

C FORMS ELEM ARRAYS 4 ADDS THEM TO SYS ARRAYS

C

c

c

C CAL -SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM, NTR, NIN. NOT. ND1 . ND2 , ND3 . ND4 C

C FORM: DATA 4 COMMON STORAGE

C REAL* 8 XY Z (NNOO , 3 ) , K (NEQN . MBAN ) , C (NEQN. MBAN) , EO (NNOD)

.

.PK<3) ,KE (4. 4) ,CE (4. 4) .EE (4) . XY (2. 4) INTEGER ID (NNOD). LM(4) LOGICAL* 1 ERR CHARACTER REG*

3

SAVE /IOLIST/

C

C REPORT ELEMENT INFORMATION TO OUTPUT DATA FILE

C WRITE (NOT. 2000) NELM. LM. PK. PC, PP. PQ. PD. REG

2000 FORMAT (6X. 14, IX, 413, IX, 7G8 . 2, A3) CC NODE ERROR TRAP

C DO 200 N-1,4

NN-LM(N)

IF (NN . LE . 0 . OR . NN . GT . NNOD ) THEN WRITE (NTM, 2010) NN

WRITE (NOT. 2010) NN ERR* . TRUE

.

RETURN

200 ENDIF

CC DETERMINE ELEMENT BANDWIDTH

C CALL MBAND (ID. LM. 4, NNOO. MBAN. MAXBAN. ERR)

C

C FORM ELEMENT COORD SUB-ARRAY XY ( 2 . 4

)

C DO 300 1*1,2

DO 300 J*1 ,

300 XY(I.J) « XYZ (LM ( J) , I) C

C FORM ELEMENT ARRAYS

C CALL IS02D (KE.CE, EE. PK. PC.PP, PQ. PD. REG. XY) C

C ADD ELEMENT CONTRIBUTION TO SYSTEM ARRAYS

C CALL ADDCM (ID. LM. KE, K. 4 , 4 .NNOD, NEQN. MBAN) CALL ADDCM (ID, LM, CE, C, 4 . 4 .NNOD. NEQN . MBAN)

DO 400 N=1 ,

400 EO (LM (N) ) - EE (N)

RETURN

2010 FORMAT ( ’ **•* ERROR: (Generated) Node *.15,’ Is out of range.')

END

c IS02D SUBROUTINE IS02D (KE, CE. EE.PK. PC. PP. PQ, PD. REG. XY) C-SUB :IS02D - 2D ISOPARMETRIC 4-NODE CONDUCTION FINITE ELEMENT

C MODIFICATION OF SUBROUTINE DEVELOPED BY R. TAYLOR

C U.C. BERKELEY FOR PROGRAM *HEAT*

C

C DICTIONARY OF VARIABLES

C

C VARIABLE

C KE ( 4 , 4 ) C CE (4,4) C EE ( 4 ) C PK (3) C PC

C PP

C PQ

C PD C REG

C XY (2,4) C IMPLICIT REAL* 8 (A-H.O-Z) CHARACTER REG*

3

REAL* 8 SG (4) ,TG(4), KE (4 , 4) , CE (4 , 4) , EE (4) , XY (2,4) , PK (3) COMMON / ISODAT/SH (3 , 4) , DV SAVE /ISODAT/

DATA SG/1 . , 1 .

. -1 .. -1 . /. TG/-1 . . 1 . , 1 .

, -1 .

/

G « 1. 0/SQRT (3.0) DO 100 1-1,4

EE (I) - 0.0 DO 100 J-1.4

CE(I.J) - 0.0

100 KE(I.J) * 0.0 PCPP*PC*PP

DO 103 L—1 ,

CALL I SOS HP (SG (L) *G, TG (L) *G. XY) RR * 0.0 DV * DV*PD

DO 101 1-1.4

101 RR * RR + XY (1 . 1) *SH (3,1) IF (REG.EQ. 'AXI ’ ) DV * DV*RR

DO 102 J-1,4 SHJ - SH(3,J)*DV

DESCRIPTION

ELEMENT CONDUCTANCE MATRIX ELEMENT CAPACITANCE MATRIX ELEMENT INTERNAL GENERATED FLUX VECTOR

ELEMENT CONDUCTIVITY TENSOR: Kxx, Kyy. Kxy ELEMENT SPECIFIC HEAT CAPACITY

ELEMENT DENSITY VOLUMETRIC INTERNAL HEAT GENERATION RATE ELEMENT THICKNESS /SUBTENDED ANGLE ELEMENT TYPE: 'PLN' OR ’AXI’ ELEMENT NOOAL COORDINATES

D -9
EE ( J) C TREE: BCOND - TREE OF SUBROUTINES TO ESTABLISH BOUNDARY CONDITIONS - EE(J) SHJ*PQ

SHJ - SRJ*PCPP

A1 - (PK(1) *SH(1. J) PK (3) *SH (2, J) ) *DV A2 - (PK (3) *SH (1, J) PK (2) *SH (2 , J) ) *DV DO 102 1-1,4 KE(I.J) - KE(I.J) A1 * SH (1,1) « A2*SH (2,1) 102 CE(I.J) - CE(I.J) SHJ*SH (3,1) 103 CONTINUE

RETURN

END

C ISOSHP SUBROUTINE ISOSHP (SS . TT, XY

)

C-SUB: ISOSHP - SHAPE FUNCTION SUBROUTINE FOR IS02D IMPLICIT REAL* 8 (A-H.O-Z) REAL* 8 SX (2,2) ,XS (2. 2) ,R (4) ,T(4) ,XY (2. 4) COMMON /ISOOAT/SH (3.4) ,DV SAVE /ISOOAT/ DATA R/-0.5, 0.5. 0.5. -0 . 5/. T/-0 . 5, -0 . 5. 0 . 5. 0 . 5/ DO 100 1-1.4

SH (3, 1 ) - (0.5+R(I)*SS)* (0. 5+T (I) *TT) SHU, I) - R (I ) • (0 . 5+T (I) *TT) 100 SH (2, 1) - T (I) • (0.5+R(I)*SS) DO 101 1-1,2 DO 101 J-1.2 XS (I, J)«0.0 DO 101 K-1,4

101 XS(Z.J) - XS(I.J) XT (I.K) *SH (J.K) DV « XS (1, 1) *XS (2. 2) - XS (1 , 2) *XS (2.1) SX (1,1) - XS (2,2) /DV

SX (2. 2) « XS (1.1) /DV SX (1, 2) - -XS (1,2) /DV SX (2. 1) - -XS (2,1) /DV DO 102 K-l. 4 CC - SH (1 , K) * SX (1,1) + SH (2 , K) *SX (2 , 1) SH (2, K) * SH (1 , K) *SX (1,2) SH (2 . K) *SX (2 . 2) 102 SH ( 1 , K) - CC RETURN

END

C VNMRT

SUBROUTINE VNMRT (ID, K, NNCO, NEON . MBAN. MAX BAN, ERR) C-SUB: VNMRT - FORMS AND ASSEMBLES VARIABLE NODE MRT ELEMENTS

C ••* PRESENTLY A STUB

RETURN

END

C MBAND SUBROUTINE MBAND (ID, LM , NDOF , NNOD . MBAN , MAXBAN , ERR) C-SUB: MBAND - DETERMINES ELEMENT BAND WIDTH

C

C CAL-SAP: DATA 4 COMMON STORAGE

cC RBCOND SUBROUTINE RBCOND (MP ID, NNOO, ERR) C-SUB: RBCOND - ROOT SUBROUTINE TO BCOND

C BCOND - ESTABLISHES BOUNDARY CONDITION FLAGS BY MODIFYING

ID (NNOD) C

C DICTIONARY OF VARIABLES

C

C VARIABLE DESCRIPTION

C VARIABLES PASSED TO SUBROUTINE

C ERR DO-WHILE TERMINATOR FLAG

C NNOD NUMBER OF NODES IN SYSTEM

C MBAN (HALF) BANDWIDTH OF SYSTEM EQUATIONS

C MPID POINTER TO ID (NNOD) IN BLANK CO***ON C ID (NNOO) EQUATION NUMBER/B.C. FLAG ARRAY;

C ABS (ID (N) ) - EQUATION NUMBER OF NODE N

C SIGN (ID (N) ) - -1 - TEMPERATURE PRESCRIBED NOOE

C SIGN (ID (N)) - +1 - FLUX PRESCRIBED NODE

C

C CAL-SAP: DATA 4 COMMON STORAGE

C COMMON MTOT.NP. IA(1) COMMON /DBSYS/ NUMA, NEXT. IDIR. IP (3) COMMON /IOLIST/ NTM, NTR, NIN , NOT . ND1 , ND2 . ND3 , ND4

CC RBCOND: DATA 4 COHION STORAGE

C LOGICAL* 1 ERR

SAVE /DBSYS/, /IOLIST/

CC--0.0 WRITE HEADER

C WRITE (NTM, 2000) WRITE (NOT, 2000) 2000 FORMAT (/' «== BOUNDARY CONDITIONS AND EQUATION NUMBERS’)

CC--1.0 FIND SEPARATOR "BOUNDARY"

C CALL FINDN (’BOUNDARY \ 8. KEY)

IF (KEY . EQ. 1 ) THEN WRITE (NTM, 2100) WRITE (NOT. 2100) RETURN

ENDIF CALL FREETY

C

cc

c

COMMON /IOLIST/ NTM, NTR, NIN, NOT, ND1 , ND2 , ND3 , ND4

MBAND: DATA

INTEGER ID (NNOD) . LM (NDOF)

CC--2.0 CALL BCOND TO DO THE WORK

C CALL BCOND (IA (MPID) . NNOD, ERR)

RETURN

MBANEL * 0 DO 10 I -1. NDOF

II - ABS (ID (LM (I) ) ) DO 10 J-l.NDOF JJ - ABS (ID (LM ( J) ) ) 10 MBANEL - MAX (MBANEL. ABS (II-JJ+1)

)

IF (MBANEL . GT . MBAN) THEN WRITE (NTM, 2000) MBANEL, MBAN

WRITE (NOT. 2000) MBANEL. MBAN

ERR-. TRUE. ENDIF

MAXBAN - MAX (MAXBAN, MBANEL)

RETURN

2000 FORMAT ( ' ***• ERROR: Bandwidth exceeded for current element.*/

. ' Bandwidth required :*,I5/ Bandwidth available: ', 15)

END

C ADDCM SUBROUTINE ADDCM (ID, LM. ELM. SYS . NDOF. NSI ZE. NNOD. NEQN . MBAN) C-SUB: ADDCM - ADDS ELEMENT MATRIX TO

C COMPACT FORM OF SYSTEM MATRIX

REAL* 8 ELM (NSI ZE, NSIZE) , SYS (NEQN, MBAN) INTEGER ID (NNOD) .LM (NSIZE)

DO 20 I -1. NDOF

II « ABS (ID (LM (I))) DO 10 J-l.NDOF JJ « ABS (ID (LM ( J) ) ) - II +1

IF ( JJ. LE. 0) GO TO 10 SYS (II , JJ) - SYS (II . JJ) ELM(I.J)

10 CONTINUE

20 CONTINUE

RETURN

END

C'

2100 FORMAT (/' — NOTE: Boundary condition data not found.'/

All nodes assumed to be flux-prescribed nodes.'/ Equation numbers = node numbers.')

END

CBCOND SUBROUTINE BCOND (ID, NNOD. ERR) C-SUB: BCOND - READS AND GENERATES BOUNDARY CONDITION FLAGS

C IMPLICIT REAL* 8 (A-H.O-Z)

C

C CAL-SAP: DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM, NTR, NIN. NOT, ND1 , ND2 , ND3 , ND4

C

C BCOND: DATA 4 COMMON STORAGE

C CHARACTER ENDFLAG*3. BC*1 DIMENSION ID (NNOD). IJK (3) LOGICAL* 1 ERR

SAVE /IOLIST/

CC—1.0 WRITE BOUNDARY CONDITIONS HEADER

C WRITE (NOT. 2100) 2100 FORMAT (/

. SX' Negative Eqtn-# - part of tenperat ure-prescrlbed boundary.'/

,6x 'Positive Eqtn-# - part of flux-prescribed boundary.*//

.SX. 'Node Eqtn-# Node Eqtn-# Node Eqtn-# Wode Eqtn-# Node Eqtn-# '

)

C

C—2.0 INTERPRET LINE OF DATA

C

20 CALL FREE

CALL FREETY

C CHECK FOR "END- CALL FREEH (' ' .ENDFLAG. 3.1) IF ( ENDFLAG . EQ . ' END ' ) GO TO 30 C DO IT CALL FREEH CC'.BC, 1,1) IF(BC.EQ. ’Q’ .OR.BC.EQ.'T') THEN

I SIGN - 1 D -1 0
IF (BC.EQ. *T* ) ISIGN « -1 CALL FREEI (’ MJK(1),3)

IF (I JK(2) .EQ. 0) IJK(2)-IJK(1)

IF(IJK(3) .EQ.O) IJK(3)-1

DO 24 N-I JK (1 ) , I JK (2 ) , I JK (3 ) IF(N.LE.O.OR.N.GT.NNOO) THEN WRITE (NTH. 2200) N WRITE (NOT, 2200) N ERR=. TRUE. GO TO 20 ENDIF

NTDOF - ABS (ID(N) ) C SET TDOF TO NEG. VALUE FOR ALL NOOES CONSTRAINED EQUAL

C TO AVOID PROBLEMS OF MULTIPLE BC SPECIFICATION USE ABS () C SO LAST REDUNDANT BC SPECIFICATION GOVERNS

WRITE (NTM, 2300) WRITE (NOT, 2300) RETURN ENDIF CALL FREETY

CC—4.0 CALL CONV TO DO THE WORK

C CALL

CONV (IA(MPID) , IA (MPXYZ) , IA (MPK) , I A (MPCH) , NNCO, NEQN , MBAN , ERR)

RETURN

2300 FORMAT (/’ — NOTE: Convection boundary data not found.')

END

DO 24 NC-l.NNOD

NCDOF-ABS (ID (NC)

)

IF (NCDOF.EQ. NTDOF) ID (NC) «ISIGN*NCDOF

24 CONTINUE

ELSEWRITE (NTM, 2210) BC WRITE (NOT, 2210) BC ERR- . TRUE

.

ENDIF GO TO 20 CC—3.0 REPORT BOUNDARY CONDITIONS IF NO ERROR ENCOUNTERED

C

30 IF (ERR) RETURN WRITE (NOT, 2300) (N . ID (N) , N-l . NNOD)

C CONV SUBROUTINE CONV (ID. XYZ. K, CH. NNOD, NEQN . MBAN. ERR) C-SUB : CONV - READS AND GENERATES CONVECTION BOUNDARY DATA

C IMPLICIT REAL* 8 (A-H.O-Z)

CC

C CAL - SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM, NTR. NIN, NOT. ND1 , ND2 , ND3 . ND4

CC CONV: DATA 4 COMMON STORAGE

C REAL* 8 K (NEQN. MBAN) . CH (NNOD) . XYZ(NNOO,3)

INTEGER ID (NNOD) , LMNEW (2) . LMDLD (2) LOGICAL* 1 ERR CHARACTER ENDFLAG* 3 . REG*

3

RETURN

2200 FORMAT (’ •••• ERROR: (Generated) Node’, 15,’ Is out of range.') 2210 FORMAT (' *••• ERROR: Boundary condition ’.Al,' not available.')

2300 FORMAT ((4X, 5(16, IX, 16, 2X) )

)

END

C TREE: CONV - TREE OF SUBROUTINES TO ESTABLISH CONVECTIVE B.C.s

SUBROUTINE RCONV (MPID, MPXYZ . MPK , MPCH. NNOD. NEQN , MBAN . ERR

)

C-SUB :RCONV - ROOT SUBROUTINE TO CONV

C

C

C-; C

ESTABLISHES CONVECTIVE BOUNDARY CONDITIONS DICTIONARY OF VARIABLES

VARIABLE DESCRIPTION--

VARIABLES PASSED TO SUBROUTINE

C ERR

C NNOD

C NEQN

C MBAN

C MPCH

C MPK

C MPXYZ

C MPID

C CH (NNOD) EXCITATION

C K (NEQN. MBAN)

C ID (NNOD) C XYZ (NNOD, 3) C

DO-WHILE TERMINATOR FLAG NUMBER OF NODES IN SYSTEM NUMBER OF (SYSTEM) EQUATIONS

(HALF) BANDWIDTH OF SYSTEM EQUATIONS POINTER TO CH (NNOD) IN BLANK COMMON POINTER TO K (NEQN) IN BLANK COMMON POINTER TO XYZ (NNOD, 3) IN BLANK COMMON POINTER TO ID (NNOD) IN BLANK COMMON

COEFS. FOR COMPUTING EFFECTIVE CONVECTIVE

SYSTEM CONDUCTANCE TRANSFER MATRIX (COMPACT FORM) EQUATION NUMBER/B.C. FLAG ARRAY; COORDINATE ARRAY (ORDERED BY NODE NUMBER) ABS (ID (N)) - EQUATION NUMBER OF NODE N SIGN (ID (N) ) - -1 * TEMPERATURE PRESCRIBED NODE

SIGN (ID (N)) - ! - FLUX PRESCRIBED NODE

C CAL-SAP : DATA 4 COMMON STORAGE

C COMMON MTOT.NP. IA (1) COMMON /DBSYS/ NUMA. NEXT. IDIR. IP (3) COMMON /IOLIST/ NTM, NTR, NIN, NOT, ND1 , ND2 . ND3 . ND4

C ICOND: DATA 4 COMMON STORAGE

C LOGICAL*! ERR

SAVE /DBSYS/, /IOLIST/

CC—0.0 WRITE HEADER

C WRITE (NTM, 2000) WRITE (NOT, 2000) 2000 FORMAT (/’ —- CONVECTIVE BOUNDARY CONDITIONS')

C

• C—1.0 DEFINE ARRAY C CALL DELETE ('CH ') CALL DEFINE ('CH MPCH. NNOD. 1

)

Cc—2.0 INITIALIZE ARRAYS

C CALL ZEROR (IA (MPCH) .NNOD, 1) CC—3.0 FIND SEPARATOR "CONVECT"

C CALL FINDNC CONVECT ’,7. KEY) IF (KEY . EQ. 1 ) THEN

SAVE

NDOF = 2

CC—1.0 WRITE CONVECTION B.C. HEADER

C WRITE (NOT, 2100) 2100 FORMAT (/

. €X ' D « Thickness for REG-PLN; Subtended angle for REG=AXI . * //

.' Surf I J H-Coef D REG') CC—2.0 GET FIRST SURFACE SEGMENT, FORM 4 MODIFY SYSTEM K C CALL FREE

CALL FREETY

CALL FREEI (' ’ .NOLD, 1) CALL FREEI ('I' ,LMOLD(l) ,2) CALL FREER ('H' .COEF, 1) D=1 .

0

CALL FREER ( 'D' .D.l) REG * ' PLN

CALL FREEH ('G'. REG. 3.1)

CALL CONVO ( ID, XYZ , NOLD, LMOLD . K. CH , COEF. REG. D. NNOD. NEQN . MBAN . ERR)

CC--3.0 GET NEXT LINE OF DATA

C 30 CALL FREE

CALL FREETY

C CHECK FOR "END- CALL FREEH ( '

' .ENDFLAG, 3.1) IF (ENDFLAG. EQ. 'END' ) THEN RETURN ENDIF

C GET NEW ELEMENT INFORMATION CALLFREEIC '.NNEW.l) CALL FREEI (

' I ' , LINEW (1 ) ,2) INCR * 1 CALL FREEI ('N' . INCR. 1) CALL FREER ('H' , COEF, 1) 0-1.0 CALL FREER ('D', D.l) REG - 'PLN' CALL FREEH CG’. REG. 3.1) C CHECK NUMERICAL ORDER

IF (NNEW . L£ . NOLD ) THEN WRITE (NTM, 2300) WNEM WRITE (NOT, 2300) NNEW

ERR-. TRUE. RETURN

ENDIF

C GENERATE MISSING SURFACE SEGMENTS

DO 34 N—NOLD+1 . NNEW- L, 1 DO 32 I -1. NDOF

32 LMOLD (I) - LMOLD (I) INCR

34 CALL CONVO (ID, XYZ. N, LMOLD, K.CH, COEF. REG, D, NNOD. NEQN, MBAN. ERR) C DO NEW SURFACE SEGMENT HOLD - NNEW

DO 36 I-l.NDOF

36 LMOLD (I) - LMJEM (I

)

CALL CONVO (ID. XYZ. NOLD, LMOLD. K. CH.COEF. REG, D. NNOD. NEQN. MBAN. ERR)

GO TO 30

2300 FORMAT (’ •••• ERROR : Surface segment number ',15,' out of order. ') D -1 1
RLSOLV

END

C CONVO SUBROUTINE CONVO (ID. XYZ . NSEG. LM. K.CH.COEF. REG. D. NNOD. NEQN. MEAN, ERR) C-SUB : CONVO - REPORTS CONVECTIVE SURFACE SEGMENT INFORMATION,

C FORMS COEFS, CH (NNOD) . FOR COMPUTING EFFECT . CONV . EXCIT

.

C MODIFIES K MATRIX FOR CONVECTION BOUNDARY CONDITIONS

C IMPLICIT REAL* 8 (A-H.O-Z) C

C CAL-SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM, NTR, NIN , NOT, ND1.ND2.ND3.ND4

C

C rORM: DATA 4 COWON STORAGE

C REAL* 8 K (NEQN , MEAN) , CH (NNOD) , XYZ (NNOD. 3) . S(2,2) INTEGER ID (NNOD) ( LM(2) LOGICAL*l ERR CHARACTER REG*

3

SAVE

CC—1.0 REPORT ELEMENT INFORMATION TO OUTPUT DATA FILE

C WRITE (NOT. 2000) NSEG. LM (1) ,LM(2) .COEF.D.REG

CC—2.0 COMPUTE COEFFICIENT. CH. FOR COMPUTING EFFECTIVE CONVECTIVE FLUX

CC ERROR TRAP FOR REGION TYPE

IF (REG. EQ. * PLN ' .OR . REG. EQ. ' AXI ' ) THEN CONTINUE

ELSEWRITE (NTM. 2200) REG WRITE (NOT. 2200) REG ERR= . TRUE

RETURN

END IF C ERROR TRAP FOR OUT-OF-RANGE NODES DO 200 N=1 ,

NN * LM (N) IF (NN . LE. 0 . OR . NN . GT .NNOD) THEN WRITE (NTM. 2210) NN WRITE (NOT, 2210) NN

ERR*. TRUE. RETURN

C ERROR TRAP FOR TEMP -PRESCRI BED NODES ELSEIF (ID (NN) .LE. 0) THEN WRITE (NTM. 2220) NN

WRITE (NOT. 2220) NN

ERR*. TRUE.

RETURN

ENDIF

200 CONTINUE

N1 * LM (1

)

N2 * LM (2) C COMPUTE SURFFACE SEGMENT LENGTH

XL * SORT ( (XYZ (Nl,l) -XYZ (N2,l))**2

. (XYZ (Nl. 2) -XYZ (N2.2) )**2)

C COMPUTE CH

IF (REG.EQ. 'PLN' ) THEN

CH (Nl ) * CH (Nl) D*COEF*XL/2.

0

CH(N2) « CH(N2) D*COEF*XL/2 .

0

ELSEIF (REG.EQ. ’AXI ’ ) THEN CH (Nl ) * CH(N1)4 (D*COEF*XL/6. 0)* (2.0*XYZ (Nl. 1) + XYZ (N2. 1 )

)

CH(N2) * CH(N2)4* (D*COEF*XL/6 . 0)* (XYZ (Nl. 1) 2 . 0*XYZ (N2. 1 ) ) ENDIF

C

C —3.0 MODIFY CONDUCTANCE TRANSFER MATRIX. K CC COMPUTE CONVECTIVE CONTRIBUTION TO K IF (REG.EQ. ’PLN' ) THEN

S(l,l) - D*COEF*XL/3 .

0

S (1 . 2) - D*COEF*XL/4.

0

S (2 , 1 ) - S (1 , 2) S (2 , 2) - S(l.l) ELSEIF (REG.EQ. ’AXI’) THEN

S(l.l) - (D*COEF*XL/4.0)* (XYZ(Nl.l) XYZ (N2 . 1 ) /3 . 0) S (1 , 2) - (D*COEF*X1/12.0)*(XYZ(N1,1) XYZ(N2,1))

S (2 , 1 ) - S(1.2)

S (2 , 2) - (D*COEF*XL/ 4.0)* (XYZ (Nl , 1 ) /3 XYZ(N2,1>) ENDIF

C ADD CONVECTIVE CONTRIBUTION TO SYSTEM ARRAYS

CALL ADDCM (ID, LM, S , K, 2, 2, NNOD, NEGN . MBAN)

RETURN

2000 FORMAT (6X,I4,2I5.5X. 2G10. 3, 3X, A3) 2200 FORMAT ( ’ ***• ERROR: Region ' , A3, * is not defined.’)

2210 FORMAT (' •••• ERROR: (Generated) Node’,15, ’ la out of range.’) 2220 FORMAT (’ •*•• ERROR: (Generated) Node’. 15,

.* la not part of flux-prescribed boundary')

SUBROUTINE RLSOLV (MPID. MPK, MPC, MPEG. MPCH . NNOD. NEQN . MBAN . ERR) C--SUB : RLSOLV - ROOT SUBROUTINE TO ALL LINEAR SOLVER SUBROUTINES

C

C DICTIONARY OF VARIABLES

C

C VARIABLE

C

C ERR

C NNOD

C NEQN

C MBAN

C IWRT

C IPRT

C NPRT

C MPNPR

C MPCH

C MPEO

C MPC

C MPK

C MPID

C NPR(NPRT)

C CH (NNOD) EXCITATION

C EO(NNOO)

C C (NEQN . MBAN) C K (NEQN , MBAN)

C ID (NNOD) Cc

c

c

C CAL -SAP : DATA 4 COMMON STORAGE

C COMMON MTOT.NP, IA(1) COMMON /DBSYS/ NUMA, NEXT. IDIR, IP (3) COMMON /IOLIST/ NTM. NTR. NIN , NOT. ND1 . ND2 , ND3 . ND4

CC RFORM: DATA 4 COWON STORAGE

C LOGICAL* 1 ERR, NOSOLV

INTEGER ITIME1. ITIME2

SAVE /DBSYS/, /IOLIST/

CC--0.0 WRITE SOLUTION HEADER 4 WRITE SYSTEM MATRICES TO <f llename> . SYS C WRITE (NTM. 2100) WRITE (NOT. 2100) 2100 FORMAT (/,’ -= SOLUTION’) CALL SYSAVE (IA(MPK) , IA (MPC) .NEQN. MBAN) NOSOLV * .TRUE.

CC—1.0 DEFINE E(NECN)

C

CALL DELETECE ’)

CALL DEFINECE ’

. MPE . NEQN. 1

)

CC—2.0 GET REPORT CONTROL INFORMATION

C CALL REPORT (MPNPR, NPRT. IPRT. IWRT, NNOD. ERR)

IF (ERR) RETURN

CC--3.0 LOOK FOR SEPARATOR "STEADY"

C CALL FINDN (’STEADY’ , 6. KEY)

IF (KEY . EQ. 0) THEN

NOSOLV*. FALSE. CALL ZEROR (IA (MPE) .NEQN, 1)

CALL TIME (ITIME1

)

CALL STEADY (IA (MPID) . IA (MPK) . IA (MPE) , IA (MPEO) .IA(MPCH)

,

IA (MPNPR) .NPRT. IWRT. NNOD, NEQN. MBAN. ERR) CALL TIME (ITIME2)

WRITE (NTM. 2300) (ITIME2-ITIME1

)

WRITE (NOT. 2300) (ITIME2-ITIME1

)

ENDIF

2300 FORMAT (/’ — NOTE: Solution time: ',15, • aeconda.’)

DESCRIPTION

DO-WHILE TERMINATOR FLAG NUMBER OF NODES IN SYSTEM NUMBER OF (SYSTEM) EQUATIONS

(HALF) BANDWIDTH OF SYSTEM EQUATIONS

WRI TE- SOLUTION -TO- FILE FLAG (l*DO IT) OUTPUT PRINT INTERVAL NUMBER OF NOOES SELECTED FOR SOLUTION OUTPUT PRINT POINTER TO NPR (NPRT) IN BLANK COMMON POINTER TO CH(NNOO) IN BLANK COWON

POINTER TO RO (NNOD) IN BLANK COWON

POINTER TO C (NECN . MBAN) IN BLANK COWON

POINTER TO K (NEQN, MBAN) IN BLANK COMMON POINTER TO ID (NNOD) IN BLANK COWON

LIST OF NODES SELECTED FOR SOLUTION OUTPUT PRINT

COEFS. FOR COMPUTING EFFECTIVE CONVECTIVE

INITIAL EXCITATION VECTOR (STORED BY NODE) SYSTEM CAPACITY MATRIX (COMPACT FORM) SYSTEM CONDUCTANCE TRANSFER MATRIX (COMPACT FORM) EQUATION NUMBER /B.C. FLAG ARRAY; ABS (ID (N) ) « EQUATION NUMBER OF NODE N SIGN (ID (N)) - -1 - TEMPERATURE PRESCRIBED NOOE

SIGN (ID (N) ) - ! - FLUX PRESCRIBED NODE

CC—4.0 LOOK FOR SEPARATOR "HARMON"

C CALL riNDNC HARMON’, 6, KEY) IF ( KEY . EQ . 0 ) THEN

NOSOLV*. FALSE. CALL ZEROR (IA (MPE) .NEON, 1)

CALL TIME (ITIME1

)

END

C»* *•*••*••»*•«•••**»•*•••••••**•••*••••*•*•**•***•• •*••***»••*•<

C TREE: LSOLV - TREE OF SUBROUTINES TO SOLVE LINEAR SYSTEM EQTNS

c

CALL HARMON (IA (MPID) . IA (MPK) . IA (MPC) . IA (MPE) . IA (MPEO) . IA (MPCH) .

. IA (MPNPR) . IA (MPK) .IA(MPE) . NPRT. IWRT. NNOD, NEQN . MBAN . ERR) CALL TIME (ITIME2

)

WRITE (NTM. 2300) (ITIME2-ITIME1

)

D -1 2
WRITE (NOT. 2300) (ITIME2-ITIME1) ENDIF C

c —5.0 LOOK FOR SEPARATOR "PREDICT"

C CALL FINDN (’PREDICT*. 7. KEY) IF (KEY . EQ. 0) THEN NOSOLV*. FALSE. CALL ZEROR(IA(MPE) .NECN.l)

CALL TIKE (ITINE1

)

CALL PREDIC (IA (MPID) . IA (MPK) . IA (HPC) . IA (MPE) . IA (KPEO) . IA (MPCH)

,

. IA(KPNPR) .NPRT.IPRT, IWRT. NNOD, NEON . MBAN. ERR) CALL TIKE (ITIME2

)

WRITE (NTH. 2300) (ITIME2-ITIKE1)

WRITE (NOT. 2300) (ITIHE2-ITIKE1) END IF CC—6.0 REPORT NO SOLUTION REQUEST ... C

IF (NOSOLV) THEN WRITE (NTH. 2600) WRITE (NOT. 2600) ENDIF

CALL FCLOSE (ND1

)

RETURN

2600 FORMAT ( / ' — NOTE: No solution control data found.')

END

C SYSAVE SUBROUTINE SYSAVE (K. C. NEQN . MBAN

)

C-SUB : SYSAVE - WRITES (K) 4 (C) TO <f llenarae> . SYS C

C CAL-SAP : DATA 4 COMMON STORAGE

C COMMON /IOLI ST/NTM, NTR . NIN . NOT . ND1 , ND2 , ND3 . ND4 C

C SYSAVE: DATA

C REAL* 6 K (NEQN , MEAN ) . C (NEQN . MBAN

)

C CALL FOPEN (ND1. ’SYS’, ’NEW’

)

WRITE (ND1) ((K(N.M) , N=1 , NEQN) . M=1 , MBAN )

,

+ ((C(N,M),N*1, NEQN ) , M=1 . MBAN)

RETURN

END

c STEADY SUBROUTINE STEADY (ID. K, E, EO. CH, NPR. NPRT, IWRT. NNOD. NEQN , MBAN . ERR) C-SUB: STEADY - FORMS (E) OF (K)(T) * (E) C CALL SOLVER AND REPORTS RESULTS C SOLUTION (T) WRITTEN OVER (E) IMPLICIT REAL* 8 (A-H.O-Z)

C

C CAL -SAP: DATA 4 COMMON STORAGE

C COMMON /IOLI ST/NTM. NTR. NIN , NOT, ND1 , ND2 , ND3 , ND4

C

C STEADY: DATA 4 COMMON STORAGE

C REAL* 8 K (NEQN. MBAN) . E (NEQN) , EO (NNOD) , CH (NNOD) INTEGER ID (NNOD). NPR (NPRT) , IJK(3) LOGICAL *1 ERR. TDOF CHARACTER ENDFLAG* 3

SAVE /IOLIST/

CC—1.0 WRITE HEADER

C WRITE (NTH, 2100) WRITE (NOT, 2100) 2100 FORMAT (/’ STEADY STATE SOLUTION’) CALL FREETY

CC—2.0 SCALE (K) FOR TEMP -PRESCRIBED NODES C

DO 20 I -1. NEQN

20 IF (TDOF (I. ID. NNOD)) K(I,1) - K (I . 1) *1 . OE+15

CC~3.0 GET T, Q. AND/OR TEX 4 FORM { E

)

C WRITE (NTM. 2300) WRITE (NOT. 2300) WRITE (NOT, 2302)

30 CALL FREE CALL FREETY

C CHECK FOR "END" CALL FREEH (’ ’. ENDFLAG. 3. 1

)

IF (ENDFLAG. EQ. ’END’ ) GO TO 40 C INTERPRET LINE AND GENERATE DATA Q - 0.0 CALL FREER ('O’ .Q. 1) T = 0.0 CALL FREER (’T* . T, 1) TEX - 0.0

CALL FREER CX* .TEX. 1) CALL FREEI ( ' ’.IJK(1),3)

IF (I JK (2) .EQ. 0) I JK (2) “I JK (1

)

IF (I JK (3) . EQ. 0) I JK (3) =1 DO 32 N«IJK(1) ,IJK(2).IJK(3)

IF (N . LE . 0 . OR . N . GT . NNOD ) THEN WRITE (NTM. 2310) N WRITE (NOT. 2310) N ERR* . TRUE

.

GO TO 30 ENDIF

NN - ID (N) NNN - ABS (NN) C TEMPERATURE-PRESCRIBED NODES

IF(NN.LT.O) THEN WRITE (NOT. 2320) N.T

E (NNN) - T*K (NNN , 1

)

C FLUX -PRESCRIBED NODES ELSEIF (NN.GT.0) THEN WRITE (NOT. 2330) N. EO (N) . Q. TEX

E (NNN) E (NNN) EO(N) + Q + TEX*CH(N) ENDIF

32 CONTINUE

GO TO 30 C

C —4 . 0 SOLVE

C

40 IF (ERR) RETURN CALL SYMBC(K.E.NECN. MBAN. NEQN. 1) CC—5.0 REPORT SOLUTION (T) C

C WRITE TO <FILENAME> . TMP OPTION

IF (IWRT . EQ. 1) THEN

TIME * 0.0 CALL FOPEN (ND1. ’TMP ’.’NEW’) WRITE (ND1) TIME, (E (N) . N-l . NEQN) CALL FCLOSE (ND1) ENDIF

C PRINTABLE OUTPUT

WRITE (NTM. 2500) WRITE (NTM. 2510) (NPR (N) . E (ABS (ID (NPR (N) ) ) ) . N»1.NPRT) WRITE (NOT. 2500) WRITE (NOT. 2510) (NPR (N) . E (ABS (ID (NPR (N) ) ) ) ,N=1.NPRT)

RETURN

2300 FORMAT (/’ — EXCITATION’)

2302 FORMAT (

/

.6X'Node Terap. Int.Flux Ext. Flux

. Ext . Temp. '

)

2310 FORMAT ( ’ •*•* ERROR: (Generated) Node'. 15,' Is out of range.')

2320 FORMAT ( 6X , I 4 , 9X , G1 0 . 3

)

2330 FORMAT (6X, 14, 19X, 3 (5X.G10.3)

)

2500 FORMAT (/ ’ — RESPONSE'//

.6X,'Node Tenp . Node Terap. Node Terap. Node

Terap. Node Tern?.') 2510 FORMAT ( (6X. 5 (I4.G11 .3) )

)

END

C HARMON SUBROUTINE HARMON (ID. K.C.E. EO. CH. NPR, KSTAR. ESTAR.

.NPRT. IWRT, NNOD. NEQN. MBAN. ERR) C-SUB: HARMON - FORMS (K*)(T*> « ( E*

)

C CALL SOLVER AND REPORTS RESULTS

C COMPLEX WRITTEN OVER DOUBLE PRECISION STORAGE CELLS:

C (K* 1 WRITTEN OVER [K] C (E*) WRITTEN OVER (E) C IMPLICIT REAL* 8 (A-H.O-Z) C

C CAL -SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/NTM.NTR.NIN.NOT.NDl.ND2.’ND3.ND4

c

C HARMON: DATA 4 COMMON STORAGE

C

REAL* 8 K (NEQN. MBAN) , C (NEQN, MBAN) , E (NEQN) , EO (NNOD) ,CH (NWOO) . +0(2) ,T (2) , TEX (2) INTEGER ID (NNOD) , NPR (NPRT) . IJK(3) COMPLEX KSTAR (NEQN. MBAN ) , ESTAR (NEQN) . QSTAR. TSTAR, TXSTAR , 1 LOGI CAL* 1 IRR, TDOF CHARACTER ENDFLAG*

3

REAL ZLAG

SAVE /IOLIST/

DATA PI/3.141592654/

CC—1.0 WRITE HEADER

C WRITE (NTM, 2100) WRITE (NOT. 2100) 2100 FORMAT (/’ STEADY HARMONIC SOLUTION’) CALL FREETY

CC--2.0 GET FREQUENCY

C

D -1 3

CALL FREE

CALL FREETY W - 0.0
2610 FORMAT ( (6X,3 (14, 2G 10.3))

)

CALL FREER ('F' ,W. 1) WRITE (NOT, 2200) W 2200 FORMAT (/ ' Frequency of excltat ion: * ,G10 .3, ' cycles /time '

)

M - W* 2 . 0*PI CC—3.0 GET (K) 4 (C), FORM (K*), 6 SCALE (K* ) FOR TEMP -PRESCRIBED

NODES

C REWIND ND1 READ (ND1 ) < (K(N,M) .N-l. NEQN) ,M-1, MBAN) ,

. (<C(N,M) .N-l.NECN) , M=1 , MBAN

)

DO 30 I-l.NEQN DO 30 J-l.MBAN

30 KSTAR (I , J) - CMP LX (K (I , J) , W'C(I,J|)

DO 32 I-l.NEQN

32 IF (TDOF (I, ID.NNOD) ) KSTAR(I.l) - KSTAR (1,1)* (1 . 0E+ 1 5 , 0.0) C

C —4.0 GET T, Q. AND/OR TEX 4 FORM (E) C WRITE (NTM, 2400) WRITE (NOT, 2400) WRITE (NOT, 2402) 40 CALL FREE CALL FREETY

C CHECK FOR "END” CALL FREEH (' ', ENDFLAG . 3. 1

)

IF (ENDFLAG . EQ. ' END ' ) GO TO 50 C INTERPRET LINE AND GENERATE DATA

C CONVERT: Z - MAG* COS (-LAG*W) * 1 MAG*SIN (-LAG*W) Q (1 ) - 0.0 Q(2) - 0.0 CALL FREER <’Q’ ,Q. 2) QSTAR - CMPLX (Q(l) *COS (-Q(2) *W) ,Q(1) *SIN (-0(2) *W)

)

T (1 ) - 0.0

T (2 ) - 0.0 CALL FREER ( 'T* , T, 2) TSTAR = CMPLX (T (1) *COS (-T (2) *W) ,T(1) *SIN (-T (2) *W)

)

TEX(l) =0.0

TEX (2) - 0.0 CALL FREER (’X’, TEX. 2) TXSTAR = CMPLX (TEX (1) *COS (-TEX(2)*W), TEX (1 ) *SIN (-TEX (2) *W)

)

CALL FREEI (

* ',IJK(1),3)

IF (IJK (2) .EQ. 0) IJK (2) =IJK (1) IF (I JK (3) . EQ. 0) IJK (3) =1 DO 42 N=IJK (1) , IJK (2) , IJK (3) IF (N.LE.O.OR.N.GT.NNOD) THEN WRITE (NTM. 2410) N WRITE (NOT, 2410) N ERR= . TRUE

.

GO TO 40 ENDIF

NN - ID (N) NNN = ABS (NN) C TEMPERATURE-PRESCRIBED NODES

IF(NN.LT.O) THEN WRITE (NOT. 2420) N,T(1),T(2) ESTAR (NNN ) - TSTAR* KSTAR (NNN . 1

)

C FLUX-PRESCRIBED NODES ELSEIF (NN.GT. 0) THEN WRITE (NOT. 2430) N. EO (N) , Q (1) . Q (2) .TEX(l) ,TEX(2) ESTAR (NNN) - ESTAR (NNN) CMPLX (EO(N)) QSTAR CH(N)*TXSTAR

ENDIF

42 CONTINUE GO TO 40 CC—5.0 SOLVE

C

50 IF (ERR) RETURN CALL SYMBCX (KSTAR. ESTAR, NEQN. MBAN. NEQN. 1

)

CC —6.0 REPORT SOLUTION { T ) C

C WRITE TO <FILENAME> . TMP OPTION

IF (IWRT.EQ.l) THEN TIME - 0.0 CALL FOPEN(NDl, 'TMP ’, , NEW’) WRITE (ND1) TIME. (ESTAR (N) , N=1 . NEQN) CALL FCLOSE(NDl) ENDIF

C PRINTABLE OUTPUT

WRITE (NOT. 2600) WRITE (NOT. 2610) (NPR (N) , ABS (ESTAR (ABS (ID (NPR (N) ) ) ) ) ,

. SLAG (ESTAR (ABS (ID (NPR (N) ) )) , W) . N-l.NPRT) WRITE (NTM, 2600) WRITE (NTM, 2610) (NPR (N) , ABS (ESTAR (ABS (ID (NPR (N) )))) ,

. ZLAG (ESTAR (ABS (ID (NPR (N) ) ) ) . W) , N=1 . NPRT)

RETURN

2400 FORMAT (/’ — EXCITATION’)

2402 FORMAT (/ .23X, f Tenperature Ext. Flux Ext. Tenp*/

. 6X, 'Node Int.Flux Arap. Lag Artp. Lag

Arap . Lag '

)

2410 FORMAT ( ' **•* ERROR: (Generated) Node', 15. * la out of range.')

2420 FORMAT (6X, 14, 10X, 2G10 . 3) 2430 FORMAT ( 6X . I 4 . G1 0 . 3 , 20X , 4G1 0.3) 2600 FORMAT (/ ' — RESPONSE’//

.13X’ Tenperature Temperature

Tenperature’

/

. 6X’Node Anp. Lag Node Anp. Lag

. Node Anp. Lag')

END

C ZLAG FUNCTION ZLAG(Z.W) C-FUN : ZLAG - COMPUTES LAG-ARG/FREQ OF COMPLEX NUMBER COMPLEX Z REAL* 8 W REAL ZLAG, PI, IM, RE DATA PI/3.141592654/

IM - AIMAG(Z)

RE - REAL (Z)

IF( (ABS (RE) .LT.1.0E-6) .OR. (W . LT . 1 . 0D-6) ) THEN ZLAG - 0.0 ELSE

ZLAG - -ATAN (IM/RE) IF (IM . GT. 0.0. AND . RE .GT .0.0) ZLAG - 2.0*PI+ZLAG

IF (RE . LT .0.0) ZLAG - PI+ZLAG

ZLAG - ZLAG/REAL (W) ENDIF

RETURN

END

C

PREDICSUBROUTINE PREDIC (ID, K. C. E, EO.CH. NPR. NPRT. IPRT, IWRT. NNOD, NEQN, MBAN, ERR) C-SUB: PREDIC - PREDICTOR -CORRECTOR 1ST O.D.E. EQUATION SOLVER

C BASED ON SOLVER IN *HEAT* BY R.L. TAYLOR - U.C. BERKELEY

C SOLVES EQUATION;

(K)(T) (C) (dT/dt ) - ( E (t )

)

FOR GENERAL EXCITATION, { E (t ) > . DEFINED BY A SET OF PIECEWISE LINEAR FUNCTIONS. METHOD BASED ON DIFFERENCE APPROXIMATION;

(T) r+l - (T) n (1-a) DT{dT/dt ) n (a) DT(dT/dt }n+l

"alpha", an Integration parameter

0 corresponds to Forward Difference method

1 corresponds to Backward Difference method

1/2 corresponds to Crank-Nichol son method

(unstable)

C DT = time step increment

CC DICTIONARY OF VARIABLES

C INTERNAL TO SUBROUTINE

C

C VARIABLE DESCRIPTION

C NEFN NUMBER OF EXCITATION FUNCTIONS DEFINED

C MPIEFN

IEFN (NNOD, 2)

: EXCITATION FUNCTION NUMBER ARRAY

C MPEFN

EFN (NEFN. 2)

: EXCITATION FUNCTION DATA

C E (NEQN)

: CURRENT (E) (ORDERED BY EQTN #) C MPT

T (NEQN)

: CURRENT (T) (ORDERED BY EQTN #) C MPTD

TD (NEQN)

: CURRENT (dT/dt) (ORDERED BY EQTN #) C MPTDD TDD (NEQN)

: INITIAL (d/dt (dT/dt ) ) TO EST TIME STEP

C IEFN (NNOD, 2) LIST OF FUNCTION NUMBERS (ORDERED BY NODE NUMBER) C IEFN (I. 1) = TEMP OR EXT. FLUX FUNCTION NUMBER

C IEFN (I, 2) - CONVECTIVE FLUID TEMP. FUNCTION NUMBER

C EFN (NEFN. 2) EXCITATION FUNCTION DATA

C EFN (1,1) « OLD VALUE

C EFN (I, 2) - NEW VALUE

C TOLD. TNEW EXCITATION FUNCTION DATA TIMES C TIME (1) START TIME

C TIME (2) END TIME C TIME (3) TIME INCREMENT

IMPLICIT REAL*8 (A-H.O-Z) C

C CAL-SAP: DATA 4 COMMON STORAGE

C COMMON MTOT.NP, IA(1) COMMON /DBSYS/ NUKA. NEXT. IDIR, IP (3) COMMON /IOLIST/ NTM. NTR . NIN . NOT . ND1 . ND2 . ND3 . ND4

CC PREDIC: DATA 4 COWON STORAGE

C REAL* 8 K (NEQN, MBAN) . C (NEQN, MBAN) . E (NEQN) , EO (NNOD) , CH (NNOD) INTEGER ID (NNOD), NPR (NPRT) CHARACTER ENDFLAG*

3

LOGICAL* 1 ERR

COMMON /LIST/ TIME (3) .ALPHA. TOLD. TNEW

SAVE /DBSYS/, /IOLIST/. /LIST/ D -1 4
C WRITE HEADER CALL FREETY

C WRITE (NOT. 2000) 2000 FORMAT (/* — DYNAMIC SOLUTION: PREDICTOR-CORRECTOR

INTEGRATION’

)

C

C —1.0 GET SOLUTION CONTROL INFORMATION

C CALL FREE CALL FREETY CALL FREER ('E'. TIME (1 ) .3) ALPHA « 0.75 CALL FREER (’A* .ALPHA. 1) IF(ALPHA. LT. 0.0. OR. ALPHA. GT. 1.0) THEN WRITE (NTH. 2100)

WRITE (NOT. 2100)

ERR-. TRUE.

RETURN

2100 FORMAT (’ ERROR: Alpha «jst be In range 0.0 to 1.0.') ENDIF NEFN - 0 CALL FREEI (’N* .NEFN.l) IF (NEFN . LE . 0) THEN WRITE (NTM.2110)

WRITE (NOT. 2110)

ERR- . TRUE

.

RETURN

2110 FORMAT ( ' **** ERROR: One or more excitation functions

. must be defined.') ENDIF

C REPORT CONTROL INFORMATION WRITE (NOT, 2120) TIME. ALPHA . NEFN

2120 FORMAT (/’ — SOLUTION CONTROL INFORMATION’/

.’ Start time ’.G10.3/

End time ' .G10.3/

.' Time step Increment '.G10.3/

.

’ Integration parameter, alpha '.G10.3/

. ' Number of excitation functions ... ’,16) CC—2.0 DEFINE ARRAYS

C CALL DELETE ('IEFN') CALL DELETE (

’ EFN ') CALL DELETE (’T ’) CALL DELETE (

' TD ’) CALL DELETE (’TDD ’) CALL DEFINE ('TDD '

, MPTDD, NEQN , 1

)

CALL DEFINE (’ TD * , MPTD, NEQN, 1

)

CALL DEFINE ('T ’. MPT. NEQN . 1

)

CALL DEFINE (’EFN ’ . MPEFN, NEFN , 2

)

CALL DEFINI (’IEFN’ . MPIEFN. NNOD, 2) NS TOR « (IDIR-NEXT-20)

IF (NSTOR.LT. 0) THEN WRITE (NTM. 2200) -NSTOR

WRITE (NOT. 2200) -NSTOR

ERR= . TRUE

.

RETURN

2200 FORMAT (' •*** ERROR: Insufficient storage to form arrays for

. solution.'/

. ' Increase IA (MTOT) by ’,16) ENDIF

C INITIALIZE CALL ZEROI (IA (MPIEFN) .NNOD. 2) CALL ZEROR (IA (MPEFN) .NEFN. 2) CALL ZEROR (IA (MPT) .NEQN. 1) CALL ZEROR (IA (MPTD) .NEQN. 1) CALL ZEROR (IA (MPTDD) , NEQN , 1

)

CC—3.0 GET EXCITATION CONTROL INFORMATION

C CALL ECNTRL (ID. IA (MPIEFN) . NNGD , ERR) IF (ERR) RETURN

CC—4.0 GET INITIAL CONDITIONS

C

C 4.1 WRITE HEADER

C WRITE (NTM, 2400) WRITE (NOT, 2400) 2400 FORMAT (/’ — INITIAL CONDITIONS')

C

C 4.2 FIND SEPARATOR "INITIAL"

C CALL FINDN (’INITIAL '.7, KEY) IF (KEY . EQ. 1) THEN WRITE (NTM. 2410)

WRITE (NOT. 2410) 2410 FORMAT (/’ — NOTE: Initial condition data not found.'/

.’ Initial temperatures assumed to be zero.’) GO TO 50 ENDIF

C 4.3 CALL ICOND TO DO THE WORK

C CALL ICOND (ID, IA (MPT) .NNOD. NEQN, ERR) CC—5.0 4 READ (K) 4 (C) FROM <f llenarae> . SYS

C 50 REWIND ND1 READ (ND1 ) ( (K (N, M) . N— 1, NEQN) . M-l , MBAN) ,

. ((C(N.M) ,N-1 .NEQN) .M-l.MBAN)

CC--6.0 INITIALIZE EXCITATION FUNCTION VALUES

C CALL FINDN ('EXCIT' . 5, KEY)

IF (KEY . EQ. 1) THEN WRITE (NTM. 2400) WRITE (NOT. 2600) 2600 FORMAT (' •••* ERROR: Separator "EXCIT" not found.’) ERR- . TRUE

.

RETURN ENDIF

DO 60 1-1.2 CALL GETEFN (IA (MPEFN) .TOLD. TNEW. NEFN. ENDFLAG)

IF (ENDFLAG.EQ. ’END’ ) THEN WRITE (NTM. 2610)

WRITE (NOT. 2610)

RETURN

60 ENDIF

2610 FORMAT (’ *•*• ERROR: Insufficient excitation data; at least two values must be given for excitation functions.’)

CC--7.0 CALL PREDI TO DO THE WORK

C CALL PREDI (ID, K, C, EO, IA (MPT) .CH.NPR. IA (MPIEFN) . I A (MPEFN) . E.

. IA (MPTD) . IA (MPTDD) , NPRT, IPRT, IWRT, NEFN. NNOD. NEQN . MBAN . ERR)

RETURN

END

CICOND SUBROUTINE ICOND (ID. T, NNOD. NEQN . ERR) C-SUB: ICOND - READS AND GENERATES INITIAL CONDITIONS

C IMPLICIT REAL* 8 (A-H.O-Z)

C

C CAL -SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/ NTM. NTR. NIN. NOT. ND1 . ND2 . ND3 . ND4 COMMON /PARC/ FIN (12) . EXT (3) CC ICOND: DATA 4 COMMON STORAGE

C REAL *8 T (NEQN) INTEGER ID (NNOD). IJX(3) CHARACTER ENDFLAG* 3. READ* 5 , FN(12)*1

LOGICAL* 1 ERR

SAVE /IOLIST/, /PARC/

CC--1.0 CHECK TO SEE IF READ -FROM- D I SK OPTION IS REQUESTED

C CALL FREE

CALL FREEH ( * ’, READ, 5,1) IF (READ. EQ. ’READ-’ ) THEN CALL FREETY

C GET FILENAME AND CREATE EXTENSION

C — STORE CURRENT FILENAME TEMPORARILY DO 10 1-1,12

10 FN (I ) - FIN (I) CALL FREEH ('O’ .FIN, 1. 12) C OPEN FILE <f 1 lenarae) . TMP , MOVE TO EOF. BACKSPACE ONE RECORD CALL FOPEN (ND1. ’TMP ’.’OLD’)

12 READ (ND1 , END-14) GO TO 12

14 BACKSPACE ND1

C READ INITIAL CONDITIONS FROM LAST RECORD READ (ND1) TIME, (T(N) ,N-1, NEQN) CALL FCLOSE(NDl) WRITE (NOT, 2100) (FIN (N) ,N«1,12) , (EXT (N) ,N«1,3)

WRITE (NOT. 2200) WRITE (NOT. 2400) (N . T (ABS (ID (N) ) ) . N-l . NNOD) C — RESTORE CURRENT FILENAME DO 16 1-1,12

16 FIN (I) - FN (I

)

RETURN

ENDIF

CC—2.0 WRITE INITIAL CONDITIONS HEADER

C WRITE (NOT, 2200)

CC—3.0 INTERPRET LINE OF DATA

C 30 CALL FREETY

C CHECK FOR "END- CALL FREEH (

’ ' .ENDFLAG. 3.1) IF (ENDFLAG.EQ. ’END' ) GO TO 40 D -1 5
C DO IT TEMP - 0.0 CALL FREER ( 'T ’. TEMP . 1

)

CALLFREEIC MJK(1) # 3) IF (I JK (2) . EQ. 0) IJK(2)-IJK(1)

IF (I JK (3) . EQ. 0) IJK(3)=1

DO 34 N-IJK(l) ,IJK(2) ,IJK(3) IF(N.LE.O.OR.N.GT.NNOD) THEN

WRITE (NTM. 2300) N

WRITE (NOT. 2300) N

ERR-. TRUE.

GO TO 36END IF

34 T (ABS (ID (N) ) ) - TEMP

C GET NEXT LINE OF DATA

36 CALL TREE GO TO 30 CC—4.0 REPORT INITIAL CONDITIONS IF NO ERROR ENCOUNTERED

C

40 IF (ERR) RETURN WRITE (NOT, 2400) (N. T (ABS (ID (N) ) ) ,N=l,NNOD)

RETURN

2100 FORMAT (/' — NOTE: Initial conditions read from last record

of: '

. 1 6A1

)

2200 FORMAT (/

. 6X, 'Node Temp. Node Temp

.

Node Temp . Node

. Temp

2300 FORMAT (' Node Temp . '

)

••*• ERROR: (Generated) Node '

. 15. ' Is out of range

2400 FORMAT ( (6X, 5 (14 ,G1 1 . 3) ) )

END

c PREDI SUBROUTINE PREDI (ID, K. C, EO. T, CH . NPR . IEFN . EFN . E. TD. TDD, NPRT, IPRT. IWRT.NEFN.NNOD.NEQN.MBAN, ERR) C-SUB: PREDI - THE KERNEL OF PREDIC

IMPLICIT REAL* 8 (A-H.O-Z)

C

C CAL-SAP: DATA 4 COMMON STORAGE C

COMMON /IOLIST/NTM, NTR, NIN. NOT. ND1 , ND2 . ND3 . ND4

C

C L1STP : DATA 4 COMMON STORAGE

C REAL* 8 K (NEON, MBAN ) , C (NEQN . MBAN ) .EO(NNOD) , T (NEQN ) , CH (NNOD) . .EFN (NEFN.2) ,E(NEQN) . TD (NEQN) .TDD (NEQN) INTEGER ID (NNOD), NPR (NPRT) , IEFN (NNOD. 2) LOGICAL*l ERR. SKIP, TDOF

COMMON /LIST/ TIME (3 ), ALPHA. TOLD, TNEW

SAVE /IOLIST/, /LIST/

WRITE (NOT. 2000) WRITE (NTM. 2000) 2000 FORMAT ( / ' — TIME STEP*)

CC—1.0 COMPUTE INITIAL TEMPERATURE RATES: (dT(0)/dt) C

C 1.1 FOR TEMP -DOF: SCALE (C]=[C)*1E15 OR SET (C]*1E15 FOR

C MASSLESS NODES

C SKIP = .FALSE. DO 10 N-l.NEQN

IF (TDOF (N. ID.NNOO) ) THEN

C (N , 1 ) - C (N, 1) *1 . 0E1S

IF (C (N, 1 ) . EQ .0.0) C (N , 1 ) - 1.0E15 ENDIF

10 CONTINUE

C

C CHECK FOR STEP CHANGE IN SPECIFIED INITIAL TEMP

C DO 12 11*1, NNOD

NN - ID (N) IFOm.LT.O) THEN

WNN * ABS (NN)

NFNT - IEFN(N.l)

TEMP - 0.0

I F (NFNT. WE. O) TEMP EFN (NFNT. 1) IF(T(NNN) .NE.TEMP) THEN

WRITE (NTM. 2100)

WRITE (NOT, 2100)

SKIP - .TRUE. ENDIF

ENDIF

12 CONTINUE

2100 FORMAT ( / ' — NOTE: Unable to estimate time step for initial step change In specified nodal temperature (s) .

'

)

C

C 1.2 FORM (E(t-0)

)

C

TM - 0.0 CALL EXCIT (ID, C, T, TD, E. CH, EO, IEFN . EFN . TM. NPRT. IPRT. IWRT, +NEFN . NNOO. NEQN . MBAN . ERR)

IF (ERR) RETURN

Cc 1.3 FORM RHS ( E) -

(

K) (T) FOR FLUX-DOF, ( dT/dt ) *DIAG [C) FOR TEMP -DOF

C CALL RHS (ID. T,TD. E. K.C.NNOO, NEQN, MBAN)

CC 1.4 SOLVE (C) (dT/dt) *= (E) C CALL SYMBC (C, TD, NEQN , MBAN , NEQN, 1

)

IF (SKIP) GO TO 30

CC—2.0 COMPUTE TIMESTEP CHECK

C

C 2.1 COMPUTE INITIAL RATE OF TEMP RATES: d/dt (dT (0) /dt

)

C CALL RHS (ID. TD, TDD, TDD. K.C. NNOD, NEQN. MBAN) CALL SYMBC (C. TDD. NEQN. MBAN. NEQN, 2) CC—2.2 COMPUTE NORMS: ||(T(0))||. | | (dT (0) /dt ) | | , | I d/dt (dT (0) /dt ) I I C

TN - 0.0

TDN - 0.0 TDDN =0.0

DO 22 N-l.NEQN

TN - TN T (N) **2

TDN * TDN TD (N) **2

22 TDDN « TDDN TDD(N)**2

TN - SQRT(TN)

TDN - SQRT (TDN) TDDN - SQRT (TDDN) C

C 2.3 EVALUATE TAYLORS EXPRESSION FOR TIME STEP ESTIMATE

C

B - 0.05

IF (TDDN. NE. 0.0) THEN DTEST = (B*TDN + SQRT (B*B*TDN*TDN 2 . 0*B*TN*TDDN) ) /TDDN WRITE (NTM, 2200) B*100.0. DTEST, TIME (3

)

WRITE (NOT, 2200) B*100.0. DTEST. TIME (3) 2200 FORMAT (/' — NOTE: Estimated time step to limit error to

. approx. ' ,F5. 2. '% is:'.G10.3,/

.

* Specified time step ls:',G10.3) ELSE

WRITE (NTM. 2210)

WRITE (NOT, 2210) 2210 FORMAT (/' — NOTE: Unable to estimate time step to limit

. error for the given system.') ENDIF

CC--3.0 READ (C) 4 (K) FROM DISK

C 30 REWIND ND1 READ (ND1 ) ( (K (N , M) . N=1 , NEQN) , M=1 , MBAN)

.

. ((C(N,M) ,N=1 .NEQN) , M=1 , MBAN)

CC--4.0 FORM LHS: ((C) aDT(K)] SCALING TEMP -DOF DIAG TERMS

C ADT = ALPHA* TIME (3) DTA - TIME (3) - ADT DO 40 N-l, NEQN

DO 40 M-l.MBAN

40 C(N.M) * C(N.M) ADT*K (N, M) DO 42 N-l.NEQN

42 IF (TDOF (N. ID.NNOO)) C(N.l) * C (N, 1) *1 . 0E15

CC—5.0 TIME STEP THRU SOLUTION

C WRITE (NOT. 2500) WRITE (NTM, 2500) 2500 FORMAT (/' — RESPONSE')

IOP - 1 ISTEP - 0 DO 500 TM*TIME (1 ) +TIME (3) ,TIME(2) .TIME (3) ISTEP - ISTEP 1 CC 5.1 FORM (E) C

CALL EXCIT (ID, C, T, TD, E. CH, EO, IEFN . EFN . TM, NPRT. IPRT, IWRT. NEFN . NNOD . NEQN . MBAN . ERR

)

IF (ERR) RETURN

C

C 5.1 PARTIAL UPDATE OF T: (T) - (T) (1-a) DT( dT/dt

)

C

DO 51 N=1 . NEQN

51 T (N) - T(N) DTA*TD (N) CC 5.3 FORM RHS: (E)-(K) (T) FOR FLUX-DOF, (dT/dt ) *DIAG (C) FOR TEMP-DOF C

CALL RHS (ID. T.TD.E. K.C. NNOD. NEQN, MBAN) C

C 5.4 SOLVE FOR (dT/dt)

C CALL SYMBC (C. TD. NEQN. MBAN, NEQN. IOP) IOP * 2 D -1 6
cC 5.5 COMPLETE UPDATE OF T: (T) * (T) aDT{dT/dt)

C DO 55 N*1.NEQN

55 T(N) - T (N) + ADT*TD(N)

C

C 5.6 REPORT RESULTS

C

IF (MOD (ISTEP, IPRT) . EQ. 0) THEN

C WRITE TO <FILENAME> . TMP OPTION

IF(IWRT.EQ.l) THEN CALL FOPEN (ND1. ’TMP * .’NEW’) WRITE (ND1 ) TM, (T (N) ,N=1 , NEQN) CALL FCLOSE (ND1

)

END IF C PRINTABLE OUTPUT WRITE (NOT, 2510) TM

WRITE (NOT, 2520) WRITE (NOT, 2530) (NPR (N) . T (ABS (ID (NPR (N) ) ) ) ,N*1 ,NPRT) WRITE (NTM, 2510) TM

WRITE (NTM. 2520) WRITE (NTM. 2530) (NPR (N) , T (ABS (ID (NPR (N) ) ) ) ,N*=1 ,NPRT) END IF

500 CONTINUE

RETURN

2510 FORMAT (/ ' Time: '.Gil. 3) 2520 FORMAT

(

.6X,'Node Ten*>. Node Temp . Node Ten*). Node

. Ten*>. Node Ten*>.’) 2530 FORMAT (( 6X, 5(14, G1 1.3))) END

c ECNTRL SUBROUTINE ECNTRL (ID. IEFN. NNOO, ERR) C-SUB: ECNTRL - FORMS EXCITATION FUNCTION NUMBER ARRAY IEFN (NNOO. 2) CC CAL-SAP : DATA 4 COMMON STORAGE

C COMMON /IOLIST/NTM. NTR . NIN, NOT. ND1 . ND2 , ND3 . ND4

CC ECNTRL: DATA 4 COMMON STORAGE

C INTEGER ID(NNOD) . IEFN (NNOO. 2) . IJK(3) LOGICAL *1 ERR CHARACTER ENDFLAG *3

SAVE /IOLIST/

CHARACTER ENDFLAG* 3

CC—1.0 UPDATE OLD VALUES

C TOLD * THEN

DO 10 I -1, NEFN

10 EFN(I.l) - EFN (1,2) CC--2.0 READ NEW VALUES

C CALL FREE

CALL FREETY

C CHECK FOR "END" CALL FREEH (' ’. ENDFLAG . 3 . 1

)

IF (ENDFLAG. EQ. 'END' ) RETURN CALL FREER CT’ .TNEW, 1) CALL FREER ( 'N ' , EFN (1,2) , NEFN)

RETURN

END

C RHS SUBROUTINE RHS (ID. T. TD. E. K. C. NNOD. NEQN. MBAN) C-SUB : RHS - FORMS RHS OF (C](dT/dt) - (E*) - (E) - (K]{T);

C

C {EMt)) * [E(t) ) - [K) (T(t) ) FOR FLUX-DOF

C (EMt)) - (dT(t) /dt)*DIAG OF [Cl FOR TEMP-DOF

C

CC <E*) IS WRITTEN OVER (TD)

IMPLICIT REAL* 8 (A-H.O-Z)

REAL* 8 T (NEQN) . TD (NEQN) , E (NEQN) , K (NEQN . MBAN ) , C (NEQN . MBAN ) INTEGER ID (NNOD) LOG I CAL *1 TDOF

DO 20 N=l, NEQN

C SCALE BY DIAGONAL FOR TEMP PRESCRIBED NODES

IF (TDOF (N. ID. NNOD)) THEN

TD (N) * TD (N) *C(N, 1) C FORM [E) -IK) (T ) WHERE (K) IS IN COMPACT STORAGE ELSE

TEMP « E(N) - K (N, 1) *T (N)

KK*=N LL=N

DO 10 1*2 . MBAN

WRITE (NOT, 2000) 2000 FORMAT (/* — EXCITATION FUNCTION NUMBERS’//

.6X'Node Temp Ext. Flux Ext . Temp .'

)

10 CALL FREE CALL FREETY

C CHECK FOR "END- CALL FREEH

(

1 ’ .ENDFLAG. 3. 1) IF (ENDFLAG. EQ. ’END’ ) RETURN

C INTERPRET LINE AND GENERATE DATA

IQ * 0 CALL FREEI CQ'.XQ.l)

IT « 0 CALL FREEI ( 'T' , IT, 1) ITEX = 0 CALL FREEI ('X’ , ITEX. 1) CALL FREEI (' * ,IJK(1),3)

IF (I JK (2) . EQ. 0) I JK (2) *IJK (1) IF (I JK (3) . EQ. 0) I JK (3) *1 DO 12 N=I JK (1 ) , I JK (2),IJK(3)

IF (N.LE.O.OR.N.GT.NNOD) THEN WRITE (NTM, 2010) N WRITE (NOT. 2010) N ERR*. TRUE

GO TO 10 2010 FORMAT ( ' **** ERROR: (Generated) Node’, 15,' la out of range.') ENDIF

NN - ID (N) NNN - ABS (NN) C TEMPERATURE-PRESCRIBED NOOES

IF(NN.LT.O) THEN WRITE (NOT, 2020) N.IT

KK * KK-1

IF(KK.GT.O) TEMP = TEMP - K (KK , I

) *T (KK)

LL * LL+1

IF (LL . LE .NEQN) TEMP « TEMP - K(N,I)*T(LL)

10 CONTINUE

TD (N) * TEMP ENDIF

20 CONTINUE

RETURN

END

C EXCIT SUBROUTINE EXCIT (ID.C.T, TD, E.CH.EO, IEFN . EFN . TM. NPRT . IPRT. IWRT, NEFN , NNOD, NEQN . MBAN . ERR

)

C--SUB : EXCIT - FOR t-TIME;

C FOR FLUX-DOF COMPUTES (E(t));

C

C (E (t ) > - (EO) (Q(t) ) (CH)T(TEX(t) ) C

C FOR TEMP-DOF COMPUTES (dT/dt(t)}

C

C (dT(t)/dt) - (l/Dt)(T(t) - T(t-Dt))

C

C FROM EXCITATION FUNCTION DATA STORED AND READ INTO EFN (NEFN. 2) C AND LIST OF FUNCTION NUMBERS

IMPLICIT REAL* 8 (A-H.O-Z)

IEFN(N.l) - IT 2020 FORMAT (CX, 14. 10X, 15) C FLUX-PRESCRIBED NOOES ELSEIF(NN.GT.O) THEN WRITE (NOT. 2030) M, IQ. ITEX

IEFN (N, 1 ) - IQ

IEFN (N, 2) - ITEX

2030 FORMAT ( 6X , 1 4 , 2 OX, 1 5 , 1 OX , 1 5

)

ENDIF

12 CONTINUE

GO TO 10

C

c

c

c

c

c

CAL -SAP : DATA 4 COMMON STORAGE

COMMON /IOLIST/NTM, NTR, NIN. NOT, ND1, ND2, ND3, ND4

EXCIT: DATA 6 COWK3N STORAGE

REAL* 0 C (NEQN , MBAN) . T (NEQN) , TD (NEQN) . CH (NNOD) , EO (NNOD) , E (NEQN)

,

EFN (NEFN, 2) INTEGER ID (NNOO). IEFN (NNOO. 2) LOGICAL* 1 ERR CHARACTER ENDFLAG*

3

END

C GETEFN SUBROUTINE GETEFN (EFN. TOLD, TNEW, NEFN, ENDFLAG) C-SUB :GETEFN - READ INPUT FILE TO UPDATE EXCITATION FUNCTION DATA

IMPLICIT REAL* 8 (A-H.O-Z)

COMMON /LIST/ TIME (3) .ALPHA, TOLD, TNEW

SAVE /IOLIST/, /LIST/

CC--1.0 UPDATE EXCITATION FUNCTION DATA IF NEEDED

C

10 IF (TM.GT.TNEW) THEN D -1 7

REAL* 8 EFN (NEFN. 2)
CALL GETEFN (EFN, TOLD, TNEW, NEFN, ENDFLAG) C--2.0 GET PRINT INTERVAL. WRITE OPTION

C IF (ENDFLAG. EQ. 'END') THEN

WRITE (NTH. 2100)

WRITE (WOT, 2100) 2100 FORMAT ( ' •••• ERROR: Insufficient excitation data.')

ERR - .TRUE. RETURN

ENDIF

GO TO 10ELSEIF (TM.LT.TOLD) THEN WRITE (NTH. 2110) TM

WRITE (NOT, 2110) TM

2110 FORMAT (• *•** ERROR: Excitation data not defined for current

. time: * .G10.3) ERR-. TRUE.

RETURN

ENDIF

C INTERPOLATION FRACTION XT (TM-TOLD) / (TNEW-TOLD) CC—2.0 FORM (E) OR (dT/dt) C CALL SEROR (E. NEQN . 1

)

DO 20 N-l.NNOD

NN « ID (N

)

NNN - ABS (NN) C FLUX-DOF: SUM INT. FLUX + EXT. FLUX + EFFECT. CONVECTIVE FLUX

IF(NN.GT.O) THEN NFNQ = I EFN (NNN. 1)

NFNC * I EFN (NNN. 2)

Q - 0.0

IF (NFNQ.NE. 0) Q - XT* (EFN (NFNQ, 2) -EFN (NFNQ, 1 ) ) 4- EFN (NFNQ. 1)

CONVE = 0.0

IF (NFNC.NE. 0) CONVE - CH (N) • (XT* (EFN (NFNC . 2) - EFN (NFNC. 1))

4- EFN (NFNC. 1)) E (NNN) - E(NNN) + EO(N) Q + CONVE

C TEMP -DOF: COMPUTE (dT/dt) FROM PRESCRIBED TEMP ELSEIF (NN.LT. 0) TEEN MFNT * IEFN(NNN.l)

TEMP -0.0

IF(NFNT.NE.O) TEMP «= XT* (EFN (NFNT, 2) -EFN (NFNT, 1) ) 4- EFN(NFNT.l)

TD (NNN) - (TEMP - T (NNN) ) /TIME (3) ENDIF

20 CONTINUE

RETURN

END

C REPORT SUBROUTINE REPORT (MPNPR, NPRT, IPRT, IWRT. NNOD, ERR) C-SUB: REPORT - READS REPORT CONTROL INFORMATION

C FORMS REPORT CONTROL ARRAY NPR (NPRT) C

e CAL-SAP : DATA & COMMON STORAGE

C CHARACTER EXT*1, FIN*1 COMMON MTOT.NP, IA(1) COMMON /DBSYS/ NUMA.NEXT, IDIR, IP (3) COMMON /IOLIST/ NTM, NTR, NIN . NOT, ND1 . ND2 , ND3 , WD4 COMMON /PARC/ FIN (12) , EXT (80) SAVE /DBSYS/. /IOLIST/, /PARC/ C

C- REPORT: DATA 4 COMMON STORAGE C INTEGER I JK (3) LOGICAL* 1 ERR CHARACTER ENDFLAG* 3, WRITEFLAG*5

SAVE /IOLIST/ CC—0.0 WRITE EEAfiER C WRITE (NTM, 2000) 2000 FORMAT (/’ — REPORT CONTROL’)

CC—1.0 FIND "REPORT"

C CALL FINDN ( ’REPORT’ , 6, KEY) IF (KEY . EQ. 1) THEN WRITE (NTM. 2100) WRITE (NOT. 2100) ERR= . TRUE

.

RETURN

ENDIF CALL FREETY

C

CALL FREE

CALL FREETY CALL FREEH (’ * , WRITEFLAG. 5 , 1

)

IWRT—

0

IF (WRITEFLAG. EQ. ’WRITE’ ) THEN

IWRT - 1 WRITE (NTM, 2200) WRITE (NOT. 2200) ENDIF

IPRT-1 CALL FREEI (

’ T’ , IPRT, 1

)

CC--3.0 FORM PRINT CONTROL ARRAY C

C DEFINE ARRAY NPR (1.1) TO GET POINTER 4 CREATE DIRECTORY ENTRY CALL DELETE ( ’NPR ’) CALL DEFINE (’NPR ’.MPNPR.1,1)

C GENERATE PRINT CONTROL ARRAY WHILE BLANK COMMON SPACE AVAILABLE NPRT—

0

30 CALL FREE

CALL FREETY

C CHECK FOR "END" CALL FREEH (’ ’, ENDFLAG, 3 . 1

)

IF (ENDFLAG. EQ. ’END’ ) GO TO 36 C DO IT CALLFREEK’ •,XJK<1),3)

IF (I JK (2) .EQ. 0) I JK (2) -I JK (1

)

IF (I JK (3) . EQ. 0) I JK (3) -1 DO 34 N-IJK (1) , I JK (2) , I JK (3) IF(N.LE.O.OR.N.GT.NNOD) THEN WRITE (NTM. 2300) N WRITE (NOT. 2300) N ERR-. TRUE. GO TO 30 ENDIF

NN - HPNPR+N-1

C CHECK BALNK COMMON SPACE

IF (NN . GT . IDIR) THEN WRITE (NTM, 2310)

WRITE (NOT. 2310)

ERR= . TRUE

.

RETURN

ENDIF NPRT - NPRT 4- 1 3 4 IA (NN) - N GO TO 30

C REDEFINE DIRECTORY FOR ACTUAL SIZE OF NPR ARRAY

36 CALL DELETE (’NPR ’) CALL DEFINE (

’ NPR ’. MPNPR , NPRT, 1

)

RETURN

2100 FORMAT ( ’ •*** ERROR : No report control Information found.’) 2200 FORMAT

(

.’ — Solution will be written to binary file <f llename> . TMP

. In form (TIME. T (N) ,N=l,NNOO) . ’

)

2300 FORMAT ( ’ ••*• ERROR: (Generated) Node, ’15* selected for printable

. output Is out of range.’) , 2310 FORMAT (

» .’ ••** ERROR: Insufficient storage to form print control array.

.* Increase IA (MTOT) by ’,16)

END

CTDOF FUNCTION TDOF (NEQ, ID. NNOD) C-FUN : TDOF - DETERMINES IF EQUATION NUMBER NEQ IS A TEMP DOF LOGICAL* 1 TDOF INTEGER ID (NNOD) TDOF « .FALSE. DO 10 N-l.NNOO

IF ( (ID (N) . LT . 0) .AND. (ABS (ID (N) ) . EQ. NEQ) ) THEN

TDOF - .TRUE. RETURN ENDIF

10 CONTINUE

RETURN

END

C* ••••*•• *••••••• •*••»**•* ft ft**************** ft •••••••••* • •• ********

C CAL - SAP LIBRARY EXTENSIONS (;**••*•****•***••••*••••*••**»••*••***••**•••••••••*•*•»•*•*•••**••**

c FREETY SUBROUTINE FREETY C—SUB:FREETY - TYPE COMMAND LINE TO SCREEN

CHARACTER* 1 LINE COMMON /ILINE/ II COMMON /CLINE/ LINE (160) COMMON /IOLIST/ NTM, NTR, NIN . NOT. ND1 . ND2 . ND3 , ND4 SAVE /ILINE/. /CLINE/. /IOLIST/

WRITE (NTM. 2000) (LINE (I) , 1*1. II) RETURN

2000 FORMAT (1X.8QA1) END D -1 8
c FINDN SUBROUTINE FINDN (SEP, NC, KEY) CHARACTER* 1 LINE. SEP (NC) COMMON /CLINE/ LINE (160) COMMON /ILINE/ II COMMON /IOLIST/ NTM, NTR, NIN , NOT, ND1 , ND2 , ND3 , ND4 SAVE /CLINE/. /ILINE/, /IOLIST/

C FIND SEPARATOR OF NC CHARACTERS IN INPUT FILE

KEY=0 REWIND NIN

50 CONTINUE

READ (NIN, 1000, ERR=200, END=200) (LINE(I) ,I=1,NC)

II * NC CALL UPPER DO 60 N-l.NC

60 IF (SEP (N) .NE. LINE (N) ) GO TO 50 GO TO 900

200 KEY=1

900 RETURN

1000 FORMAT (80A11 END

C SAVE SUBROUTINE SAVE C—SUB: SAVE - WRITES INCORE DATA BASE TO <FILENAME> . DBS CHARACTER EXT*1. FIN*1 COMMON MTOT. NP., IA (10000) COMMON /DBSYS/ NUMA, NEXT. IDIR, IP (3) COMMON /IOLIST/ NTM. NTR. NIN . NOT. ND1 , ND2 , ND3 , ND4 COMMON /PARC/ FIN (12) . EXT (80) SAVE /DBSYS/, /IOLIST/. /PARC/

CALL FOPEN (ND1, 'DBS • .'NEW') WRITE (ND1) NUMA. NEXT, IDIR. MTOT. IP WRITE (ND1) (IA (I) , 1=1 , NEXT) . (IA ( J) . J=IDIR, MTOT) CALL FCLOSECNDl)

RETURN END

c ZEROI SUBROUTINE ZEROI (IA. NR. NC) DIMENSION IA (NR, NC) DO 10 1=1, NR

DO 10 J=1 , NC

IA (I, J) = 0 10 CONTINUE

RETURN

END

C ZEROR SUBROUTINE ZEROR (A. NR. NC) REAL* 8 A (NR,NC) DO 10 1=1, NR DO 10 J=1 , NC A (I , J) = 0.0

10 CONTINUE

RETURN

END

c SYMBC SUBROUTINE SYMBC (C . B. MAX , MB . ND . IOP

)

C—SUB: SYMBC - SYWETRIC COMPACT-FORM EQUATION SOLVER

C SOLVES SYMMETRIC BANDED EQUATIONS;

C

C [C) (X) = ( B) C

C IN-CORE WITH [C] STORED IN COMPACT FORM C (ND, MB)

.

C SOLUTION X (ND) IS WRITTEN OVER B (ND) . C

C MAX = MAXIMUM DOF FOR PARTIAL L*D*U DECOMPOSITION

C MB = (HALF) BANDWIDTH

C ND = NUMBER OF DEGREES OF FREEDOM

C IOP = 1 : COMPLETE SOLVE

C IOP * 2 : RESOLVE

C

IMPLICIT REAL* 8 (A-H.O-Z) DIMENSION C (ND, MB) , B (ND)

C-l —PERFORM THE L*D*U DECOMPOSITION OF C AND FORM Y NBM - MB

DO 300 M*2 , MAX

N - M-l

ir(C(N,l) .EQ.0.0) OO TO 300 WN - B(N) NNB - MIN ( (N+NBM) .MAX) KJ - NNB - N + 1 II - 1 DO 200 I - M, NNB

XX XX 1 IF (C (N, II) .EQ.0.0) GO TO 200 GIN - C (N, II) /C (N, 1) B (I) = B (I ) - GIN*WN

IF (IOP . EQ. 2) GO TO 200 J = 0 DO 100 K-II.KJ

J = J + 1 100 C(I,J) = C(I,J) - C (N , K) *GIN

200 CONTINUE

300 CONTINUE

C-2 --PERFORM BACKSUBSTITUTION

M » MAX

400 IF (C (M, 1) .NE.0.0) B (M) «B(M)/C(M.l)

M = M - 1 IF(M.LE.O) GO TO 600 KJ = MIN (MB. MAX-M+1) J = M DO 500 N=2 , KJ

J = J + 1 500 B(M) = B (M) - C (M, N ) *B ( J) GO TO 400

600 RETURN

END

C SYMBCX SUBROUTINE SYMBCX (C. B. MAX. MB. ND, IOP) C—SUB: SYMBCX - SY METRIC COMPLEX COMPACT-FORM EQUATION SOLVER

C SOLVES SYMMETRIC BANDED COMPLEX EQUATIONS;

C

C (C*J (X*) = (B*) cC IN-CORE WITH (C*) STORED IN COMPACT FORM C (ND. MB)

C SOLUTION X* (ND) IS WRITTEN OVER B* (ND) . C

C MAX » MAXIMUM DOF FOR PARTIAL L*D*U DECOMPOSITION

C MB - (HALF) BANDWIDTH

C ND » NUMBER OF DEGREES OF FREEDOM

C IOP - 1 : COMPLETE SOLVE

C IOP - 2 : RESOLVE

C

IMPLICIT COMPLEX (A-H.O-Z) DIMENSION C (ND. MB) , B (ND) C-l—PERFORM THE L*D*U DECOMPOSITION OF C AND FORM Y NBM = MB -

DO 300 M=2 , MAX

N = M-l IF (C (N , 1 ) . EQ . ( 0 . 0 , 0 . 0 ) ) GO TO 300

WN = B (N) NNB = MIN( (N+NBM) .MAX) KJ = NNB - N + 1 II = 1 DO 200 I = M. NNB

XI - XX + 1 IF(C(N.II) .EQ. (0.0. 0.0) ) GO TO 200 GIN = C(N.II)/C(N,1)

B (I ) = B (I ) - GIN*WN

IF (IOP . EQ. 2 ) GO TO 200 J = 0 DO 100 K=II , KJ

J « J + 1 100 C(I.J) = C(I,J) - C (N, K) *GIN

200 CONTINUE

300 CONTINUE

C-2- -PERFORM BACKSUBSTITUTION

M = MAX

400 IF(C(M. 1) .NE. (0.0,0. 0) ) B (M) =B(M)/C(M,1)

M = M - 1 IF(M.LE.O) GO TO 600 KJ = MIN (MB. MAX-M+1)

J - M DO 500 N-2.KJ

J = J + 1 500 B (M) = B (M) - C (M, N) *B ( J) GO TO 400

600 RETURN

END

CROUT SUBROUTINE CROUT (A, B. N, ND, LD. KEY) C

C GENERAL EQUATION SOLVER - M.I. HOIT 1/17/83

C SYMMETRIC AND NON- SYMMETRIC

C

C KEY - 0 TRIANGULARIZE AND SOLVE

C KEY = 1 TRIANGULARIZE ONLY

C KEY - 2 FORWARD AND BACK SUBSTITUTION ONLY

C IMPLICIT REAL* 8 (A-H. O-Z) DIMENSION A (ND, ND) , B (ND, LD) C

II - 1 IF (A(l.l) .NE.0.0) GO TO 10 PAUSE ' •*** ERROR: SUB: CROUT: Attenpt to solve singular system. GO TO 900 10 IM = 1 IF (KEY. EQ. 2) GO TO 200

IF (N . EQ . 1 ) GO TO 200 DO 140 J-2.N

JM - J - 1 A ( J, 1) - A ( J, 1 ) /A (1 , 1) “

IF ( J .EQ. 2) GO TO 130 DO 120 1-2. JM

IM - I - 1 DO 100 K=1 , IM

C FORM L (I, J) 4 U(I.J)

A(J.I) = A ( J, I) - A ( J, K) *A (K, I) 100 A(I.J) = A (I , J) - A (I , K) *A (K. J) D = A (I , I) C CHECK IF A ZERO IS ON DIAGONAL

IF(D. NE.0.0) GO TO 120 PAUSE * *•** ERROR: SUB: CROUT: Attempt to solve singular system. D -1 9
GO TO §00

120 A(J.I> - A(J,I> /D 130 DO 140 K-1,JH

140 A(J.J) - A(J.J) - A ( J, K) • A <K, J) C FORWARD AND BACK SUBSTITUTE

200 IF (KEY . EO . 1 ) GO TO 900 DO 400 L-l.LD

C FORM T (1* L) IF(N.EQ.l) GO TO 300 DO 210 1-2.

IN - I - 1 DO 210 K-l.IH

210 B(I,L> - B(I.L) - A (I , K) *B (K, L) C FORM X (X« L) 300 B(N,L) - B (N.L) /A(N.N) IFOI.EO.l) GO TO 400 IN - * - 1 DO 330 1-2.

JH « IN 1 DO 310 J « JM.N

310 B (IN, L) - B(IH.L) - A (IM. J) *B ( J. L) 320 B(IH.L) - B (IN. L) /A (IN. IN) 330 IN - IN - 1 C

400 CONTINUE

900 RETURN END
NBS-114A irev. 2»bc) U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

4. TITLE AND SUBTITLE

PUBLICATION OR 2 REPORT NO.

. Performing Organ. Report No. 3. Publication Date NISTIR 88-3868 OCTOBER 1988

A Discrete Thermal Analysis Method (DTAM) for Building Energy Simulation with

DTAM1 Users Manual

5. AUTHOR(S)

James Axley

6. PERFORMING ORGANIZATION (If joint or other than NBS. see instructions) 7. Contract/Grant No. national bureau of standards

DEPARTMENT OF COMMERCE 8. Type of Report & Period Covered

WASHINGTON, D.C. 20234

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

National Institute of Standards and Technology

Gaithersburg, MD 20899

U.S. Department of Energy

Washington, D.C.

10 . supplementary notes

|

I Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here) This document includes a report that describes the theoretical basis of the program

DTAM 1 and a users manual for the program. DTAM1 is a general purpose building energy

simulation based upon- discrete analysis techniques, including, but not limited to, the

Finite Element Method, used in orther fields of physical simulation. It is the product

of a first phase of development of Discrete Thermal Element Analysis Techniques for

Building Energy Simulation that are expected to provide a means to unify existing

building energy simulation theory. DTAM1 provides a library of discrete thermal

elements, that may be assembled to model thermal systems idealized to have constant

material and heat transfer properties (i.e., linear idealizations), including:

- ID two-node thermal resistance elements

- single-node lumped capacitance elements

- two-node fluid flow loop element

- ID two-to-four node isoparametric conduction Finite Elements

- 2D four-node isoparametric conduction Finite Elements (planar and axisymmetric)

Equations defining a variable node mean radiant temperature element are also presented

in the report. Steady state and transient analysis capabilities are included.

Temperature, heat flow rate, and convective boundary conditions may be modeled and

s-stem temperature variables may be constrained to be equal so that mixed assemblages

of ID and 2D elements may be employed.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons

)

building energy simulation, building dynamics, computer simulation techniques, discrete

analysis techniques, discrete thermal elements, dynamic simulation, finite element

analysis, DTAM1

.

13. AVAILABILITY 14. NO. OF

PRINTED PAGES

fin Unlimited

I 1 For Official Distribution. Do Not Release to NTIS 79

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

• 20402. 15. Price

fxl Order From National T echnical Information Service (NTIS), Springfield, VA. 2216 1 $13.95

U S COMM- D C 6 04 3»P 8

0
.

'

.