NATL INST. OF STAND & TECH R.I.C.

A 111 DM 7131TE

NfST

PUBLICATIONS

\ NIST SPECIAL PUBLICATION 400-96

U.S. DEPARTMENT OF COMMERCE/Technology Administration

National Institute of Standards and Technology

QC

100

U57

JO. 400-96

995

Semiconductor Measurement Technology:

HOTPAC: Programs for Thermal Analysis

Including Version 3.0 of the TXYZ Program,

TXYZ 30, and the Thermal MultiLayer

Program, TML
Jhe National Institute of Standards and Technology was established in 1988 by Congress to "assist industry

in the development of technology . . . needed to improve product quality, to modernize manufacturing processes,

to ensure product reliability . . . and to facilitate rapid commercialization ... of products based on new scientific discoveries." NIST, originally founded as the National Bureau of Standards in 1901, works to strengthen U.S. industry's

competitiveness; advance science and engineering; and improve public health, safety, and the environment. One of the agency's basic functions is to develop, maintain, and retain custody of the national standards of measurement, and

provide the means and methods for comparing standards used in science, engineering, manufacturing, commerce,

industry, and education with the standards adopted or recognized by the Federal Government.

As an agency of the U.S. Commerce Department's Technology Administration, NIST conducts basic and applied

research in the physical sciences and engineering, and develops measurement techniques, test methods, standards, and

related services. The Institute does generic and precompetitive work on new and advanced technologies. NIST's

research facilities are located at Gaithersburg, MD 20899, and at Boulder, CO 80303. Major technical operating units and their principal activities are listed below. For more information contact the Public Inquiries Desk, 301-975-3058.

Office of the Director

• Advanced Technology Program

• Quality Programs

• International and Academic Affairs

Technology Services

• Manufacturing Extension Partnership

• Standards Services

• Technology Commercialization

• Measurement Services

• Technology Evaluation and Assessment

• Information Services

Materials Science and Engineering

Laboratory

• Intelligent Processing of Materials

• Ceramics

• Materials Reliability'

• Polymers

• Metallurgy

• Reactor Radiation

Chemical Science and Technology

Laboratory

• Biotechnology

• Chemical Kinetics and Thermodynamics

• Analytical Chemical Research

• Process Measurements

• Surface and Microanalysis Science

• Thermophysics2

Physics Laboratory

• Electron and Optical Physics

• Atomic Physics

• Molecular Physics

• Radiometric Physics

• Quantum Metrology

• Ionizing Radiation

• Time and Frequency'

• Quantum Physics'

Manufacturing Engineering Laboratory

• Precision Engineering

• Automated Production Technology

• Intelligent Systems

• Manufacturing Systems Integration

• Fabrication Technology

Electronics and Electrical Engineering

Laboratory

• Microelectronics

• Law Enforcement Standards

• Electricity

• Semiconductor Electronics

• Electromagnetic Fields'

• Electromagnetic Technology'

• Optoelectronics 1

Building and Fire Research Laboratory

• Structures

• Building Materials

• Building Environment

• Fire Safety

• Fire Science

Computer Systems Laboratory

• Office of Enterprise Integration

• Information Systems Engineering

• Systems and Software Technology

• Computer Security

• Systems and Network Architecture

• Advanced Systems

Computing and Applied Mathematics

Laboratory

• Applied and Computational Mathematics2

• Statistical Engineering2

• Scientific Computing Environments2

• Computer Services

• Computer Systems and Communications2

• Information Systems

'At Boulder, CO 80303.

2 Some elements at Boulder. CO 80303.
Semiconductor Measurement Technology:

HOTPAC: Programs for Thermal Analysis

Including Version 3.0 of the TXYZ Program,

TXYZ 30, and the Thermal MultiLayer

Program, TML

John Albers

Semiconductor Electronics Division

Electronics and Electrical Engineering Laboratory

National Institute of Standards and Technology

Gaithersburg, MD 20899-0001

August 1995

U.S. DEPARTMENT OF COMMERCE, Ronald H. Brown, Secretary

TECHNOLOGY ADMINISTRATION, Mary L. Good, Under Secretary for Technology

NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, Arati Prabhakar, Director
National Institute of Standards and Technology Special Publication 400-96

Natl. Inst. Stand. Technol. Spec. Publ. 400-96, 92 pages (Aug. 1995)

CODEN: NSPUE2

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1995

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, DC 20402-9325
This Special Publication is dedicated to Frank F. Oettinger

on the occasion of his retirement from Federal service.

Frank was instrumental in beginning the

NBS/NIST effort on the Kokkas model and the TXYZ codes.

He has also been supportive of the body of work presented here as well as

its integration into other Semiconductor Electronics Division activities.

iii

Semiconductor Measurement Technology:

HOTPAC: Programs for Thermal Analysis Including

Version 3.0 of the TXYZ Program, TXYZ30,

and the Thermal MultiLayer Program, TML

Table of Contents

Page

Abstract 1

Part 1. General Elements of Thermal Model 3

1.1 - Introduction 3

1.2 - Steady-State Heat Flow: Single Rectangular Layer 3

1.3 - Steady-State Heat Flow: Multilayer Rectangular Structure 8

1.4 - Kokkas Three-Layer Model 10

1.5 - Kokkas Three-Layer Model: Form of the Function U(n,m) 16

1.6 - Kokkas Three-Layer Model: Behavior of Fourier Coefficients

for Small Values of the Argument 17

1.7 - Kokkas Three-Layer Model: Behavior of Fourier Coefficients

for Large Values of the Argument 19

1.8 - Special Case of Power Source Covering Top Surface 23

1.9 - Effect of Upper Summation Limits on Temperature 24

1.10 - The TXYZ and TXYZ20 Codes 24

Part 2. Analytic Evaluation of Line and Area Averages 25

2.1 - Introduction 25

2.2 - Calculation of Average Temperatures 26

2.3 - Numerical Implementation and TXYZ30 30

2.4 - Total Area Averages and the Thermal Resistance 31

2.5 - Applications of Area Averages 31

Part 3. Recursion Relation Solution of Multilayer Model 31

3.1 - Introduction 31

3.2 - The N-Layer Electrical Problem 33

3.3 - The N-Layer Thermal Problem 38

3.4 - The Thermal Recursion Relation 38

3.5 - Relation to Previous Solutions 39

3.6 - Numerical Implementation of Thermal MultiLayer Code, TML 41

3.7 - Generalized One-Dimensional Thermal Resistance 42

3.8 - Application to Buried Oxide Structure 42

3.9 - Summary and Conclusions 44

Availability of HOTPAC Software Package 45

Acknowledgments 46

References 47

v
Page

TXYZ30 I/O File Listing - txyz30io.l 50

TXYZ30 I/O File Listing - txyz30io.2 51

TXYZ30 I/O File Listing - txyz30io.3 52

TML I/O File Listing - tmlio.l 53

TML I/O File Listing - tmlio.2 54

TML I/O File Listing - tmlio.3 55

TML I/O File Listing - tmlio.120 56

TML I/O File Listing - tmlio.220 58

TML I/O File Listing - tmlio.320 60

Appendix A - TXYZ30 Listing 62

Appendix B - TML Listing 77

List of Figures

1. Geometry of the multilayer structure in which the steady-state

temperature is calculated 7

2. Geometry of the structure in which the line average or area average

temperature is calculated 27

3. Top view of the surface of the structure shown in figure 2 28

4. Schematic representation of the geometry used in the Schumann and Gardner

multilayer Laplace equation analysis 34

5. Schematic representation of the geometry used in multilayer analysis

of the steady-state heat flow problem 37

6. Example of the small increase in computation time associated with

increasing the number of layers in the thermal structure 40

7. Example of the application of the thermal recursion relation to the calculation

of the surface temperature of an SOI (Silicon-On-Insulator) structure 43

VI
Semiconductor Measurement Technology:

HOTPAC: Programs for Thermal Analysis Including

Version 3.0 of the TXYZ Program, TXYZ30,

and the Thermal MultiLayer Program, TML

John Albers

Semiconductor Electronics Division

National Institute of Standards and Technology

Gaithersburg, MD 20899

Abstract

This report presents and discusses a number of recent developments in the steady-state

thermal analysis of multiple-layer structures. These include: 1) analytical evaluation of

line and area average temperatures and 2) a recursion relation technique for calculating

the steady-state surface temperature of a multilayer structure with an arbitrary number

of layers. The application of the analytic averaging to the TXYZ code is incorporated in

the updated code, TXYZ30, while the multilayer recursion relation solution along with

analytic averaging are included in the Thermal MultiLayer code, TML. Both of these are contained in the HOTPAC software package.

The first part of this report contains a discussion of the general elements of the multiplelayer thermal model. This is presented in some detail for the case of the Kokkas three-layer

problem and the associated TXYZ code. The previous update of the code, TXYZ20, is also

discussed. This incorporates more flexible handling of input data, assignment of positive or negative noninteger weights to the various heat sources or heat sinks, and improved

evaluation of limiting forms in the code.

The second part of the report presents the analytical evaluation of line and area averages

which can be calculated directly. The analytical calculation of the area average temperature should provide for a more direct and convenient connection to the experimentally

measured temperature which tends to average over the measurement area. The line and

area averaging are incorporated into the TXYZ code to produce the TXYZ30 update.

The third portion contains a detailed discussion of the calculation of the surface temperature of a multilayer structure with an arbitrary number of layers. This is based on a

recursion relation technique previously employed in electrical spreading resistance analysis

to determine the surface potential from the multilayer electrical Laplace equation. The

line and area averaging techniques can be applied directly to the recursion relation solution

and are included in the Thermal MultiLayer, TML, code.

The appendices contain the listing of the annotated, internally documented FORTRAN

source codes for TXYZ30 and TML. These codes as well as several sample input and output

data files are available in ASCII format on DOS-formatted floppy disks. The sample input

and output data files are included so that the user can check the programs for proper

operation as well as become familiar with the setup and use of the codes. Users of the

1
previous versions of the TXYZ code should find the TXYZ30 and TML codes easy to use and should benefit from the wider range of problems which they can be used to address.

Key words: average temperature; FORTRAN programs; Fourier analysis; Laplace equation; multilayer model; semiconductor devices; semiconductor materials; steady-state heat

flow; thermal conductivity.

2
PART 1. GENERAL ELEMENTS OF THERMAL MODEL

1.1 - INTRODUCTION

The operation of a semiconductor device relies on the passage of current through portions

of the structure. This is accompanied by power dissipation and heating which gives rise

to a temperature distribution. In time, thermal stresses may arise and lead to device

degradation and possible failure. In addition to the effects of thermal stresses on device

reliability, the modeling of the steady-state thermal response of a system is also useful

since the temperature distribution may have an effect upon the local mechanical, thermal,

electrical, or optical properties of the system. Consequently, an accurate physical model

of the temperature distribution under the power condition of actual operation is of great

importance. It is also possible to ascertain the relative effects of composition (thermal

conductivity) and geometry (layer thickness). By carefully optimizing composition and

geometry, it may be possible to minimize the thermal stress and hence ensure optimal

device lifetime.

The physical and mathematical model used here is taken from the work of Kokkas which

has been used previously to construct the TXYZ programs. The TXYZ computer program

has been used for a number of years for the thermal analysis of semiconductor devices and

packages. This program makes use of the closed form, Fourier series solution of the steadystate heat flow equation for the general case of a rectangular three-layer structure with

multiple heat sources on the top surface. TXYZ provides for the calculation of the temperature at any set of points in this structure and has proven useful for the determination

of the steady-state temperature distribution of semiconductor chips and packages.

The purpose of this work is to report on several recent advances which have been made

in the thermal analysis. These include the analytical evaluation of line and area averages

as well as the recursion relation solution of the steady-state surface temperature of a

mulitlayer structure with an arbitrary number of layers.

1.2 - STEADY-STATE HEAT FLOW: SINGLE RECTANGULAR LAYER

Consider a material of uniform thermal conductivity, K\ , in the form of a rectangular box of

lateral dimensions Lx ,L y

, and thickness L\ . The problem is to determine the temperature,

T(z, y, z), inside the material. The temperature is assumed to satisfy the steady-state heat

flow equation [1]

V

2 T(z, *,,.,) =0. (1)

As this equation is second order in the three coordinates, there are six boundary conditions.

Four of these will be provided by the lateral boundary conditions. In the present problem,

all four of the lateral boundary conditions are provided by the assumption that the lateral

surfaces are adiabatic; i.e., there is no heat flow out of the lateral boundaries of the material;

3
dT(x,y,z)

dx x=0,L,

dT(x,y,z)

dy = 0. (2)

y=0,L l

The remaining two boundary conditions will be provided by the vertical boundary conditions (in z). These vertical boundary conditions will not be specified at the present time

as the intent of the present section is to obtain a general solution of the one-layer problem

where only the lateral boundary conditions are specified. As the above equation is formulated in Cartesian coordinates, it is convenient to use Fourier analysis techniques to solve

the x and y portion of the equation. The Fourier transform with respect to the variables x and y is used, remembering that the geometry is constrained to Q,LX , and 0,Ly

. This

is defined as [2]

UxJy ,z)= l I T(x,y,z)exp(-2iri(xfx + yfy ))dxdy,

Jo Jo

(3)

where fx ,fy are the Fourier transform variables which are conjugate to the variables x,y.

The inverse Fourier transform is defined as

/• +oo />+oo

/ T{fx ,fy ,z)exp(27ri(xfx + yfy ))dfx dfy

. -oo J —oo

(4)

The requirement that there is no heat flow out of the sides of the structure, i.e., dT(x, y, z)/dx

and dT(x,y, z)/dy are zero when x and y are either equal to zero or to Lx or Ly , respectively, leads to a consideration of the expression

dT(x,y,z)

dx

/ / r(fx ,fy

, z)exp(2Kiyfy )2Trfx l -sin(27rfx x) + icos(27rfx x)>dfx dfy , (5)

J —oo J —oo

and a similar expression for dT(x,y, z)/dy. If this is to be zero at the origin, this would

require that the cosine term be removed. Next, consider the resulting expression at the

other lateral boundary, i.e., at x = Lx where it is supposed to be zero. The only way that

this could be the case is if the argument of the sine function is an integer times 7r, or that

the Fourier transform variable is of the form

4
The same argument applies to the y-dependent portion. The Fourier representation of the

temperature with the above lateral boundary conditions may then be written as

/+OO /'+00

/ T(n,m,z)cos(n7cx/Lx )cos(mTry/Ly)dfx dfy

. (7)

-oo J —oo

The Fourier transform equation is now written as

r(n,m,2)= / / T(x, y, z) cos(mrx/Lx ) cos(rmry/ Ly )dxdy. (8)

Jo Jo

As the system is of finite size, it is convenient to write the integral in eq (7) as a sum over

the Fourier cosine terms which fit into the rectangular geometry. Further, as the cosine

function is symmetric around the origin, the sums may be written over only the positive

values of the indices. It is important to keep in mind that the terms corresponding to

m,n=0 do not have the factor of 2 coming from the symmetry of the cos function. In

addition, the differentials may be written as

» n n + 1 n 1 dfx = A— = — — = —. (9)

±j x Li x JU X JU X

The Fourier representation of the temperature may be written as

7Y x _ V" V" 4T(n,m,z)cos(mrx/Lx )cos(miTy/Ly)

hho (*»o + l)(«mo + l)Lx Ly

' (1U)

where 8 nn> is the Kronecker delta and is equal to unity if n — n' and zero otherwise. By

substituting eq (10) into eq (1) and using

„2

( d

2 d2 d2 \

and

d

2 ——- cos(nirx /

L

x ) = —(nw/Lx ) 2 cos(nirx/

L

x ), (12)

ox1

and the same relation for the y-dependence, it is straightforward to show that

5
4 cos( mrx / Lx ) cosimiry/

v ) { . ._ , 2

. ,_ d

2

1 . .

n=0 m=0 v / v / » v (13)

As the sum is zero for arbitrary values of the variables x and y (and the cosine terms are, in general, nonzero), then a necessary and sufficient condition that eq (13) is satisfied is that

{-(nn/Lx ) 2 - (rmr/Ly ) 2 + J^}r(n,m,z) = 0. (14)

This differential equation may be solved analytically using elementary methods. If the

variable, 7, is defined as H(D2+

(^)T- (i5)

eq (14) may be rewritten as

d

2 r(n, m, z) — 7

2 t(ti, m, z) = 0. (16)

<9z 2

The solution of this equation is

r(n, m, z) = a cosh(7^) + /? sinh(72), (17)

where the coefficients a and /?, which may be functions of 7, are determined from the two

2-dependent boundary conditions.

The above equation is the general solution for the z-dependent Fourier expansion coefficients for a single rectangular layer. It provides a convenient basis set for the discussion

of the problem of a rectangular multilayer structure where all of the layers have the same

lateral dimensions. For the multilayer case, the solution in each of the layers can be ex- pressed in the form of the above equation where the coefficients are to be determined from

the two z-dependent boundary conditions appropriate to each of the layers. This is used

in the next section where the multilayer problem is discussed and then specialized to the

three-layer case.

C
Figure 1. Geometry of the multilayer structure in which the steady-state temperature is calculated.

7
1.3 - STEADY-STATE HEAT FLOW: MULTILAYER RECTANGULAR STRUCTURE

The basic problem considered in this section is the general approach to the calculation of

the temperature in a multilayer structure depicted in figure 1 . The layers are characterized

in terms of the thermal conductivities and thicknesses Kl: L{(i = 1,2, . . .,iV), where the

numbering begins at the bottom layer and proceeds up to the top or N-th layer. Then,

the bottom layer will be layer number 1 with temperature T\(x,y, z). The next layer up

will be layer number 2 with temperature T2(x, y, z). This continues to the top layer which

will be layer N with temperature Tw(x,y, z).

The equations are discussed in general for the N-layer problem and then specialized to the

three-layer problem as originally described by Kokkas [3,4].

This will allow for the discussion of the Kokkas three-layer solution and the TXYZ code

while providing the necessary development for the discussion of the recursion relation in

the third part of the report.

The thicknesses and thermal conductivities of these layers are of considerable importance

in the dissipation of heat generated by the power sources on the surface of the top layer.

These power sources are typically the regions at or near the surface of the device where

currents are passed into the device during normal operation. Consequently, the generation

of heat in the device is one of the unavoidable side effects of device operation.

The mathematical formulation of this problem is based upon the following set of assumptions:

1) the lateral dimensions of all the layers in the structure are equal while the thicknesses

may be different;

2) each layer is of uniform, isotropic, temperature-independent thermal conductivity;

3) there is no heat loss from the lateral surfaces due to either radiation or convection

(adiabatic surfaces) - heat flow in the structure takes place by conduction;

4) there is no input power density inside the structure - heat is generated only on the

top surface;

5) there are no heat losses due to interconnections to the top layer; and

6) the heat sink, which is in contact with the bottom layer, is ideal and has a temperature equal to ambient.

The temperature in each of the layers is assumed to satisfy the steady-state heat flow

equation,

v

2 r(^,t/,5) = o. (is)

8
The general one-layer solution provides a convenient and useful basis for the multilayer

solution. Then, the temperature in each layer may be written as

Ti(x,y,z)

The Fourier coefficients for the solution of the steady-state heat flow equation in the i-th

layer of an iV-layer structure may be written as

r,(n, ra, z) = T{(jz) — oti cosh(7z) + sinh(72), (i = 1, iV), (20)

where the expansion coefficients, a; and are determined from the z-dependent boundary

conditions.

Notice that the assumption that there is no heat flow out of the lateral boundaries,

-0, (i = l,2, . . .,i\0, (21)

is already contained in eq (19) as discussed in the previous single-layer analysis.

The 2N ^-dependent boundary conditions used to solve the set of equations in eqs (19)

and (20) may be expressed as follows.

First, heat enters (or leaves) the structure through the portions of the top layer where

power is applied. This is expressed as

= P(*,y), (22)

2 =0

where kn is the thermal conductivity of the top layer, and where P(x,y) is the power

function. This is expressed as P(x,t/) = PoU(x,y), where Pq is the power density (W/cm2

), and U(x,y) describes the surface geometrical distribution of the power sources. Next, the

bottom-layer boundary condition is provided by the requirement that the temperature is continuous across the interface between the bottom layer and the heat sink and equal to

the ambient which is taken as zero,

oo

. 4Tj(n, m, z) cos(mrx /

L

x ) cos(rmry/

L

y )

07^0 ' («nO+l)(«mfl+lM

'

dTt (x,y,z)

dx

dTi{x,y,z)

x=0.L, dy

dTN (x,y,z)

Ti(x,y,z) = Ta = 0,

9

(23)
where the function is evaluated at the interface. Finally, the remaining 2(7V — 1) conditions

are provided by the requirement that the temperature and the normal component of the

heat flow are continuous across the internal interfaces. These are expressed as

Ti(x, y, z) = T? _j (x, y, z), (24)

oz oz

where the functions and their derivatives are to be evaluated at the interfacial boundaries.

For the case of an TV-layer structure, the substitution of eqs (19) and (20) into the boundary

conditions given by eqs (22 - 25) gives rise to a set of 2N equations in 2iV unknowns

({a;(7)}, {fliij)}, i — l,...,iV). The analytic solution of this system of equations requires

the use of matrix algebra involving Cramer's rule [5] and the Laplace method [5] for the

evaluation of the resulting determinants. Clearly, this also can become rather tedious,

especially since the expansion coefficients are functions of the continuous variable, 7.

For cases up to i = 3, Kokkas [3,4] was able to work out the system of equations for not

only the surface temperature but also the temperature inside the three-layer structure.

The Kokkas three-layer solution is reviewed in the next section and the N-layer problem

is returned to in Part 3.

1.4 - KOKKAS THREE-LAYER MODEL

(25)

The Kokkas three-layer model [3,4] involves the solution of the multilayer equations, eqs

(19, 20, 22 - 25), for the case where N=3. This can be carried out analytically as the

systems of equations are still tractable. Specializing the multilayer equations to the case

of three layers yields the Fourier coefficients

Tj(n, m, z) = Ti(jz) = at cosh(7.?) + fa sinh(72), (i = 1,2, 3). (26)

The corresponding boundary conditions are

*3

dTN (x,y,z)

z=0

T3 (x,y,z) =T2 (x,y,z)\

T2 (x,y,z)

z=-L

z=-(L3 +L2 ) = Ti(x,y,z)

z= — L3

z=-(L3 + L2 )

(27)

(28)

(29)

10
«3

dT3 (x,y,z)

dz z= —Lz «2

dT2 (x,y,z)

8z

(30)

Z= — L;

*2

dT2 (x,y,z)

dz = «1

z=-(L3 + L 2 )

dTx {x,y,z)

dz

(31)

z=-(L3 + L 2 )

Ti(x,y,z)

z— — L, = Ta = 0, (32)

where all temperatures are measured relative to the ambient heat sink temperature and

Lz — L$ + L2 + L\.

It is important to note that the origin of the depth scale is at the surface of the top layer

and that all vertical distances are negative.

The solution of the system of equations is by standard matrix algebra techniques. However,

instead of using this method explicitly here, the Fourier coefficients for each of the layers are presented and are shown to satisfy the heat flow equation and the appropriate boundary

conditions. In particular, the Fourier coefficients in the three layers are presented in

reference [4]. Specializing these to the steady-state situation, they are

r3 (n,m,2) = Aj£cosh(7(L3 + z)) + Csinh(7(L3 + z))j, (33)

where

r2 (n,m,z) = A^D cosh(7(L3 + L2 + z)) + E sinh(7(L3 + L2 + z)) }, (34)

Ti(n,m, z) = Asinh(7(L3 + L2 + L\ + z)), (35)

A =

U(n, m)Po

k37 ( B sinh( 7L3 ) + C cosh(7L3 ) }

(36)

B = Dcosh(7L2 ) + £sinh(7L2 ), (37)

(38)

11
D = sinh(7l/i), (39)

E = — cosh(7jLi), (40)

«2

and

'L>x r^y

U(n,m) = / / U(x,y)cos(nirx/Lx )cos(miry/

y )dxdy, (41)

Jo Jo

is the double Fourier cosine transform of the power density uniformity function.

Now, it is shown that the above are the solutions of the z-dependent part of the steady-state

heat flow equation (see eq (16))

d

2 Ti(n,m,z) - j

2 Ti(n,m,z) = 0, (42) d

where the subscript i takes on the values of 3,2,1. This may be easily shown to be the case

as

d

2 cosh(7(L + z)) = 7

2 cosh(7(£ + z)) , (43) 2 dz

and

d

2

dz<

sinh(7(i: + z)) = 7

2 sinh(7(L + *)), (44)

where L is a constant and is equal to L3, L3 + X2, or £3 + £2 + L\ in eqs (33 - 35). Hence,

eqs (33 - 35) satisfy the z-dependent differential equation. Next, it is shown that these

Fourier coefficients satisfy the appropriate boundary conditions. The first of these is that

drgfu, m, z)

«3-

dz = P(n, m) = U(n, m)P0 . (45)

2 =0

Using the Fourier coefficient given by eq (33) for the top layer, this may be evaluated as

12
«3

dr3 (n,m,z)

2 =0

ac 3 7^4< Z?sinh(7.L3 ) + C cosh(7.L3 )

«37

U{n, m)P0 K37 [ 5sinh(7L3 ) + Ccosh(7L3 ) #sinh(7L3 ) + C cosh(7i:3 ) )

f7(n, m)Po (46)

Hence, the top-layer boundary condition is satisfied by the T3 (n,m,z). Next, consider the

bottom-layer boundary condition, i.e.,

Ti(n,m,z)

z=-(L3 + L 2 +Ll ) (47)

Making use of eq (35), this may be readily evaluated as

ri(n,m,z)

z=-(L3 + L 2 +L 1 ) = A sinh(7(L3 + L2 + £1 - L3 - L2 - £1)) = 0. (48)

Hence, the last boundary condition is satisfied. The final set of boundary conditions to be

verified are the ones which pertain to the interface boundary conditions. These are

T3 (n,ra,2) =r2 (n,m,2)

z=-L z= — L3 (49)

T2 (n,m,z)

z=-(L3 +L2 ) Ti(n,m,z)

z=-(L3 + L7 ) (50)

«3

dT3 (n,m,z)

dz z= —L3

K2

dT2 (n,m, z)

z— — L3 (51)

«2

dT2 (n,m,z) = «i

z=-(L3 +L2 )

dr\{n, m, z)

z=-(L3 + L2 )

(52)

In order to verify these equations, it is simplest to calculate the right- and left-hand sides of

the equations and then compare them directly. The left-hand side of the first temperature

continuity equation may be evaluated as

13
T3 (n,m,z)

z= —L3 A^Bcosh(j(L3 - L3 )) + Csinh(7(L3 - Ls))}

= AB — a|d cosh(7L2 ) + E sinh(7i:2

) }

(53)

The right-hand side of the equation may be evaluated as

r2 (n,m,z)

2= —

s = A| JDcosh(7(L3 + L2 -L3 )) + Esinh(7(Z-3 + L2 -£3))}

= A{D cosh(7£2 ) + E sinh(7i:2

) } • (54)

Hence, the first of the temperature continuity equations satisfies the boundary condition.

Next, consider the second temperature continuity equation,

T2 (n,m,z)

z=-(L3 +L2 ) = TX (n,m,z)z=-(L3 + L2 ) (55)

The left-hand side of the equation may be evaluated as

r2 (n,m,z) = A\ D cosh(j{L3 +

2 -3 -2 )) + E sinh(y(L3 +2 -

3 -2 ))\ z—-(L3 + L2 ) { >

= AD = Asinh(7Li). (56)

The right-hand side is

ri(n, m, z)

2=-(L3 +L2 ) Asmh(j(L3 + L2 + - L3 - L2 )) = Asinh^Li). (57)

Hence, the second temperature continuity equation is satisifed. The next equation to be

verified is the first heat-flow continuity equation; i.e.,

«3

<9r3 (n,m,z)

3;

«2

dT2 (n,m,z)

2= — L? (58) z= —L3

Evaluating the left-hand side leads to

14
«3

5r3 (n,m, z)

dz z=-L: = «37^{^smh(7(JL3 -L3 )) +Ccosh(7(L3 -L3 ))}

= KzjAC = K2jA^Dsmh(-)L2 ) + £cosh(7£2 )}. (59)

The right-hand side may be evaluated as

«2

dr2 (n,m, 2:) 5; = K2 7^{^sinh(7L2 ) + £cosh(7L2 )} (60)

The final equation to be evaluated is that for the heat flow continuity between the second

and third layers; i.e.,

«2

dT2 (n,m,z)

dz

«1

dri(n, m, 2)

2=-(L3 +L2 ) z=-(L3 + L 2 )

The left-hand side of this equation may be evaluated as

«2-

<9r2 (rc,m, z)

dz 2 =_(L3 +L2 )

K 2 7A|Dsinh(7(i:3 +L2 - L3 - L2 )) + E cosh(7(£3 + L2 - L3 -£2 )) j = ac 2 7A£i = K17A cosh(7Z/i ).

The right-hand side is

(61)

(62)

dri(n, m, z)

dz = /c 1 7Acosh(7(L3+L2 +Li-L3 -JL2 ) = KijAcosh(jLi). (63)

2 =_(L3 +L2 )

Now that it has been shown that the Fourier coefficients satisfy the steady-state heat flow

problem and the appropriate boundary conditions, there are several points to be considered

before getting into the body of the program. These include: (1) evaluation of the function

J7(n,m) for a uniform power source of given size and (2) simplification of the Fourier

coefficients for .subsequent numerical analysis. The latter point is necessary as the limits

of small 7 and large 7 may give rise to overflow or underflow problems when the program

is constructed.

15
1.5 - KOKKAS THREE-LAYER MODEL: FORM OF THE FUNCTION U(n,m)

The first thing to be considered is the form of the function U(n,m) for an arbitrary number

of heat sources. This function is defined as

l-L X pLy

U(n,m) = I I U(x, y) cos(mtx/Lx ) cos(mny/ Ly )dxdy. (64)

Jo Jo

The analysis can most easily be accomplished in terms of a single uniform heat source. The

case of an arbitrary number of heat sources follows by summing the results of each heat

source. Further, if any of the heat sources are nonuniform, their effects can be constructed

by suitably overlapping a number of uniform heat sources. In the coordinate system being

used, consider a single heat source denoted by the index i with a corner at the location

(xi,yi) and lengths along the x- and y-directions given by (lx l Jyi). Over the area of the

heat source, U(x,y) is assumed to be uniform and equal to unity. Away from the area

of the heat source, U(x,y) is assumed to be equal to zero. Then, U(x,y) may be viewed

as being a unit step function over the surface of the power source. Consequently, the

contribution from this single heat source may be written as

I I U(x, y) cos(nirx/

x ) cos(miry/Ly )dxdy =

Jo Jo

rxi + lxi ryi +lyi

I I cos(nirx /

x ) cos(mny/

y )dxdy. (65)

The integrals can be simply evaluated to give the result that

TT Lx Ly f . /mr(xi + lx{)\ . /mrxi\}

Ui(n, m) = - - y - < sin —^— - sin —— > v

' ; (n7r)(m7r) \ V Lx ) V Lx )

J J. (mir(yi + lyi)\ . /rmryl \\ X

l

Sm

l Ly )- Sm

{-Ly-)j-

^

Similar expressions may be written for each of the heat sources and then summed to give

the cumulative heat source effect.

Before turning to the small 7 and large 7 behavior of the Fourier coefficients, it is important

to consider the behavior of the function U(n,m) for either n = 0, m = 0, or both. This

is important in the numerical implementation of the solutions as the program will have to

calculate the double Fourier cosine transform over the range of n,m required by the sum

in eq (19). Once the value of n or m is zero, there will be problems with most machines as

far as evaluating the seeming divergence. This can be circumvented by investigating the

16
behavior of the function for n = 0. This may be readily carried out by first noting that the

function ?7(n, m) (for the particular form of U(x, y)) is a product of two terms. This may

be simply written as U(n, m) = U(n)U(m). Then, consider the n (from the x integration)

contribution to the function which is given by

U[n) = _|sm ^ _ }_ sm (__jj. (67)

There is an apparent divergence or infinity if n is simply set equal to zero. This is the way

in which a computer would look at the expression. However, this infinity is not real as can be seen by using the expansion of the sine function for small values of the argument. In

particular,

x

3 sm(x) = x - — -\ . (68)

Making use of this expression for the sine function, it is straightforward to show that

v tt< \ v L* f

• fnir(xi +lxi)\ fnnxt \\ hm U ( n ) — urn —< sin — sin > =

n _0 n—0 717T [ V Lx / V Lx J

J hstl

lx—) ~ Kir))

= lXi

- (69)

The same conclusion holds for the y-dependent portion, i.e., U(m). This must be specially

coded to bypass any overflow problem. The specific coding of the heat source may be

found in the program listing in the function UZERO.

In general, the function U(n,m) is oscillatory and does not approach zero sufficiently fast

for large values of n or m. In particular, eq (66) shows that Ul {n,m) —> 0 like 1/nm as m,n —> oo. In addition, the cosine terms in eq (19), i.e., cos(nnx/

L

x ) cos(m7ry/L 2/ ), do

not approach a definite limit for large values of the argument. Consequently, the product

U(n, m) cos(mvx /

L

x ) cos(rmry/

L

y

) tends to zero slowly.

1.6 - KOKKAS THREE-LAYER MODEL: BEHAVIOR OF FOURIER COEFFICIENTS

FOR SMALL VALUES OF THE ARGUMENT

As has been seen in the treatment of the function U(n, m), care must be taken for the case where both n and m are equal to zero (or 7 = 0). The same considerations must be carried

out for the Fourier coefficients in the three layers. In the solutions in the three layers, the

summation indices n,m appear. The summation over these variables is of influence in

the variable 7 according to eq (15). Also, the Fourier coefficients contain the hyperbolic

functions which depend upon 7. For large values of 7, the sinh and cosh functions grow

17
exponentially. This can present special numerical problems when the summation variables

approach the upper limits which may be required for the case of very small heat sources. Hence, special care must be taken to study the behavior of the Fourier coefficients for small

7 and large 7 to remove any potential numerical overflow problems. Once this is properly

taken care of, the Fourier coefficients and the solutions will be numerically well behaved.

First, consider the small 7 behavior of the Fourier coefficients. This is done by considering

the small 7 behavior of the eqs (33 - 40) and the small argument behavior of the hyperbolic

functions. In the following discussion as well as the discussion of the large 7 behavior, the

term U(n, ra)Po

/

k3 will be removed for convenience. This term will henceforth be included

explicitly in the sum in eq (10) as the Fourier coefficients for all three layers contain this as a common factor through A (see eqs (33 - 36)). Then, for small 7,

E 11

«2

(70)

D « jLi , (71)

C^M+ -L (72)

K3 I K 2 I

B^ 1 L1 + —7L2 , (73)

K 2

and, remembering that the factor U(n,m)Po / k$ has been included explicitely in eq (19),

A«i-3

-. (74)

7 «1

Making use of these expressions and the small argument behavior of the hyperbolic functions, it is straightforward to investigate the small 7 behavior of the solutions. In particular,

. . 1 K$ [ _ K\ K2 f 9 , _ K\ 1 , _ ,1

T3 (n,m,z) w <^ 7L1 + —7Z2 + —WL2 LX + — WL3 + z) } j K-i { k2 k3

I k2 ) J w —i Li + —L2 + —(£3 + z)

«1 [ «2 «3

Then,

18
^3 ^3

lim r3 (n, m, z) = (L3 + z) H 1-2 H £i- (75) 7^0 «2 /t x

Next,

r2 (n,m,z) w - —I 7^ + —7(^3 + L2 + 2)},

7 Ki I «2 J

Then,

]im T2 (n,m,z) = —L1 + —(Ls + L2 + z). (76) 7^0 Kl K2

And finally,

1 ac 3 Ti(rc,m,z) « <[ 7(^3 + L2 + L2 + z)

7 AC!

Then,

^3

limri(n,m,2) = —(L3 + L2 + Li + 2). (77) 7^0 K\

These special forms of the Fourier coefficients (in the limit as 7 —> 0) are necessary in the

code to bypass overflow problems for small values of the argument.

1.7 - KOKKAS THREE-LAYER MODEL: BEHAVIOR OF FOURIER COEFFICIENTS

FOR LARGE VALUES OF THE ARGUMENT

As the Fourier coefficients have been investigated for small 7 and have been shown to be

well-behaved when properly written, what remains is to write these coefficients in a form

which is amenable for investigating their large 7 behavior. As noted before, the hyperbolic

functions, sinh and cosh, grow exponentially for large values of the argument. On the

other hand, the hyperbolic tanh approaches unity for large values of the argument. With

this in mind, let us investigate the form of the Fourier coefficients, written as much as possible in terms of the tanh, which takes care of this potential numerical difficulty. To

this end, it is convenient to introduce the shorthand notation for the hyperbolic functions,

c(x) = cosh(x), s(x) = sinh(x), and t(x) — tanh(z). Making use of this shorthand

notation, the Fourier coefficients may be written as

r3 (n,m,z) = a[Bc(7(L, + z)) + Cs(7(Ls + 2))}, (78)

19
r2 (n, m, *) = a{Dc(j(L3 + L2 + «)) + £*(7(L3 + ^2 + *)) }, (79)

where,

T1 (n,m,z) = As(j{Li + L2 + Li + z), (80)

A

7 t^(7L3) + CC(TL3 )i'

(81)

5 = Dc(iL2 ) + Es{jL2 ), (82)

C= —{^(7X2)4-^(7X2)}, (83)

D = j( 7Ii), (84)

E='^c{1L1 ). (85)

«2

As in the investigation of the small 7 behavior of the Fourier coefficients, the factor

U(n,m)Po I K3 has been deleted from eq (81) for convenience. This factor may simply

be included in the Fourier representation of the temperature, eq (19), as it is common to

all three layers.

Now, the above equations (eqs (78 - 85)) will be rewritten by making use of the definition

of the hyperbolic tanh; i.e., t(x) — s(x)/c(x). First, consider the coefficient C.

C = —{Ds(7L2 ) + Ec(7L2 ) «3 }

C = ^

s(7Li)s(7L2 ) + ~c(7£i)c(7L2 ) I,

«3 I «2 J

c = ^c^l^lA ti-yLiX-YU) + -

«3 I K2 (86)

20
Next, the coefficient B may be written as

B = Dc(jL2 ) + Es(jL2 ),

B = s{ 1L1 )c{lL2) + —0(7^)5(7X2),

«2

B = c(TL1 )c(7L2 )|<(7L1 ) + ^t{jL2 )

Then,

Bs(yL3 ) + Cc(1L3 ) =

c{ 1L1 )c{lL2)s{iLAi{1L1 ) + —*( 7X2 )} + 0(7X00(7X2)0(7X3) —(tfrl^frLa) + —j

- c(7X1 )c(7X2 )C(7L3 )(<(7iiW7^3) + —t(7L3 )t(7L2 ) + ^frXj^X,) + —j. (88)

[ «2 ^3 «3 J

Then, the coefficient A may be written as

A =

1 (89) 70(7X^(7X3)0(7X3)

{

t{lLl)t{lLz) + S,( 7X3)<

( 7X2) + S^WTX,) + g

It is convenient to define the function ^(7) as

0(7) - 7 T~s (90)

|<(7X1 )*(7X3 ) + %t(<yL3 )t(7L2 ) + g*(7X1 )*(7X2 ) + ^}

which is well behaved for all values of 7. Then, the coefficient A may be written as

A= . "' 7 > .. . (91)

70(7X1)0(7X2)0(7X3)

21
By making use of the above procedure, it is relatively straightforward to show that the

Fourier coefficients as given by eqs (78 - 80) may be written as

r3 (n,m,z) =

n(7 )c(7(£. + »)) f +

«, + <( + z))

M /

} + «.) 1 (92)

70(7X3) [ K2

AC3 V «2/ J

0(7)^(7(^3+^2+^)) f v , , r

, , vx\ , Q,v r2 (n,m,z) = < H7^i) + —%(£s +^2 + 2

) f,

93

70(7X3)0(7^2) t ^2 J

and

0(7)^(7(^ + ^2+^+2))

n(n,m,2) = -. (94) 7c(7L3 )c(7L2 )c(7Li)

In eqs (92 - 94), the function 0(7) and the terms inside the curly brackets are well behaved

for all values of the variable 7. The sinh and cosh terms which remain may still give rise to

numerical overflow problems for large values of the argument. However, as both of these

functions grow exponentially for large values of the argument and they appear in both the

numerator and the denominator of the Fourier coefficients, there will be cancellation. This

cancellation for large values of the argument will not be worked out in detail here but is contained in the FORTRAN listing of the program in the function FUNZ.

These Fourier coefficients are used in the equation

„ 4U(n, m^n, m,z)cos(mrxLx )cos(m'Ky/Ly)

Ti(x,y,z) = P0 > > 77——T77-

, 1U T 95)

for i = 1,2,3 to obtain the solutions in each of the three layers. In the above, the term

PoU(n,m)/K3 has been written out explicitly and is no longer contained in the Fourier

coefficients. It is important in obtaining the solution in x,y, z to use the appropriate layer

equation. This is automatically taken care of in the program as the depth z is compared

with the various thicknesses and the corresponding layer Fourier coefficient is used. The

user does not have to specify which layer is to be used.

There are several parts of the above equation which can be identified. First, there is the z-dependent part of the solution in each layer, Tj(n,m,z). The two cosine functions,

cos(mrx /

x ) cos(rmry/

y ), arise from the x- and y-dependent parts of the steady-state

equation. Their specific form arises from the requirement of the lateral boundary conditions. Finally, the requirement that heat enters or leaves the top surface only through the

22
heating elements introduces the function PoU(n,m). It should be noted that Pq appears

only in front of the sum. It is a scale or multiplicative factor which for convenience may

be set equal to unity. Its effect is to scale the temperature in a direct linear fashion.

Performing the calculation for Pq = 1 for a given structure then provides the temperature

distribution for arbitrary Pq which can then be obtained by multiplying by the Pq which

is used.

An interesting exercise left to the reader is to show when the three thermal conductivities are equal that the one-layer solution is obtained. Also, another exercise is to show that

the one-layer solution is obtained when the thicknesses of the second and third layers are

set equal to zero.

1.8 - SPECIAL CASE OF POWER SOURCE COVERING TOP SURFACE

As a special case of eq (95), consider the situation of a single power source completely

covering the top surface; i.e., there is a single heat source with lateral dimensions equal

to that of the three-layer structure. In this particular case, it is shown that the above

equation reduces to the familiar thermal resistance equation. The easiest way to proceed

with the analysis is to consider the specific form of the function U(n,m). From eq (66),

Ui(n, m) =

(96)

For the particular case of uniform surface coverage, X\ = y\ = 0, Ixi — Lx , and ly\ — Ly

. Upon substituting these values into the equation, the function reduces to

U\(n,m) = * y |sin(n7r) sin(m7r) 1. (97)

(n7r)(m7r) ^ J

This is zero when the indices are nonzero. For the case where both of the indices are zero, the use of the expansion of the sine function gives rise to the result that

Ui(n,m) = Lx LySn0 Smo. (98)

If this form of the U(n, m) function is substituted into eq (95) for the temperature in each

of the three layers and the form of the Fourier expansion coefficients as 7 —> 0 (eqs (75 - 77)) is used, it is readily shown that the temperatures in the three layers may be written

as

I*3( W) = Po{^ +

^

+ H (99)

23
L3 + L2 + z Li

T2 {x,y,z) = P0 \ + —

K2 Ki

(100)

rp (

v p f U+L2 +LX +z \ 2l(x,y,z) = P0 | —

J

(101)

As Pq is the power density per unit area, these equations give rise to the usual results of

the one-dimensional calculations of the thermal resistance.

1.9 - EFFECT OF UPPER SUMMATION LIMITS ON TEMPERATURE

It is important to keep in mind that the Fourier coefficients are functions of the variables n

and m. As discussed in the previous sections, the power density function, U(n,m), tends

to zero very slowly for large values of the argument. Also, the cosine terms in eq (95)

do not tend to any limit as the arguments approach infinity. It is the Fourier coefficients

which are responsible for the convergence of the sum. This is especially the case for the

surface (z = 0) value of the temperature. This consideration comes into play when setting

the upper summation limit.

Consider the case of a structure with lateral dimension of L and a heat element of lateral

dimension of A. For the Fourier series to begin "seeing" this element, the cosine function

in the basis set must have at least one complete cycle in A. As the element is rectangular,

one cycle is certainly not enough, as the cosine is a poor representation of the rectangle.

Consequently, higher "frequencies" are required to assure the adequacy of the representation. The rule of thumb is that the number of terms should be at least on the order of

L/A. A stronger rule of thumb would require several times L/A. In general, the stronger

rule of thumb should be applied to the smallest heat element to achieve better accuracy.

1.10 - THE TXYZ AND TXYZ20 CODES

The above equations have been coded in the TXYZ program [6]. Over time, various

updates were included in the TXYZ20 program [7]. In the original version of TXYZ, all heat sources had equal unit weights. The updated TXYZ20 program allows for noninteger

weights to be assigned to each heating element. This weight may be positive or negative

and hence represent the effects of either a heat source or a heat sink.

The reason for the relatively compact nature of the original TXYZ code and its numerical

efficiency is that great care was exercised in the investigation of the functions for the small

and large argument regimes. For small values of the argument, the evaluation of limiting

forms led to expressions which were correct and free of artificial numerical singularities.

For large values of the argument, special care was required. Many of the functions involved

contained the hyperbolic sinh and cosh functions. These approach the exponential function

for large values of the argument. These provide the possibility for numerical infinities or

overflows. Careful investigation of the forms of the functions showed that these numerical

24
infinities were removable by other numerical infinities. For example, an exponential in the

numerator and denominator of a function would give rise to two overflows when evaluated

by a computer but would cancel each other and give rise to a finite result when analytic

evaluation was properly applied. The construction of the function, FUNZ, in the TXYZ

code was based upon such evaluations and their numerical implementations. With the

passage of time and the use of the TXYZ code, one place where such cancellation was

found to be incomplete was in the case of a fairly thick middle layer. For most cases

of interest to the modeling of the steady-state thermal response of semiconductor device

structures, these layers are usually thin. However, this situation has now been carefully

investigated and the numerically stable limiting form has been determined and used in the

code in the middle-layer portion of the FUNZ function. The introduction of the nonsingular

Hmiting form is contained in the TXYZ20 code.

The original TXYZ program and the TXYZ20 update have been used for a number of

investigations. Some of these include: the thermal evaluation of VLSI packages using test chips [8], the understanding of thermal resistance measurements [9], the investigation of

the thermal interaction between electromigation test structures [10], and the benchmarking

of other codes as, for example, in the modeling of MMIC devices for the determination of

the channel temperatures during life tests [11]. The codes have been run on a variety of

machines. They were originally programmed and run on a VAX 11/785. They now are run on PCs and Sun SPARC10 Workstations [12].

PART 2. ANALYTIC EVALUATION OF LINE AND AREA AVERAGES

2.1 - INTRODUCTION

Electrical and optical measurements are often used to determine the local steady-state

temperature of semiconductor device structures. For example, the electrical resistance of a segment of a metallization stripe may be used to measure the average temperature of

that segment of the stripe [13]. Measurements of the temperature of a metal fine on a thin

layer of silicon dioxide deposited on a silicon substrate have been performed previously.

These have been used in conjunction with the calculation of the average temperature of

the metal fine in order to determine the thermal conductivity of the silicon dioxide layer

[14].

Another example is found in an electrical technique for the measurement of the peak junction temperature of power transistors [15]. This has recently been applied to the use of

the gate voltage in the measurement of the average channel temperature of a power GaAs

MESFET (MEtal-Semiconductor Field Effect Transistor) [16,17]. The work in reference

[5] employs the calculated average channel temperature and the measured average channel

temperature to determine a scale factor which is used to extract the peak channel temperature. It is the peak channel temperature which plays a central role in estabhshing

operating conditions for reliable device performance.

The temperature calculations involved in the above determinations of the thermal conductivity of a Si02 film [14] and the peak channel temperature of power GaAs MESFETs [17]

25
make use of the Kokkas model for the steady-state temperature in a rectangular device

structure. In particular, the point function temperature, T(x,t/, z), was calculated by us- ing the numerical implementation of the Kokkas model [4] as presented in the TXYZ codes

[6,7]. However, as indicated in the above discussion, the measured temperature is actually

an average over the area over which the electrical measurement is performed. In order

to make the connection between the measured average temperature and the calculated

average temperature, references [14] and [17] use a set of assumed representative points

at which T(x,y,z) is calculated. These values are then used to construct what might be

called a point-by-point evaluation of the average.

A number of questions need to be addressed in this point-by-point evaluation of the average. First, for a given heat source or active area geometry, how does one choose a set of

representative points? Second, how many points are sufficient to adequately represent the

effects of variations of T(x, y, z) over the area? Third, what are the effects of any residual

oscillations in T(x,y,z) which are inherent in the Fourier analysis involved in the model?

Fortunately, the Kokkas model and the TXYZ numerical implementation provide for an

analytic evaluation of the averages over arbitrary areas as well as fine segments. This gives

the averages directly and thus obviates consideration of questions related to the choice of a set of representative points as well as the number of points. In addition, the calculation

yields a uniformly convergent value of the average, thus bypassing any residual oscillations

of individual values of T(x, y, z). The averages take less time to compute than the pointby-point evaluations and provide a better representation of the measured averages.

2.2 - CALCULATION OF AVERAGE TEMPERATURES

From eq (95), the point function temperature, T{(x,y,z), is

rp / x „ V" \ - ±U{n,m)Ti{n,m,z) cos(nnx/Lx ) cos(miry/Ly ) Ti{x,y,z) = PQ > V 77—7-7777- —777 T

• (102)

{OnO + l)(om0 + l)Lx LyK3 n=0 m =

As indicated above, the point-by-point evaluation of the area average temperature makes

use of a series of point values given by eq (102).

Alternatively, the average temperatures may be calculated directly. The development applies to both line and area averages. Line averages provide a certain amount of information,

but it is the area average which is related to the experimentally measured temperature.

Consider an arbitrary fine segment from (xj) to (xj + Ixj) and an arbitrary rectangular area from (xj,yj) to (xj + lxj,yj + lyj). These are pictured in figure 2 as well as in the

top view shown in figure 3. It is important to note that these are shown on the surface

(z — 0) of the structure for the purposes of clarity and illustration. The arbitrary fine

segment and arbitrary area (in the x-y plane) may be located on or inside the structure

(0 < z < -(Lj + L2 + Ls )).

26

Line for average

TOP VIEW

Figure 3. Top view of the surface of the structure shown in figure 2. The line for averaging

and area for averaging are separated from the heat source for the purpose of clarity.

28
The z- dependence is explicitly expressed in terms of the averages which are defined below.

The average of T{(x, y, z) over the arbitrary line segment is defined as

{Ti(y,z))i = ~ I

3+ 3 Ti(x,y,z)dx, (103)

and the average of T{(x,y,z) over the arbitrary rectangular area

(Ti(z)) a = - — / / Tt (x,y,z)dxdy. (104)

ixj lyj j x . j y

.

Notice that the x and y dependence of Ti(x,y,z) is completely explicit. This means that

the integrals in the line and area averages may be analytically evaluated. It is convenient

to define the functions, Ax (n) and Ay (m), by

, . Lx I . (nitix j+lxj) \ /mrxj\ 1 A' (n) -^(sm

( l, )- sm (-ir)[ < 105>

and

vm,^{s^r^^)-5I„(^)}. doc,

By substituting eq (102) into eqs (103) and (104) and making use of the definitions of the

functions, Ax (n) and Ay (m), given by eqs (105) and (106), it is straightforward to evaluate

the fine and area averages.

Then, the average of Ti(x,y,z) over the arbitrary fine segment is

irnt w -fo 4U(n,rn)Ti(n,m,z)Ax (n)cos(rrnry/L y ) IXj ^-J

^

(6n0 + l)(0m0 + l)Lx Ly K3

and the average of Ti(x, y, z) over the arbitrary rectangular area is

(107)

Pq 4JJ(n,m)Ti(n,m, z)Ax (n)Ay (m)

{ l[z))a=

JX~ly~J Z^

Q l^Q (6nO + l)(6mO+l)Lx Ly K3

' ^

29
Several observations about the line and area average temperature defined by eqs (107) and

(108) are important to note. First, the averages may be calculated by a single evaluation

of the double sum. This means that the averages may be calculated about as quickly as the evaluation of the point value of the temperature, T;(x, y, z). Second, the averages may

include any part or all of the x-y plane defined by the structure. If the area is on the

top surface, it also may contain none, part, or all of the heat source(s). This flexibility

provides for the evaluation of area average temperatures over a portion of the heat source

or over a portion of the surface away from the heat source(s). Third, the averages may be

calculated for the surface (z = 0) as well as any depth (z < 0) inside the structure. Area

averages calculated for different areas and as a function of depth provide information on

the temperature spreading and heat transport inside the structure.

2.3 - NUMERICAL IMPLEMENTATION AND TXYZ30

The basic TXYZ program has been used as the framework in which to construct the

programs to numerically evaluate the equations for the line and area averages in eqs (107)

and (108). The resulting FORTRAN77 code has been named TXYZ30. This code has been

constructed to allow for the calculation of point functions, line averages, or area averages

by means of an integer variable switch which is set in the beginning of the input data file.

It is important to keep in mind that the fines and areas over which the averages are computed are arbitrary. They may be on the top surface (z — 0) or inside the structure

{z < 0). In addition, the averages may be easily computed not only over part or all of

the heat sources but also at locations away from the heat sources. Most applications of

these codes would probably make use of area averaging. Line averaging may be useful in

examining the fine average temperatures. Typical evaluations of the averages take about

the same CPU time as is required to calculate a single point value of T(x,y,z).

The direct calculation of the area average temperature vastly simplifies the analysis of

the data when investigating the thermal conductivity of thin surface layers. The direct

calculation of the averages provides for a faster and more accurate technique for the case of

the determination of the peak channel temperature from the ratio of the measured average

channel temperature and the calculated average channel temperature. This is especially

relevant when there may be a significant variation of T(x,y,z), both across and along the

active area. The point-by-point evaluation of the average channel temperature must take

this into account to provide an adequate representation of the variation to ensure suitable

accuracy. The use of eq (108) is not encumbered by this requirement and provides the

averages directly.

The TXYZ30 code is contained in the first appendix. The FORTRAN source code is also contained in the HOTPAC software package. In addition, there are several sample

files. These are txyz30io.l, txyz30io.2 and txyz30io.3 where the io refers to input/output.

The file extension refers to the example of a point function (1), a line average (2) and an

area average (3). These files contain the input data, the annotated input data and the

corresponding output data for each case. The user can electronically cut out the input

30
data for the input.dat file needed to run TXYZ30.

2.4 - TOTAL AREA AVERAGES AND THE THERMAL RESISTANCE

An interesting example of the use of eq (108) is provided by the situation where the average

is performed over the entire top surface (z = 0) area of the structure. For the case of a

single heat source of dimensions Ix by ly, it is straightforward to show that the average

over the total area (average surface temperature) is (ftP'^Sr. (io9)

which is just the one- dimensional (thermal resistance) result scaled by the ratio of the heat

source area and the total area.

2.5 - APPLICATIONS OF AREA AVERAGES

The evaluation of the averages is simple and direct and easily performed for any number of

arbitrary fine segments or areas. The computed averages over part or all of the heat source area should give a more direct representation of the measured values and hence provide a more accurate determination of the thermal properties as well as peak temperatures. In

addition, the averages are not restricted to the areas where heat enters the structure (over

part or all of the heat source or sources). The averages may also be computed for various

depths inside the structure. Both nonheat source averages and depth averages can provide

information on the effects of structure geometry and thermal properties on the heat flow.

Infrared imaging techniques are also used to experimentally measure the surface temperature of semiconductor device structures. The temperature is an average over the spot size

of the imager. Application of the calculation of the area averages for this situation will

speed the analysis of the average temperatures being measured. This should substantially

improve the understanding and application of the techniques.

Before turning to the thermal multilayer problem, it is important to notice that the application of the averaging technique is general and depends upon the explicit x- and in- dependence of the temperature expressed in eq (102). The thermal multilayer problem will

have exactly the same kind of explicit x- and i/- dependence so that the averaging can be

applied directly. This is discussed in the next part of the report.

PART 3. RECURSION RELATION SOLUTION OF MULTILAYER MODEL

3.1 - INTRODUCTION

It has been two decades since Kokkas' original work on the thermal analysis of multiplelayer structures [3,4]. This is based on the solution of the heat flow equation for a multilayer

rectangular structure. While the interface boundary conditions for the continuity of the

31
temperature and the heat flow are for a general multilayer structure, the equations have

been solved in detail and evaluated for the case of up to a three-layer structure. The

results of the three-layer steady-state heat flow analysis are exact for this situation. This

three-layer solution has been numerically implemented in the TXYZ programs [3,4]. This

model and the codes have found wide use in a variety of steady-state problems.

During this same 20-year time span, there have been a number of advances which have

taken place in the theoretical analysis of the analogous electrical problem of two-probe

spreading resistance. The starting point of the work on the electrical problem is the

multilayer Laplace equation analysis of Schumann and Gardner [18,19]. Originally, the

numerical analysis required a mainframe computer to do the rather cumbersome matrix

algebra which develops quickly once the calculations get past two or three layers. This

limitation coupled with difficulties involved in the evaluation of the accompanying oscillatory integrals served to severely limit the applicability of the technique. However, a major advance came with the introduction of a recursion relation by Choo and coworkers

[20] for the construction of the solution of the iV-layer problem from the solution of the

(N — l)-layer problem. The recursion relation had previously been introduced in the geophysical literature [21,22]. The book by Koefoed [21] provides a comprehensive discussion

of Laplace's equation and its solution, the multilayer approach applied to geoelectric resistance measurements, and the derivation of the recursion relation. It should be noted that

the mathematical description of spreading resistance and geoelectric resistance are exactly

the same even though the length scales are vastly different.

The introduction of the recursion relation eliminated the cumbersome matrix manipulations and allowed for the relatively straightforward algebraic evaluation of the functions

required for the multilayer analysis. The form of the recursion relation has been clarified

by Berkowitz and Lux [23] and has provided one of the cornerstones of the routine calculation of the spreading resistance from the resistivity structure as well as the extraction

of the resistivity structure from the spreading resistance. The remaining cornerstone was

introduced by inventive quadrature techniques which allows for the rapid evaluation of the

necessary integrals. Calculations which would have been prohibitive if not impossible 20

years ago can now be performed routinely on PC- based systems. The model, the numerical

analysis, and the FORTRAN codes have been summarized recently [24].

As the electrical problem and the thermal problem are similar from the mathematical

viewpoint, the question arose as to the applicability of the recursion relation technique,

with possible and appropriate modification, to the steady-state thermal problem. If this

could be accomplished, then it would provide for an exact solution of the steady-state

surface temperature for a multilayer structure with any number of layers. The present

work describes the modified recursion relation and its application to the thermal problem.

The salient features of the recursion relation are presented along with a comparison of its results with those of the Kokkas three-layer solution. The use of the recursion relation

for extension of the calculations beyond three layers is discussed, with emphasis on the

small increase in calculation time for each added layer. The computational efficiency of

the recursion relation is also discussed. This new technique should prove extremely useful

32
in the understanding of the thermal effects of multilayer structures.

As the similarity to the spreading resistance problem is part of the motivation, a brief

discussion of the analysis involved in this problem is presented. This is followed by the

corresponding analysis for the multilayer thermal problem. The connection of these two

problems is then made, along with a discussion of the need for a modified recursion relation

for the thermal case. Its form is then determined and the thermal recursion relation is presented. Finally, the implications of this thermal recursion relation are discussed.

3.2 - THE N-LAYER ELECTRICAL PROBLEM

For the case of the spreading resistance problem, the multilayer solution of the Laplace

equation provides for the calculation of the resistance between the contacts on the top

surface of a nonuniform conductivity (resistivity) structure. The fundamental assumption

of the multilayer analysis is that Laplace's equation is satisfied in each "layer" in the

material. The problem is set up in cylindrical coordinates to emulate the symmetry of a

single circular contact on the top surface of the material. The potential is assumed to be

independent of the angular variable in this coordinate system. The Laplace equation may

then be written as

V

2 V(r, z) = —V(r, z) + - —V(r, z) + —V(r, z) = 0, (110) or1 r or ozL

where V(r, z) is the potential, r is the radial coordinate, and z is the depth coordinate.

This equation may be solved by means of separation of variables, with the result that a

particular solution is

V(r,z) =exp(-Az)J0 (Ar)+exp(+Az)Jo(Ar), (111)

where Jo(Xr) is the Bessel function and A

2 is the separation of variables constant. The

boundary condition on the r part of the solution is that V(r, z) approaches zero as r tends

to infinity. This is satisfied by the above for all values of A. Then, the general solution is an integral of the particular solution with a weighting factor and is of the form

V(r,z) = J{(l + 6(\))exv(-\z) + j>(\)exp(+\z)}j0 (\r)d\, (112)

where the weighting functions, 0(A) and ^(A), are determined from the z-dependent boundary conditions. The above is the general solution of Laplace's equation in cylindrical co- ordinates for a single layer. This provides the framework in which the multilayer solution

may be addressed.

33
vN

dN-i Vn-i

•

•

•

•••

•

-

•

•

d2 V2

di «1 Vi

Figure 4. Schematic representation of the geometry used in the Schumann and Gardner

multilayer Laplace equation analysis. In this figure, di, ai, and Vi are the thickness,

electrical conductivity, and potential, respectively, in the i-ih layer.

34
The depth-dependent portion of the multilayer geometry is presented in figure 4. Each

layer is described by a thickness and electrical conductivity. Charge neutrality is assumed

in each layer so that the potential in each layer satisfies Laplace's equation. The general

one-layer solution given by eq (112) provides a convenient and useful basis for the required

layer solution. Then, the solution of Laplace's equation in the z'-th layer in an TV-layer

structure may be written as

Vi(r,z) = J

{(l+^(A))exp(-Az) + Vz(A)exp(+A^)}j0 (Ar)dA. (113)

The boundary conditions used to solve the system of equations (i.e., determine {0t (X)}

and {V'j(A)}, i = 1, ...,N) are provided by conditions on the top surface, the intermediate

interfaces, and the bottom surface. On the top surface, current flow takes place only

through the probe which is modeled as a circular plate of radius a. Then the top surface

boundary condition is expressed as

dVN(r,z)

-<?N o = J{r): (114) oz

where is the electrical conductivity of the top layer, J(r) is the probe current density

[18,19], and the derivative is evaluated at z = 0. On the bottom surface, the potential is assumed to be well behaved and, more specifically, is assumed to approach zero. This is usually expressed as

lim Vi(r,z) = 0. (115)

Equations (114) and (115) provide 2 of the 2N conditions required to solve the iV-layer

system of equations. The remaining 2(iV — 1) conditions are provided by requirements

at the interfaces between the layers where the potentials and the current densities are assumed to be continuous [18,19]. These are expressed as

Vi (r,z) = Vi -1 (r,z), (116)

and

dV,(r,z) dV.-^z) = (U7)

where the functions and their derivatives are to be evaluated at the interfacial boundaries.

35
For the case of an iV-layer structure, the substitution of eq (113) into the boundary conditions given by eqs (114 - 117) gives rise to a set of 2iV equations in 2N unknowns ({#2 (A)},

{^i(A)}, i — l,...,iV). The analytic solution of this system of equations requires the use

of the Cramer's rule of linear algebra. Clearly, this can become rather tedious especially

since the expansion coefficients are functions of the continuous variable, A. It is possible

to show that the potential on the top surface of the iV-layer structure may be written as

[18,19]

Vn ^°) = 2^J0 A

dK (U8)

where the kernel function, An(\) = 1 + 2#/v(A), depends upon the electrical conductivities

and thicknesses of the "layers" in the multilayer structure (through the solution of the

above system of simultaneous equations). Schumann and Gardner were able to work out

the system of equations for the cases up to i = 3. Cases beyond three layers presented a major stumbling block.

The development of the recursion relation by Koefoed [21] and its use by Choo et aJ.

[20] provided a major breakthrough in the evaluation of eq (118) and effectively removed

the numerical difficulty associated with matrix inversion. The philosophy and utility of a

recursion relation is that the kernel function, A/v(A), for an ./V-layer structure can be easily

generated from the kernel of an (N — 1) layer structure by means of an algebraic relation.

This algebraic relation represents a reduction of the matrix algebra and the subsequent

manipulations with determinants. The reader may wish to refer to Koefoed's book for a

detailed discussion. In particular, if the kernel is known for the (N — 1) layer case, then

the kernel for the N layer case is given by

<T N - 1 u(\) + aNAN- 1 (\) AN {A) = —— —, (119)

(TN-l + <TN«j{X)AN -i(A)

where

=

l-exp(-2A^)

v

' 1 +exp(-2Ac?N )

v '

and d,N is the thickness of layer N. In practice;, the recursion relation is begun with the

one-layer case from which the two-layer kernel is generated. Then the three-layer kernel

is generated from the two-layer kernel. This sequence is repeated until the iV-layer kernel

is determined. Notice how the conductivity and thickness of the iV-th layer enter in the

recursion relation.

Berkowitz and Lux [23] have shown that the recursion relation for the kernel function could

be expressed as

36
Ln

LN-1 KN-1 T|SJ-1

•

•

•

•

•

•

•

•

• K2

Li *1 Ti

Figure 5. Schematic representation of the geometry used in multilayer analysis of the

steady-state heat flow problem. In this figure, Li, K{, and Tt are the thickness, thermal

conductivity, and temperature, respectively, in the 2-th layer.

37
AN (\)

<jn-i tanh(AJ/v) + CAr^TV-i(A)

cryv-i + &N An -i (A) tanh( Xd^

(121)

The recursion relation is most commonly used in this form. This is in part due to the

hyperbolic tangent function which is bounded between (0,1), making the numerical implementation relatively easy. This coupled with the quadrature techniques introduced by

Berkowitz and Lux has moved the analysis into everyday applicability.

3.3 - THE N-LAYER THERMAL PROBLEM

The depth-dependent portion of the multilayer geometry for the thermal problem is shown

in figure 5. Figures 4 and 5 are presented to reinforce the connection between the electrical

and thermal problems and the subsequent development of the thermal recursion relation.

From the development of Part 1, it is possible to show that the temperature on the top

surface of the iV-layer structure may be written as

rr i n\ p 4U(n, m)TN(n, m, 0) cosjnwx Lx ) cos( miry/Ly ) TN {x,y,0) = P0 2^ 71—7~7vx

.i\r r

1 ^ l22 > t^On^O (6n0 +l)(dm0 +l)Lx Ly K NJ

The key elements in eq (122) are the Fourier coefficients, tjv(t), which depend upon the

thermal conductivities and thicknesses of the layers in the multilayer structure (through

the solution of the above system of simultaneous equations).

As discussed previously, the solution of eq (122) for cases beyond three layers presented a major difficulty. The analogous difficulty with the solution of the multilayer Laplace

equation and its removal by means of the recursion relation provides the impetus for the

introduction of a recursive scheme in the solution of multilayer thermal equation.

It is important to note that the factor of 7 has been explicitly written in the denominator of

eq (122). This has been done to simplify the discussion of the recursion relation described

in the next section. This recursion relation will be used to take into account the multilayer

thicknesses and thermal conductivities.

3.4 - THE THERMAL RECURSION RELATION

The strong resemblances of eqs (22) and (114), eqs (24) and (116), and eqs (25) and (117)

lead to the possibility of expressing the Fourier expansion coefficient for the surface temperature of an iV-layer structure in the form of eq (121) with the appropriate transcription

from electrical conductivities to thermal conductivities, A to 7, An(X) to 77^(7), and

to Ln- Then the thermal recursion relation would take the form

38
kn-1 tanh^Lw) + «jyrjy-i(7)

kn-1 + ^NTN-iil) tanh(jLN

)

'

(123)

As the electrical and thermal bottom-layer boundary conditions are of different form, it is necessary to evaluate the corresponding ^(7). The form of ^(7) may be obtained from

eq (92) by setting 2 = 0 and then setting the thickness of the two top layers equal to zero. The result is that

Equations (123) and (124) are the central results of the determination of the thermal

recursion relation. They provide for the determination of the Fourier coefficients of the

surface temperature for an TV-layer structure from the set of thicknesses and thermal conductivities (Ln Kj-, i = 1, TV) by repeated use of the thermal recursion relation. This process

begins with the one-layer Fourier coefficient as given by eq (124). This is substituted into

eq (123) to generate the two-layer Fourier coefficient. The two-layer Fourier coefficient is substituted into eq (123) to determine the three-layer Fourier coefficient. This process is repeated until the Fourier ('06151cient for the desired number of layers is generated. Note

that there is no restriction on the number of layers. This means that the exact steady-state

surface temperature can now be determined for a multilayer structure with an arbitrary

number of layers. It is important to note that this exact steady-state surface temperature

satisfies the heat flow equation as well as the necessary boundary conditions. It does not

need verification with existing layered solutions.

The recursion relation technique should complement and possibly extend some of the other

efforts [25-29] aimed at problems beyond three layers.

It is instructive to compare and connect the results of the thermal recursion relation with

those presented in the literature. The principal one is based upon the case of a threelayer structure where the result of the thermal recursion relation is compared with that

obtained from the Kokkas model and the TXYZ code. Using eqs (123) and (124), it is straightforward to show that the Fourier coefficient for the three-layer structure is given

ti(7) = tanh(7£i). (124)

3.5 - RELATION TO PREVIOUS SOLUTIONS

by

Mi)

K2K3 tanh(7Li) + kiK3 tanh(7X2 ) + «i«2 tanh(7Ls) + k\ tanh(7Li) tanh(7L2 ) tanh(7L3 )

39
1 2 3 4 5 6 7

NUMBER OF LAYERS

Figure 6. Example of the small increase in computation time associated with increasing

the number of layers in the thermal structure. The slope represents about a 5-percent per

layer increase.

40
It is also straightforward to show that this is in exact formal agreement with the Kokkas

and TXYZ expressions for z = 0 for the three-layer case. This does not serve as a limited

verification of the recursion relation but rather shows that the two are in agreement, as expected. In addition to this exact formal agreement, numerical simulations have been

performed with the TXYZ code and a preliminary code based upon the recursion relation.

As expected, these are in numerical agreement. One interesting result is that the recursion

relation code actually runs about 15 percent faster. This is due to the more compact

nature of the recursion relation.

It is interesting to note in eq (125) that if kx = k2 = «3, then ts(j) is independent of

the thermal conductivities and reduces to tanh(7(Xi + L2 + L$)) which is the result for a

single layer of thickness equal to L\ + L2 + L$. This line of reasoning can be extended to

the iV-layer equation. For the case where the thermal conductivities of the layers become

equal, repeated use of the addition formula for the hyperbolic tangent shows that the form

of the Fourier coefficient becomes tanh{j(Li + —h Ln)) as expected for a single layer of

thickness equal to L\ + • • • -f Ln3.6 - NUMERICAL IMPLEMENTATION OF THERMAL MULTILAYER CODE, TML

The numerical implementation of the thermal recursion relation solution is contained in

the Thermal MultiLayer, TML, code. The FORTRAN source code is listed in the appendix

and is also contained in the HOTPAC software package.

For the sake of illustrating some of the features of the code, the TML program has been used

in calculations of the surface temperature for structures of up to seven layers. The case of seven layers is used to provide a convenient range and does not represent a limitation of the

recursion relation. Any number of layers is possible and can be easily achieved on modern

computing systems. These temperature calculations were performed for the purpose of

testing the computational efficiency of the recursive scheme. The results are presented

in figure 6 and were obtained on a VAX 11/785. They are shown for the purpose of

illustration of the utility of the recursion relation. There is only about a 5 percent increase

in computation time for each layer added to the structure. This is strong evidence for the

effectiveness of the algebraic nature of the recursion relation.

These test cases were also run for the situation where all the thermal conductivities were

equal. In keeping with the above discussion, the results reduced to those for a single layer

of thickness equal to L\ -\ VL^. This provided further support of the recursion relation

and its numerical implementation.

For the purposes of acquainting the user with the TML code, a number of examples are contained in the files. The particular files, tmlio.l, tmlio.2, tmlio.3, contain the three-layer

problem contained in the corresponding TXYZ30 related files, txyz30io.l, txyz30io.2, and

txyz30io.3. The TML code, like the TXYZ30 code, contains the point function, line average

and area average options. These are illustrated by the above.

For the purpose of showing the use of TML for real multilayer cases, the files tmlio.120,

41
tmlio.220, and tmlio.320 are included in the software package. These illustrate the case

of the TML calculation for a 20-layer case in the point function, line average and area average modes.

3.7 - GENERALIZED ONE-DIMENSIONAL THERMAL RESISTANCE

The special case of uniform surface coverage by a single heat source for a three-layer

structure has been presented in section 1.8 of this report. There it was shown that only

the 7 —> 0(n = 0, m = 0) term needs to be considered. This same line of reasoning applies

to the full TV-layer case. Investigation of the recursion relation for small values of 7 shows

for the case of uniform surface coverage by a single heat source that

N

T 2^(0) = P0 £ (126)

which is just the generalized one-dimensional thermal resistance result.

In section 2.4, the average surface temperature (over the entire top surface (z = 0) area

of the structure) was calculated for the case of a single heat source of dimensions Ix by

ly. The result is presented in eq (109) and is just the one-dimensional thermal resistance

scaled by the ratio of the heat source area and the total area. Using the analysis used to

obtain eq (126), it is easy to show that the average surface temperature for the N-layer

structure is

,TM / ^ Ix ly v^. L,

(TN (0)) a =Po-r-rY,-'

1

y t=i

which is just the generalized one-dimensional thermal resistance scaled by the ratio of the

heat source area and the total area.

3.8 - APPLICATION TO BURIED OXIDE STRUCTURE

The thermal recursion relation should be especially useful in the understanding of the

effects of multilayer structure thickness and thermal conductivities on the surface temperature. The possibility of using the surface temperature as a probe of the thermal

conductivities of the structure is particularly intriguing. An example of the use of this

technique is found in the case of the steady-state surface temperature of a buried oxide

structure. Recent trends in semiconductor fabrication make use of the buried oxide as a way of electrically isolating the top active region. Probably the most popular version of

this is found in SIMOX (Separation by IMplanted OXygen) structures where the buried

oxide is formed by high-dose implantation of oxygen followed by high-temperature annealing [30]. The depth of the buried oxide is controlled by the energy of the incident oxygen

ion beam.

(127)
0 200 400 600 800 1000

INVERSE THERMAL CONDUCTIVITY

OF Si02 (cm K / W)

Figure 7. Example of the application of the thermal recursion relation to the calculation

of the surface temperature of an SOI (Silicon-On-Insulator) structure. The insert (upper

left) depicts the structure where the layer thicknesses and heater dimensions are not to

scale. The buried SiC>2 is 0.00001 cm thick. Curves A, B, and C are for the cases where the

surface layers of Si are 0.00001 cm, 0.00002 cm, and 0.00003 cm thick, respectively. Notice

how the temperatures do not follow a one-dimensional interpretation in keeping with the

two-dimensional heat flow in the more thermally conductive surface silicon layer.

43
When processed correctly, the implantation damage should be removed and a buried layer

of SiO-2 should be formed. Electrical isolation is certainly achieved but the thermal conductivity of the buried oxide may impose a price of thermal isolation which may give rise

to higher-than-desired operating temperatures in the device. These kinds of structures are prime candidates for the investigation of self-heating effects in multiple-layered structures

[31,32].

As an example of the recursion relation technique, preliminary calculations were performed

on a series of buried oxide structures. These were done to ascertain the possibility of using

surface temperature measurements as a probe of the thermal conductivity of the buried

oxide. The basic structure was modeled as a chip (rectangular structure) with lateral dimensions of 1.0275 cm by 0.7737 cm. The vertical structure consisted of a surface layer of

Si over a buried layer of SiC>2 over a substrate layer of Si. The buried oxide was 0.00001 cm

thick and the substrate Si layer was 0.03189 cm thick. Three thicknesses of the surface Si

layer were used in the calculations. These included 0.00001 cm, 0.00002 cm, and 0.00003

cm, and are typical of SIMOX structures. The thermal conductivities of the surface and

substrate Si layers were taken as the room temperature value of 1.55 W/cm°K and the

thermal conductivities of the buried oxide were taken in the range from 0.001 to 0.015

W/cm°K . A small resistive heat source on the top surface was modeled as a rectangular

element of lateral dimension of 0.001005 cm by 0.0825 cm. The temperature (normalized

by the power density, T/ Pq) was calculated in the area of the heat source for the ranges of

surface Si thickness and buried oxide thermal conductivities. This was done to ascertain:

1) the effects of the surface silicon thickness and 2) the sensitivity of the calculated temperature on the buried oxide thermal conductivity. The results are presented in figure 7 where

the calculated surface temperatures (normalized by Pq) in the area of the small surface

heat source are plotted against the inverse thermal conductivity of the buried oxide for

the three surface silicon thicknesses. These results indicate the important dependence on

the thickness of the surface silicon layer. In addition, the calculated temperature decreases

with increasing thickness of the surface silicon layer. If the heat flow were one-dimensional

and described by eq (126), then the calculated temperature should increase with increasing

thickness of the surface silicon layer. The fact that the trend is in the opposite direction

points to the two-dimensional heat flow caused by the more thermally conductive surface

silicon layer. This necessitates careful interpretation of the temperature data on these

structures.

3.9 - SUMMARY AND CONCLUSIONS

Using the strong mathematical similarity of the multilayer Laplace equation of spreading

resistance and the multilayer steady-state heat flow equation, a recursion relation has been

developed for the multilayer thermal problem. There are no restrictions on the number of

layers which can be used in the calculations. This makes the results truly multilayer. This

recursion relation has been shown to produce the Kokkas-TXYZ results for the three-layer

case. It also gives rise to a generalized one-dimensional thermal resistance result for the

case of uniform surface coverage. A preliminary program based on this recursive scheme

has been shown to provide very good computational speed. As an example, it has been

44
applied to investigate the possible determination of the buried oxide thermal conductivity

in SIMOX-type structures. This recursion relation technique is simple, elegant, and powerful and should be extremely useful in the investigation and understanding of the effects

of layer thicknesses and thermal conductivities on the steady-state surface temperatures

of multilayer structures [33, 34].

AVAILABILITY OF HOTPAC SOFTWARE PACKAGE

The source codes for the HOTPAC software package are written in FORTRAN77 and have

been run on a VAX 11/785 minicomputer system as well as a Sun SPARC10 workstation.

These codes should run on other systems with little, if any, need for modification.

There are two programs contained in the HOTPAC software package. The FORTRAN77

source code and sample input and output data files are available in ASCII format on DOSformatted floppy disks. This package is self-contained and is straightforward to run once

the FORTRAN is compiled and linked by the user-supplied software. The sample input

and output data files are included so that the user can check the programs for proper

operation as well as to become acquainted with the setup and use of the codes.

For more information or to receive a copy of the code and the report, please contact:

John Albers

NIST, Bldg. 225, Room B-310

Gaithersburg, MD USA 20899

Telephone: 1-301-975-2075

FAX Number: 1-301-948-4081

Netmail (Internet) address: Albers@sed.eeel.nist.gov

45
ACKNOWLEDGMENTS

Over the years, the author has had the opportunity and privilege of communicating and/or

collaborating with a number of people in the fields of thermal modeling, numerical simulation and analysis. The author would like to thank all these people for their interest and

input:

Frank F. Oettinger for originally suggesting the work on the Kokkas model and for his

continued interest during the course of that and subsequent work. His relentless use of the

original TXYZ program during its developmental stage brought many of the early bugs to

the surface.

Stephen E. Ross for changing the original TXYZ code into a portable FORTRAN form.

His assistance in getting the three-dimensional plots for the original TXYZ report is also

appreciated.

Harry Schafft for his continued use of the code and the collaboration on the electromigration

problem which fostered the variable weight idea and the use of the code for sinks as well as sources. Subsequent work with Harry on the problem of the determination of the thermal

conductivity of thin, surface layers of Si02 led to the idea of the analytical evalauation of

line and area averages.

Dave Blackburn for using the code for thick middle layers which showed the need to

determine a numerically stable form for the FUNZ function for this situation.

Harry L. Berkowitz, previously of Fort Monmouth and Solid State Measurements, for a number of past collaborations on the multilayer Laplace equation description of electrical

spreading resistance. These ultimately led to the investigation of the recursion relation

technique for the thermal multilayer problem.

Jim Wilson and colleagues at Texas Instruments for showing how the TXYZ code could

be used as an exact benchmark for other programs.

AU of the TXYZ/TXYZ20 users for making use of the code. He hopes that TXYZ/TXYZ20

users become HOTPAC (TXYZ30/TML) users.

Finally, the author wishes to thank E. Jane Walters for her care and patience in the

preparation of the report.

This work was funded by the National Institute of Standards and Technology's Advanced

Technology Program in support of the American Scaled-Electronics Consortium.

46
REFERENCES

H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Second Edition (Oxford

University Press, 1959).

R. V. Churchill, Fourier Series and Boundary Value Problems, Second Edition (McGrawHill, 1963).

A. G. Kokkas, Analysis and Design of Electrothermal Integrated Circuits, Ph.D. Thesis, MIT, 1972.

A. G. Kokkas, Thermal Analysis of Multiple-Layer Structures, IEEE Trans. Electron

Devices ED-21, pp. 674-681, 1974.

F. Ayres Jr., Matrices (Schaum Publishing Co., 1950).

J. Albers, Semiconductor Measurement Technology: TXYZ: A Program for Semiconductor IC Thermal Analysis, NBS Special Publication 400-79, 1984.

J. Albers, Semiconductor Measurement Technology: Version 2.0 of the TXYZ Thermal

Analysis Program: TXYZ20, NIST Special Publication 400-89, 1992.

F. F. Oettinger, Thermal Evaluation of VLSI Packages Using Test Chips - A Critical

Review, SoUd State Technology (February 1984).

F. F. Oettinger and D. L. Blackburn, Semiconductor Measurement Technology: Thermal Resistance Measurements, NIST Special Publication 400-86 (July 1990).

H. A. Schafft and J. Albers, Thermal Interactions Between Electromigration Test

Structures, IEEE Proceedings on Microelectronic Test Structures, Vol. 1, No. 1,

pp. 132-137 (February 1988).

J. L. Wright, B. W. Marks and K. D. Decker, Modeling of MMIC Devices for Determining MMIC Channel Temperatures During Life Tests, Proceedings of the Seventh

IEEE SEMI-THERM Symposium, pp. 131-139 (1991).

Certain commercial equipment, instruments, or materials are identified in this paper

in order to specify the procedure adequately. Such identification does not imply rec- ommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best

available for the purpose. In spite of the author's experience that the programs performs correctly with every set of data which has been tried, there can be no assurance

that the program will perform equally well on all (possibly anomalous) data. Therefore, both the author and NIST assume no liability for any possible losses resulting

from the use of this program.

[13] H. A. Schafft and J. S. Suehle, The Measurement, Use and Interpretation of the

47
Temperature Coefficient of Resistance of Metallizations, Solid State Electronics 35,

No. 3, pp. 403-410 (1992).

[14] H. A. Schafft, J. S. Suehle, and P. G. Mirel, Thermal Conductivity Measurements of

Thin-Film Silicon Dioxide, Proceedings IEEE International Conference on Microelectronic Test Structures 2, No. 1, pp. 121-125 (March 1989).

[15] D. L. Blackburn, An Electrical Technique for the Measurement of the Peak Junction Temperature of Power Transistors, Proceedings of the 13th IEEE International

RehabUity Physics Symposium, pp. 142-150 (1975).

[16] C. Z. Lu, M. Z. Wang, X. Gui, G. B. Gao, and B. J. Yi, Fast Measurement of the Peak

Junction Temperature of Power Transistors Using Electrical Method, Proceedings of

the Seventh IEEE SEMI-THERM Symposium, pp. 39-43 (1991).

[17] M. Z. Wang, C. Z. Lu, Z. H. Cheng, Z. Wang, S. W. Feng, G. Y. Ding, X. X. Li, and

G. B. Gao, Electrically Measuring the Peak Channel Temperature of Power GaAs

MESFET, Proceedings of the Ninth IEEE SEMI-THERM Symposium, pp. 112-116

(1993).

[18] P. A. Schumann and E. E. Gardner, Application of Multilayer Potential Distribution

to Spreading Resistance Correction Factors, Journal of the Electrochemical Society

116, pp. 87-91, 1969.

[19] E. E. Gardner and P. A. Schumann, Spreading Resistance Correction Factors, SohdState Electronics 12, pp. 371-375, 1969.

[20] S. C. Choo, M. S. Leong, and K. L. Kuan, On the Calculation of Spreading Resistance

Correction Factors, Sohd-State Electronics 19, pp. 561-565, 1976.

[21] O. Koefoed, The Apphcation of the Kernel Function in Interpretation of Geoelectrical

Resistivity Measurements, Berlin- Stuttgart, Gebriider Borntraeger, 1968.

[22] E. D. Sunde, Earth Conduction Effects in Transmission Systems, New York, Dover,

1968.

[23] H. L. Berkowitz and R. A. Lux, An Efficient Integration Technique for Use in the

Multilayer Analysis of Spreading Resistance Profiles, Journal of the Electrochemical

Society 128, pp. 1137-1141, 1981.

[24] J. Albers and H. L. Berkowitz, Semiconductor Measurement Technology: A Collection

of Computer Programs for Two-Probe Resistance (Spreading Resistance) and FourProbe Resistance Calculations, RESPAC, NIST Special Publication 400-91, 1993.

[25] G. N. Ellison, Thermal Computations for Electronic Equipment New York, Van Nostrand Reinhold, 1984.

48
[26] G. N. Ellison. TAMS - A Thermal Analyzer for Multilayer Structures, ElectroSoft 1,

pp. 85-97. 1990.

[27] G. N. Ellison, Extensions of a Closed Form Method for Substrate Thermal Analyzers to Include Thermal Resistances from Source-to-Substrate and Source-to- Ambient,

Proceedings of the Seventh IEEE SEMI-THERM Symposium, pp. 140-145, 1991.

[28] C. C. Lee, J. Min, and A. L. Palisoc, An Integration Algorithm for the Temperature

Solution of the Four-Layer Infinite Plate Structure, Proceedings of the Fifth IEEE

SEMI-THERM Symposium, pp. 48-51. 1989.

[29] C. C. Lee and D. H. Chien. Thermal and Package Design of High Power Laser Diodes,

Proceedings of the Ninth IEEE SEMI-THERM Symposium, pp. 75-80, 1993.

[30] G. K. Celler and A. E. White. Buried Oxide and Suicide Formation by High-Dose

Implantation in Silicon, Materials Research Society Bulletin, pp. 40-46, 1992.

[31] L. T. Su. K. E. Goodson, D. A. Antoniadis, M. I. Flik, and J. E. Chung, Measurements

and Modeling of Self-Heating Effects in SOI NMOSFETS. Proceedings of the IEEE

International Electron Devices Meeting, pp. 357-360, 1992.

[32] P. R. Ganci. J-J. J. Hajjar. T. Clark, P. Humphries, J. Lapham, and D. Buss, SelfHeating in High Performance Bipolar Transistors Fabricated on SOI Substrates, Proceedings of the IEEE International Electron Devices Meeting, pp. 417-420, 1992.

[33] J. Albers. An Exact Solution of the Steady-State Surface Temperature for a General

Multilayer Structure. Proceedings of the Tenth IEEE SEMI-THERM Symposium,

pp. 129-137. 1994.

[34] J. Albers. An Exact Recursion Relation Solution for the Steady-State Surface Temperature of a General Multilayer Structure, IEEE Transactions on Components, Packaging and Manufacturing Technology - Part A, 18, No. 1, pp. 31-38, 1995.

49
TXYZ30 I/O FILE LISTING - txyz30io.l

110 10

1 .1 .5 1.00

.4 .5 50 50

111 0.0 1.

1

1 5. 0.0

1

1 0.0 0.0

1 1.0

1 4.0 2.0 4.0 2.0

The above is the input used in the input.dat file for the example of a txyz30 point function calculation. This input is annotated below.

INPUT DESCRIPTION OF INPUT

1 itype (=1 for point function calculation)

10 10 x and y dimensions of rectangular structure

1 0.1 thickness and thermal conductivity of top layer (3) .5 1.00 thickness and thermal conductivity of middle layer

.4 .5 thickness and thermal conductivity of bottom layer

50 50 upper summation limits for n and m summations

1 iedgex (=1, then read in three values on next line)

11 0.0 1. number of values, first point, increment

1 iedgey (=1, then read in three values on next line)

1 5. 0.0 number of values, first point, increment

1 iedgez (=1, then read in three values on next line)

1 0.0 0.0 number of values, first point, increment

1 1.0 number of sources and power density

1 4.0 2.0 4.0 2.0 weight, x, length along x, y, length along y

Below is the output contained in the output.dat file which is calculated

and written by txyz30 using the above input.

STEADY-STATE THERMAL CALCULATION KOKKAS ANALYSIS

POINT FUNCTION EVALUATION OF T (X, Y, Z)

X Y z T (X, Y, Z)

0 00000 5. 00000 0 00000 0. 02264

1 00000 5. 00000 0 00000 0. 03674

2 00000 5. 00000 0 00000 0. 15234

3 00000 5. 00000 0 00000 0. 65010

4 00000 5. 00000 0 00000 4. 30401

5 00000 5. 00000 0 00000 7 . 56167

6 00000 5. 00000 0 00000 4. 30401

7 00000 5. 00000 0 00000 0. 65010

8 00000 5. 00000 0 00000 0. 15234

9 00000 5. 00000 0 00000 0. 03674

10 00000 5. 00000 0 00000 0. 02264

10.00 LY= 10.00

THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L3= 1.00000 L2= 0.50000 Ll= 0.40000

K3= 0.10000 K2= 1.00000 Kl= 0.50000

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

NUMBER OF HEAT SOURCE S= 1 POWER DENSITY= 1.000000

WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

50
TXYZ30 I/O FILE LISTING - txyz30io.2

210 10

1 .1

.5 1.00

.4 .5

50 50

111 0.0 10.

1

1 5. 0.0

1

1 0.0 0.0

1 1.0

1 4.0 2.0 4.0 2.0

14.0 2.0

The above is the input used in the input.dat file for the example of a txyz30 line average calculation. This input is annotated below.

INPUT

2

10 10

1 .1

.5 1.0

.4 .5

50 50

111 0.0 1

1

1 50. 0. 0

1

1 0.0 0. 0

1 1.0

1 4.0 2. 0 14.0 2. 0

DESCRIPTION OF INPUT

itype (=2 for line average calculation)

x and y dimensions of rectangular structure

thickness and thermal conductivity of top layer (3)

thickness and thermal conductivity of middle layer (2)

thickness and thermal conductivity of bottom layer (1)

upper summation limits for n and m summations

iedgex (=1, then read in three values on next line)

number of values, first point, increment

iedgey (=1, then read in three values on next line)

number of values, first point, increment

iedgez (=1, then read in three values on next line)

number of values, first point, increment

number of sources and power density

4.0 2.0 weight, x, length along x, y, length along y

nline (=1 for one line) as itype=2

x, length along x for line element

Below is the output contained in the output.dat file which is calculated

and written by txyz30 using the above input.

STEADY-STATE THERMAL CALCULATION KOKKAS ANALYSIS

LINE AVERAGE EVALUATION OF TEMPERATURE

LINE # XLINE LXLINE Y Z AVE TEMP

1 4.00000 2.00000 5.00000 0.00000 6.71515

LX= 10.00 LY= 10.00

THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L3= 1.00000 L2= 0.50000 Ll= 0.40000

K3= 0.10000 K2= 1.00000 Kl= 0.50000

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

NUMBER OF HEAT SOURCES= 1 POWER DENSITY= 1.000000

WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

NUMBER OF LINES= 1

1 4.00000 2.00000

51
TXYZ30 I/O FILE LISTING - txyz30io.3

310 10

1 .1

.5 1.00

.4 .5

50 50

111 0.0 10.

1

1 5. 0.0

1

1 0.0 0.0

1 1.0

1 4.0 2.0 4.0 2.0

14.0 2.0 4.0 2.0

The above is the input used in the input.dat file for the example of a txyz30 area average calculation. This input is annotated below.

INPUT DESCRIPTION OF INPUT

3 itype (=3 for area average calculation)

10 10 x and y dimensions of rectangular structure

1 .1 thickness and thermal conductivity of top layer (3)

.5 1.00 thickness and thermal conductivity of middle layer (2)

.4 .5 thickness and thermal conductivity of bottom layer (1)

50 50 upper summation limits for n and m summations

1 iedgex (=1, then read in three values on next line)

11 0.0 1. number of values, first point, increment

1 iedgey (=1, then read in three values on next line)

1 5. 0.0 number of values, first point, increment

1 iedgez (=1, then read in three values on next line)

1 0.0 0.0 number of values, first point, increment

1 1.0 number of sources and power density

1 4.0 2.0 4.0 2.0 weight, x, length along x, y, length along y

1 narea (=1 for one area) as itype=3

4.0 2.0 4.0 2.0 x, length along x, y, length along y for area

Below is the output contained in the output.dat file which is calculated

and written by txyz30 using the above input.

STEADY-STATE THERMAL CALCULATION KOKKAS ANALYSIS

AREA AVERAGE EVALUATION OF TEMPERATURE

AREA # XAREA LXAREA YAREA LYAREA Z AVE TEMP

1 4.00000 2.00000 4.00000 2.00000 0.00000 5.97551

LX= 10.00 LY= 10.00

THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L3= 1.00000 L2= 0.50000 Ll= 0.40000

K3= 0.10000 K2= 1.00000 Kl= 0.50000

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

NUMBER OF HEAT SOURCE S= 1 POWER DENSITY= 1.000000

WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

NUMBER OF AREAS= 1 1 4.00000 2.00000 4.00000 2.00000

52
TML I/O FILE LISTING - tmlio.l

110 10

3

1 .1

.5 1.00

.4 .5 50 50

111 0.0 1.

11 5. 0.0

1 1.0

1 4.0 2.0 4.0 2.0

The above is the input used in the input.dat file for the example of a tml point function calculation. This input is annotated below.

INPUT DESCRIPTION OF INPUT

1 itype (=1 for point function calculation)

10 10 x and y dimensions of rectangular structure

3 number of layers in the structure

1 .1 thickness and thermal conductivity of layer 3 .5 1.00 thickness and thermal conductivity of layer 2 .4 .5 thickness and thermal conductivity of layer 1 50 50 upper summation limits for n and m summations

1 iedgex (=1, then read in three values on next line)

11 0.0 1. number of values, first point, increment

1 iedgey (=1, then read in three values on next line)

1 5. 0.0 number of values, first point, increment

1 1.0 number of sources and power density

1 4.0 2.0 4.0 2.0 weight, x, length along x, y, length along y

Below is the output contained in the output.dat file which is calculated

and written by tml using the above input.

STEADY-STATE THERMAL MULTILAYER CALCULATION

USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS

POINT FUNCTION EVALUATION OF SURFACE T(X,Y)

X Y T(X,Y)

0.00000 5.00000 0.02264

1.00000 5.00000 0.03674

2.00000 5.00000 0.15234

3.00000 5.00000 0.65010

4.00000 5.00000 4.30401

5.00000 5.00000 7.56167

6.00000 5.00000 4.30400

7.00000 5.00000 0.65010

8.00000 5.00000 0.15234

9.00000 5.00000 0.03674

10.00000 5.00000 0.02264

LX= 10.00 LY= 10.00

NUMBER OF LAYERS IN STRUCTURE= 3 THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L 3= 1.00000 K 3= 0.10000

L 2= 0.50000 K 2= 1.00000

L 1= 0.40000 K 1= 0.50000

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

POWER DENSITY= 1.000000

NUMBER OF HEAT SOURCE S= 1 WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

53
TML I/O FILE LISTING - tmlio.2

210 10

3

1 .1

.5 1.00

.4 .5

50 50

111 0.0 1.

1

1 5. 0.0

1 1.0

1 4.0 2.0 4.0 2.0

14.0 2.0

The above is the input used in the input.dat file for the example of a tml line average calculation. This input is annotated below.

INPUT DESCRIPTION OF INPUT

2 itype (=2 for line average calculation)

10 10 x and y dimensions of rectangular structure

3 number of layers in the structure

1 .1 thickness and thermal conductivity of layer 3 .5 1.00 thickness and thermal conductivity of layer 2

.4 .5 thickness and thermal conductivity of layer 1 50 50 upper summation limits for n and m summations

1 iedgex (=1, then read in three values on next line)

11 0.0 1. number of values, first point, increment

1 iedgey (=1, then read in three values on next line)

1 5. 0.0 number of values, first point, increment

1 1.0 number of sources and power density

1 4.0 2.0 4.0 2.0 weight, x, length along x, y, length along y

1 nline (=1 for one line) as itype=2

4.0 2.0 x, length along x for line element

Below is the output contained in the output.dat file which is calculated

and written by tml using the above input.

STEADY-STATE THERMAL MULTILAYER CALCULATION

USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS

LINE AVERAGE EVALUATION OF SURFACE TEMPERATURE

LINE # XLINE LXLINE Y AVE TEMP

1 4.00000 2.00000 5.00000 6.71515

LX= 10.00 LY= 10.00

NUMBER OF LAYERS IN STRUCTURE= 3 THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L 3= 1.00000 K 3= 0.10000

L 2= 0.50000 K 2= 1.00000

L 1= 0.40000 K 1= 0.50000

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

POWER DENSITY= 1.000000

NUMBER OF HEAT SOURCES= 1 WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

NUMBER OF LINES= 1

1 4.00000 2.00000

54
TML I/O FILE LISTING - tmlio.3

3

10 10

3

1 .1

.5 1.00

.4 .5

50 50

111 0.0 1.

11 5. 0.0

1 1.0

1 4.0 2.0 4.0 2.0

14.0 2.0 4.0 2.0

The above is the input used in the input.dat file for the example of a tml area average calculation. This input is annotated below.

INPUT DESCRIPTION OF INPUT

3 itype (=3 for area average calculation)

10 10 x and y dimensions of rectangular structure

3 number of layers in structure

1 . 1 thickness and thermal conductivity of layer 3 .5 1.00 thickness and thermal conductivity of layer 2 .4 .5 thickness and thermal conductivity of layer 1 50 50 upper summation limits for n and m summations

1 iedgex (=1, then read in three values on next line)

11 0.0 1. number of values, first point, increment

1 iedgey (=1, then read in three values on next line)

1 5. 0.0 number of values, first point, increment

1 1.0 number of sources and power density

1 4.0 2.0 4.0 2.0 weight, x, length along x, y, length along y

1 narea (=1 for one area) as itype=3

4.0 2.0 4.0 2.0 x, length along x, y, length along y for area

Below is the output contained in the output.dat file which is calculated

and written by tml using the above input.

STEADY-STATE THERMAL MULTILAYER CALCULATION

USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS

AREA AVERAGE EVALUATION OF SURFACE TEMPERATURE

AREA # XAREA LXAREA YAREA LYAREA AVE TEMP

1 4.00000 2.00000 4.00000 2.00000 5.97551

LX= 10.00 LY= 10.00

NUMBER OF LAYERS IN STRUCTURE= 3 THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L 3= 1.00000 K 3= 0.10000

L 2= 0.50000 K 2= 1.00000

L 1= 0.40000 K 1= 0.50000

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

POWER DENSITY= 1.000000

NUMBER OF HEAT SOURCES= 1 WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.000.00 4.00000 4.00000 2.00000 2.00000

NUMBER OF AREAS= 1

1 4.00000 2.00000 4.00000 2.00000

55
TML I/O FILE LISTING - tmlio.120

110 10

20

i a a a a a

. lUUUUU 1 A A A A A . lUUUUU

i a a a a a

. lUUUUU T T O O A O

i a a a a a

. lUUUUU 1 O C O Q ^

i a a a a a

. 1UUU UU , 141z D4

i a a a a a .100000 T C O A O A .100000 .177828

.100000 . 199526

.100000 . 223872

.100000 .251189

.100000 OAT O O O .281838

-1 r\ A A A A

. 100000 O 1 £ o o o

. 316228

.100000 .354813

i r\ r\ r\ r\ r\

. 100000 . i ybio / .100000 .446683

.100000 .501187

.100000 .562341

i r\n Ann B lUUUUU . Do u y O /

1 A A A A A - lUUUUU . / U / y 1 3

1 /"\ A A A A .100000 n n a o o o .100000 .891250

.0 1.

0.0

.0

0 2.0 4 0 2.0

The above is the input used in the input.dat file for the example of a tml point function calculation for a 20 layer structure.

This input is annotated below.

INPUT DESCRIPTION OF INPUT

1 itype (=1 for point function calcu

10 10 x and y dimensions of rectangular

20 number of layers in the structure

.100000 .100000 thickness and thermal conductivity

.100000 .112202 thickness and thermal conductivity

.100000 .125893 thickness and thermal conductivity

.100000 .141254 thickness and thermal conductivity

.100000 .158489 thickness and thermal conductivity

.100000 .177828 thickness and thermal conductivity

.100000 .19952 6 thickness and thermal conductivity

.100000 .223872 thickness and thermal conductivity

.100000 .251189 thickness and thermal conductivity

.100000 .281838 thickness and thermal conductivity

.100000 .316228 thickness and thermal

.100000 .354813 thickness and thermal

.100000 .398107 thickness and thermal conductivity

.100000 .446683 thickness and thermal conductivity

.100000 .501187 thickness and thermal conductivity

.100000 .562341 thickness and thermal conductivity

.100000 .630957 thickness and thermal conductivity

.100000 .7 07 945 thickness and thermal conductivity

.100000 .7 94328 thickness and thermal conductivity

.100000 .891250 thickness and thermal conductivity

lation)

structure

conduct ivit

y

conductivity

of layer 20

of layer 19

of layer 18

of layer 17

of layer 16

of layer 15

of layer 14

of layer 13

of layer 12

of layer 11

of layer 10

of layer 9 of layer 8 of layer 7 of layer 6 of layer 5 of layer 4 of layer 3

of layer 2 of layer 1

56
TML I/O FILE LISTING - tmlio.120

50 50 upper summation limits for n and m summations

I iedgex (=1, then read in three values on next line)

II 0.0 1. number of values, first point, increment

1 iedgey (=1, then read in three values on next line)

1 5. 0.0 number of values, first point, increment

1 1.0 number of sources and power density

1 4.0 2.0 4.0 2.0 weight, x, length along x, y, length along y

Below is the output contained in the output.dat file which is calculated

and written by tml using the above input for a 20 layer structure.

STEADY-STATE THERMAL MULTILAYER CALCULATION

USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS

POINT FUNCTION EVALUATION OF SURFACE T (X, Y)

X Y I (X,Y)

0.,00000 5 .00000 0. 01703

1.,00000 5 .00000 0. 02522

2.,00000 5 .00000 0. 10517

3.,00000 5 .00000 0. 42075

4..00000 5 .00000 3. 24066

5..00000 5 .00000 5. 81674

6.. 00000 5 . 00000 3. 24066

7..00000 5 .00000 0 . 42075

3,.00000 5 . 00000 0. 10517

9..00000 5 .00000 0. 02522

10,.00000 5 .00000 0. 01703

10.00 LY= 10.00

NUMBER OF LAYERS IN STRUCTURE=20

THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L20= 0 .10000 K20= 0.,10000

L19= 0 . 10000 K19= 0.,11220

L18= 0 .10000 K18= 0.,12589

L17 = 0 .10000 K17 = 0.,14125

LI 6= 0 .10000 K16= 0.,15849

L15= 0 . 10000 K15= 0 .. 17783

L14= 0 . 10000 K14= 0 , .19953

L13= 0 . 10000 K13= 0..22387

L12= 0 . 10000 K12= 0..25119

Lll = 0 . 10000 Kll = 0..28184

L10 = 0 . 10000 K10= 0 . .31623

L 9= 0 .10000 K 9= 0..35481

L 8= 0 .10000 K 8= 0..39811

L 7 = 0 . 10000 K 7 = 0..44668

L 6= 0 .10000 K 6= 0..50119

L 5= 0 . 10000 K 5= 0..56234

L 4= 0 . 10000 K 4= 0 . 63096

L 3= 0 .10000 K 3= 0 .70794

L 2= 0 .10000 K 2= 0 .79433

L 1= 0 .10000 K 1 = 0 .89125

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

POWER DENSITY= 1.000000

NUMBER OF HEAT SOURCES= 1 WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

57
TML I/O FILE LISTING - tmlio.220

210 10

20

100000 . 100000

100000 . 112202

100000 _i_ \y \j w \j kj . 125893

. 100000 . 141254

. 100000 . 158489

. 100000 . 177828

. 100000 .199526

. 100000 .223872

. 100000 .251189

. 100000 .281838

. 100000 .316228

. 100000 . 354813

.100000 .398107

100000 446683

.100000 .501187

. 100000 .562341

.100000 .630957

. 100000 .707945

. 100000 .794328

100000 . 891250

.0 1.

0.0

.0

0 2.0 4.0 2.0

14.0 2.0

The above is the input used in the input.dat file for the example of a tml line average calculation for a 20 layer structure.

This input is annotated below.

INPUT

210 10

20

DESCRIPTION OF INPUT

itype (=2 for the line average calculation)

x and y dimensions of rectangular structure

number of layers in the structure

.100000 .100000 thickness and thermal conductivity of layer 20

.100000 .112202 thickness and thermal conductivity of layer 19

.100000 .125893 thickness and thermal conductivity of layer 18

.100000 .141254 thickness and thermal conductivity of layer 17

.100000 .158489 thickness and thermal conductivity of layer 16

.100000 .177828 thickness and thermal conductivity of layer 15

.100000 .199526 thickness and thermal conductivity of layer 14

.100000 .223872 thickness and thermal conductivity of layer 13

.100000 .251189 thickness and thermal conductivity of layer 12

.100000 .281838 thickness and thermal conductivity of layer 11

.100000 .316228 thickness and thermal conductivity of layer 10

.100000 .354813 thickness and thermal conductivity of layer 9 .100000 .398107 thickness and thermal conductivity of layer 8 .100000 .446683 thickness and thermal conductivity of layer 7 .100000 .501187 thickness and thermal conduct ivit

y

of layer 6 .100000 .562341 thickness and thermal conductivity of layer 5 .100000 . 630957 thickness and thermal conductivity of layer 4 .100000 .707945 thickness and thermal conductivity of layer 3

.100000 .794328 thickness and thermal conductivity of layer 2

.100000 .891250 thickness and thermal conductivity of layer 1

58
TML I/O FILE LISTING - tmlio.220

50 50

111 0.0

115. 0

1 1.0

1 4.0

14.0 2.0

0.0

2.0

1.

4.0 2.0

upper summation limits for n and m summations

iedgex (=1, then read in three values on next line)

number of values, first point, increment

iedgey (=1, then read in three values on next, line)

number of values, first point, increment

number of sources and power density

weight, x, length along x, y, length along y

nline as itype=2 (=1 for one line, etc)

x, length along x for line element

Below is the output contained in the output.dat file which is calculated

and written by tml using the above input for a 20 layer structure.

STEADY-STATE THERMAL MULTILAYER CALCULATION

USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS

LINE AVERAGE EVALUATION OF SURFACE TEMPERATURE

LINE # XLINE LXLINE Y AVE TEMP

1 4.00000 2.00000 5.00000 5.19144

LX= 10.00 LY= 10.00

NUMBER OF LAYERS IN STRUCTURE=2 0 THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L2 0= 0..10000 K20= 0 .10000

L19= 0..10000 K19= 0 .11220

L18= 0..10000 K18= 0 . 12589

L17 = 0..10000 K17 = 0 .14125

LI 6= 0..10000 K16= 0 . 15849

L15= 0..10000 K15= 0 .17783

1,14= 0..10000 K14= 0 .19953

L13= 0..10000 K13= 0 .22387

L12= 0..10000 K12= 0 .25119

Lll= 0..10000 Kll= 0 .28184

L10= 0..10000 K10= 0 .31623

L 9= 0..10000 K 9= 0 .35481

L 8= 0..10000 K 8= 0 .39811

L 7 = 0..10000 K 7= 0 .44668

L 6= 0..10000 K 6= 0 .50119

L 5= 0.,10000 K 5= 0 .56234

L 4= 0..10000 K 4= 0 .63096

L 3= 0..10000 K 3= 0 .70794

L 2= 0..10000 K 2= 0 .79433

L 1= 0..10000 K 1 = 0 .89125

UPPER SUMMATION LIMITS NUP= 50 MUP= 50

POWER DENSITY= 1.000000

NUMBER OF HEAT SOURCES= 1 WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

NUMBER OF LINES= 1

1 4.00000 2.00000

59
TML I/O FILE LISTING - tmlio.320

310 10

20

50 50

111 0

1

t rt r\ r\ r\ rt lUUUUU .lUUUUU

1 r\ r\ r\ r\ r\ 100000 i i o o n o

. llzzOz

100000 .125893

lUUUUU . 141zo4

i n r\ rt r\ f\ 1U0U0U t c o yi o a lUUUUU . 1 / / OZ O

ilUUUUU n o n n n

i n n n n n x u u u u u . Z Z O O / Z -LUUUUU i nnnnn . Z 3 X X O _7

i-LUUUUU nnnnn . ^ U 1 U JO

-LUUUUU J 1 U/iii o

iX nnnnn u u u u u • J J4010

100000 .398107

ix nnnnn u u u u u . If! ODOJ

100000 .501187

100000 .562341

100000 .630957

lUUUUU . / u / y4o

lUUUUU 9 Q /I "3 9 O

i f\ r\ r* rs lUUUUU . ayizou

0 1.

0.0

0

2.0 4 0 2.0

5.

1.

4.0

0 2.0 4.0 2.0

The above is the input used in the input.dat file for the example of a tml area average calculation for a 20 layer structure.

This input is annotated below.

INPUT

3

10 10

20

. 100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

.100000

100000

,112202

,125893

,141254

,158489

,177828

,199526

,223872

,251189

,281838

.316228

,354813

,398107

,446683

.501187

,562341

.630957

,707945

.794328

,891250

DESCRIPTION OF INPUT

itype (=3 for the line average calculation)

x and y dimensions of rectangular structure

number of layers in the structure

thickness and thermal conductivity

thickness and thermal conductivity

thickness and thermal conductivity

thickness and thermal conductivity

thickness and thermal conductivity o thickness and thermal conductivity

thickness and thermal conductivity

thickness and thermal conductivity

thickness and thermal conductivity

thickness and thermal conductivity

of layer 20

of layer 19

of layer 18

of layer 17

of layer 16

of layer 15

of layer 14

of layer 13

of layer 12

of layer 11

of layer 10

of layer 9 of layer 8 of layer 7 of layer 6 of layer 5 of layer 4 of layer 3

of layer 2

of layer 1

60
TML I/O FILE LISTING - tmlio.320

50 50

111 0.0 1.

1

1 5. 0.0

1 1.0

1 4.0 2.0 4.0 2.0

14.0 2.0 4.0 2.0

upper summation limits for n and m summations

iedgex (=1, then read in three values on next line)

number of values, first point, increment

iedgey (=1, then read in three values on next line)

number of values, first point, increment

number of sources and power density

weight, x, length along x, y, length along y narea as itype=3 (=1 for one area, etc)

x, length along x, y, length along y for area element

Below is the output contained in the output.dat file which is calculated

and written by tml using the above input for a 20 layer structure.

STEADY-STATE THERMAL MULTILAYER CALCULATION

USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS

AREA AVERAGE EVALUATION OF SURFACE TEMPERATURE

AREA # XAREA LXAREA YAREA LYAREA AVE TEMP

1 4.00000 2.00000 4.00000 2.00000 4.64315

LX= 10.00 LY= 10.00

NUMBER OF LAYERS IN STRUCTURE=20

THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS

L2 0= 0 10000 K20= 0 10000

LI 9= 0 10000 K19= 0 11220

L18= 0 10000 K18= 0 12589

L17= 0 10000 K17= 0 14125

LI 6= 0 10000 K16= 0 15849

L15= 0 10000 K15= 0 17783

L14= 0 10000 K14= 0 19953

L13= 0 10000 K13= 0 22387

L12= 0 10000 K12= 0 25119

Lll = 0 10000 Kll= 0 28184

L10= 0 10000 K10= 0 31623

L 9= 0 10000 K 9= 0 35481

L 8= 0 10000 K 8= 0 39811

L 7 = 0 10000 K 7= 0 44668

L 6= 0 10000 K 6= 0 .50119

L 5= 0 10000 K 5= 0 56234

L 4= 0 10000 K 4= 0 63096

L 3= 0 10000 K 3= 0 70794

L 2= 0 10000 K 2= 0 79433

L 1= 0 10000 K 1= 0 89125

UPPER SUMMATION LIMITS NUP=

POWER DENSITY= 1.000000

NUMBER OF HEAT SOURCES= 1

50 MUP= 50

WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES

HEAT SOURCE WTSOUR XSOUR YSOUR LXSOUR LYSOUR

1 1.00000 4.00000 4.00000 2.00000 2.00000

NUMBER OF AR£AS= 1

1 4.00000 2.00000 4.00000 2.00000

61
APPENDIX A - TXYZ30 LISTING

q* *****************************************************************************

C TXYZ30 - TXYZ VERSION 3.0 - VERSION DATE 06/27/95

C

C THIS IS VERSION 3 . 0 OF THE TXYZ THERMAL ANALYSIS CODE. THIS VERSION

C RETAINS THE CHANGES MADE IN VERSION 2.0 (ASSIGNMENT OF ARBITRARY WEIGHTS,

C BOTH + AND -, TO THE HEAT SOURCES ON THE SURFACE AND THE FIX OF A PROBLEM

C WITH THICK MIDDLE LAYERS) AND INCORPORATES THE FOLLOWING ADDITIONS:

C 1) ANALYTIC AVERAGING OF THE TEMPERATURE OVER ARBITRARY LINE SEGMENT (S),

C 2) ANALYTIC AVERAGING OF THE TEMPERATURE OVER ARBITRARY AREA ( S )

C 3) RENUMBERING OF LAYERS TO BE IN KEEPING WITH THE NOTATION TO BE USED

C IN THE MULTILAYER THERMAL PROBLEM (SEE TML CODE)

C THE LINE AND AREA AVERAGING ARE PERFORMED ANALYTICALLY FOR UNIFORM

C WEIGHTING OF THE AVERAGE. THE AREA AVERAGE SHOULD PROVIDE A CLOSER

C LINK WITH EXPERIMENT WHERE THE MEASUREMENT DOES AN AVERAGING OVER

C SOME AREA ELEMENT (ELECTRICALLY OR OPTICALLY)

C

C THIS PROGRAM CALCULATES THE STEADY-STATE TEMPERATURE DISTRIBUTION

C FOR A RECTANGULAR THREE LAYER STRUCTURE WITH AN ARBITRARY NUMBER

C OF RECTANGULAR HEAT SOURCES/SINKS ON THE TOP SURFACE.

C THE CALCULATION FOLLOWS FROM THE INPUT OF THE THICKNESSES AND THERMAL

C CONDUCTIVITIES OF THE THREE LAYERS.

C THE TEMPERATURE (RELATIVE TO THAT OF THE BOTTOM HEAT SINK) MAY BE

C CALCULATED AS A POINT FUNCTION, A LINE AVERAGE, OR AN AREA AVERAGE

C ANYWHERE ON OR INSIDE THE THREE-LAYER STRUCTURE.

C IT IS IMPORTANT TO EMPHASIZE THAT THE CALCULATION IS GENERAL FOR THE

C THREE LAYER STRUCTURE AND THE APPLICATION TO SEMICONDUCTOR STRUCTURES

C IS A SPECIAL CASE.

C THE STARTING EQUATIONS USED ARE GIVEN IN EQUATIONS (13)-(23), WITH S=0

C (ZERO FREQUENCY, STEADY-STATE CONDITION), IN THE PAPER BY KOKKAS (BELOW).

C THE LINE AND AREA AVERAGES FOLLOW FROM ANALYTICAL EVALUATION OF THESE

C EQUATIONS

C

C REFERENCES: THE ORIGINAL MATHEMATICAL ANALYSIS OF THE THREE-LAYER

C STRUCTURE WAS PERFORMED IN THE PAPER "THERMAL ANALYSIS

C OF MULTIPLE-LAYERED STRUCTURES" BY ACHILLES G. KOKKAS,

C IEEE TRANS. ELEC . DEV. VOL. ED-21, NO. 11, 674-681 (1974).

C THIS PAPER WAS DRAWN FROM HIS PHD THESIS: A. G. KOKKAS,

C "ANALYSIS AND DESIGN OF ELECTROTHERMAL INTEGRATED CIRCUITS,

C PH.D. THESIS, MIT, 1972.

C

C THE ORIGINAL FORTRAN IMPLEMENTATION OF THE STEADY STATE

C KOKKAS EQUATIONS IS CONTAINED IN THE TXYZ CODE AND

C WAS PRESENTED IN THE REPORT "SEMICONDUCTOR MEASUREMENT

C TECHNOLOGY: TXYZ: A PROGRAM FOR SEMICONDUCTOR IC THERMAL

C ANALYSIS" BY JOHN ALBERS, NBS SPECIAL PUBLICATION 400-7 6 C (APRIL 1984)

C

C VERSION 2.0 OF THE TXYZ CODE WAS PRESENTED IN THE REPORT

C "SEMICONDUCTOR MEASUREMENT TECHNOLOGY: VERSION 2 . 0 OF THE

C TXYZ THERMAL ANALYSIS PROGRAM: TXYZ2 0" BY JOHN ALBERS,

C NIST SPECIAL PUBLICATION 400-89 (JUNE 1992)

C

C VERSION 3.0 OF THE TXYZ CODE AND THE ACCOMPANYING THERMAL

C MULTILAYER CODES ARE DISCUSSED IN THIS REPORT

C "SEMICONDUCTOR MEASUREMENT TECHNOLOGY: HOTPAC : PROGRAMS

C FOR THERMAL ANALYSIS INCLUDING VERSION 3 . 0 OF THE TXYZ

C PROGRAM, TXYZ30, AND THE THERMAL MULTILAYER PROGRAM, TML"

C BY JOHN ALBERS, NIST SPECIAL PUBLICATION 400-96.

Q* ******************************************************** ^ ^ ^ * ^ * ^ ^ * ^* **** *^ *

62
APPENDIX A - TXYZ30 LISTING

C IN THE PRESENT FORM, THE PROGRAM ALLOWS UP TO 500 TERMS TO BE INCLUDED

C IN BOTH THE N SUM (ALONG X) AND THE M SUM (ALONG Y) . TO GO BEYOND

C THIS NUMBER, SUBSTITUTE THE FOLLOWING TWO LINES WITH THE APPROPRIATE

C VALUES OF NX AND MY FOR THE FIRST TWO DIMENSION STATEMENTS - MAKE SURE TO

C REMOVE THE COMMENTS FROM THE NEW LINES AND COMMENT OUT THE REPLACED LINES,

C DIMENSION X( 100) , Y(100), Z(100), COSYT(MY)

C DIMENSION ARUZER (NX, MY) , ARFUNZ (NX, MY)

C ALSO REPLACE THE TESTING LINE ABOVE THE LINE LABELLED 110 WITH

C IF (NUP . GT . NX . OR . MUP . GT . MY ) GO TO 3999

Q* ***************************************************************************** DIMENSION X (100) , Y(100), Z(100), COSYT(500)

DIMENSION ARUZER(500, 500) , ARFUNZ (500, 500)

DIMENSION WTSOUR (50) , XSOUR(50), YSOUR(50)

DIMENSION XLINE(30), XAREA(30), YAREA (30)

REAL LXLINE(30), LXAREA (30) , LYAREA (30)

REAL LXSOUR(50), LYSOUR(50), K3, K2, Kl, LX, LY, L3, L2, LI

COMMON K3, K2, Kl, LX, LY, L3, L2, LI

COMMON NSOUR, WTSOUR, XSOUR, YSOUR, LXSOUR, LYSOUR

C INPUT DATA IS READ FROM 10 AND OUTPUT GOES TO 12.

C THESE ARE WRITTEN IN LOWER CASE. MANY OPERATING SYSTEMS TAKE

C UPPER AND LOWER CASE AS EQUIVALENT. THE UNIX OPERATING SYSTEM

C VIEWS THE UPPER AND LOWER CASE NAMED FILES AS DIFFERENT. UNIX

C DEFAULTS TO THE LOWER CASE WHICH IS USED HERE.

open ( unit=10, file=' input . dat' , status = 'unknown')

open ( unit=12, file=' output .dat' , status = 'unknown')

PI=3 . 14159265

STEADY-STATE THERMAL CALCULATION KOKKAS ANALYSIS'

)

THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS'

)

2 FORMAT (IX

3 FORMAT (IX

4 FORMAT (IX

5 FORMAT (IX

6 FORMAT (IX

1 ' MUP='

7 FORMAT (IX

8 FORMAT (IX

9 FORMAT (IX

1 ' LXSOUR'

10 FORMAT (4X

11 FORMAT (IX

27 FORMAT (IX

22 FORMAT (IX

31 FORMAT (IX

32 FORMAT ( IX

33 FORMAT (IX

34 FORMAT (IX

35 FORMAT (IX

40 FORMAT (IX

41 FORMAT (IX

42 FORMAT (IX

43 FORMAT (6X

44 FORMAT (IX

45 FORMAT (IX

l'AVE TEMP

46 FORMAT (IX

47 FORMAT (IX

1 ' LYAREA'

48 FORMAT (IX

11X,F10.5)

5)

5)

'K3= ',F10.5,' K2= ',F10.5,' Kl= ',F10

'L3= ',F10.5,' L2= ' ,F10.5, ' Ll= ',F10

'UPPER SUMMATION LIMITS ' , 2X, ' NUP=' , 15,

15)

' NUMBER OF HEAT SOURCES=' , 15)

'WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES')

' HEAT SOURCE' , 4X, ' WTSOUR' , 8X, ' XSOUR' , 8X, ' YSOUR' , 7X,

7X, ' LYSOUR' ) 13, 5X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5)

'POWER DENSITY=' , Fll . 6)

'LX= ',F7.2,3X, ' LY= ',F7.2)

F12 . 4, 2X,F12 . 4, 2X, F12 . 4, 2X,F12 . 4)

617)

'NUMBER OF LINES=',I3)

13, 3X,F10.5, 3X,F10.5)

'NUMBER OF AREAS=' , 13)

13, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5)

'POINT FUNCTION EVALUATION OF T(X,Y,Z)')

'LINE AVERAGE EVALUATION OF TEMPERATURE')

'AREA AVERAGE EVALUATION OF TEMPERATURE')

' X' , 12X, ' Y' , 12X, ' Z' , 9X, ' T (X, Y, Z) '

) F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5)

' LINE #' , 5X, ' XLINE' , 7X, ' LXLINE' , 8X, ' Y' , 12X, ' Z' , 9X,

)13, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3x, flO

' AREA #' , 5X, ' XAREA' , 7X, ' LXAREA' , 8X, ' YAREA' , 7X,

6X, ' Z' , 7x, 'AVE TEMP'

)

13, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X, FlO. 5, IX, FlO

5)

5,

63
APPENDIX A - TXYZ30 LISTING

51 FORMAT (IX, 14, 3X, 14

52 FORMAT (IX, F10 . 5, 3X, F10 . 5)

53 FORMAT (IX, 11)

54 FORMAT (IX, 14, 3X,F10.5, 3X,F10.5)

55 FORMAT ( IX, Fl 0.5)

56 FORMAT (IX, 12, 3X,F10.5)

57 FORMAT (IX, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 5)

88 FORMAT (IX, 'YOUR UPPER LIMIT OF SUMMATION IS TOO LARGE. TRY AGAIN')

q* *****************************************************************************

INPUT SECTION

=1 FOR POINT FUNCTION, T(X,Y,Z), EVALUATION

=2 FOR LINE AVERAGE, <T(Y, Z)>, EVALUATION

=3 FOR AREA AVERAGE, <T(Z)>, EVALUATION

c

c

c THE FOLLOWING

c

c I TYPE

c

c

c

c LX

c LY

c L3

c K3

c L2

c K2

c LI

c Kl

c NUP

c MUP

c

c IEDGEX

cc

c

c

c

c

cc

c

c

c IEDGEY

cc

c

c

c

cc

cc

c

THICKNESS OF TOP LAYER

THERMAL CONDUCTIVITY OF TOP LAYER

THICKNESS OF BOTTOM LAYER

THERMAL CONDUCTIVITY OF BOTTOM LAYER

UPPER LIMIT OF N SUM, X DIRECTION

UPPER LIMIT OF M SUM, Y DIRECTION

INDEX FOR GENERATING THE VALUES OF X TO BE USED

=1 TO READ DATA FOR FIXED INCREMENT X VALUES

=2 TO READ IN ARRAY OF X VALUES OF NONFIXED INCREMENT)

IF IEDGEX=1 THEN READ THE THREE VARIABLES (ON SAME LINE)

ILX THE NUMBER OF X VALUES TO BE USED

XI THE VALUE OF THE FIRST POINT IN X

STEPX (THE INCREMENT IN X)

IF IEDGEX=2 THEN READ THE VARIABLE AND ARRAY (ONE PER LINE)

ILX THE NUMBER OF X VALUES TO BE USED

X(I) THE ARRAY OF X VALUES (1=1, ILX)

=1 TO READ DATA FOR FIXED INCREMENT Y VALUES

=2 TO READ IN ARRAY OF X VALUES OF NONFIXED INCREMENT)

IF IEDGEY=1 THEN READ THE THREE VARIABLES (ON SAME LINE)

ILY THE NUMBER OF Y VALUES TO BE USED

Yl THE VALUE OF THE FIRST POINT IN Y

STEPY (THE INCREMENT IN Y)

IF IEDGEY=2 THEN READ THE VARIABLE AND ARRAY (ONE PER LINE)

ILY THE NUMBER OF Y VALUES TO BE USED

Y(I) THE ARRAY OF Y VALUES (1=1, ILY)

64
APPENDIX A - TXYZ30 LISTING

C IEDGEZ INDEX FOR GENERATING THE VALUES OF Z TO BE USED

C =1 TO READ DATA FOR FIXED INCREMENT Z VALUES

C =2 TO READ IN ARRAY OF Z VALUES OF NONFIXED INCREMENT)

C IF IEDGEZ=1 THEN READ THE THREE VARIABLES (ON SAME LINE)

C ILZ THE NUMBER OF Z VALUES TO BE USED

C Zl THE VALUE OF THE FIRST POINT IN Z C STEPZ (THE INCREMENT IN Z)

C IF IEDGEZ=2 THEN READ THE VARIABLE AND ARRAY (ONE PER LINE)

C ILZ THE NUMBER OF Z VALUES TO BE USED

C Z(I) THE ARRAY OF Z VALUES (1=1, ILZ)

C

C NOTE: ENTER THE Z-RELATED VARIABLES (Zl, STEPZ OR THE Z(I) ARRAY)

C AS ZERO OR POSITIVE QUANTITIES. THE PROGRAM CONVERTS THE FINAL

C Z(I) ARRAY TO ZERO OR NEGATIVE QUANTITIES AS THE CALCULATION

C TAKES THE Z VARIABLE TO BE ZERO OR NEGATIVE.

C

C IMPORTANT: THE POINT FUNCTION EVALUATION OF THE TEMPERATURE IS

C THE MOST ELEMENTAL VERSION. IN ORDER TO SIMPLIFY THE

C CONSTRUCTION OF THE DATA FILES FOR LINE AND AREA

C AVERAGES, THE PROGRAM EXPECTS TO SEE THE ABOVE IEDGEX,

C IEDGEY AND IEDGEZ DATA. THIS IS READ EVEN IF IT IS NOT

C USED FOR THE AVERAGE VERSIONS. HOWEVER, THE LINE AND

C AREA INFORMATION MAY THEN BE SIMPLY APPENDED TO THE END

C OF THE DATA FILE IN ORDER TO RUN THESE VERSIONS. SEE THE

C SAMPLE I/O FILES FOR AN EXAMPLE OF THIS.

C

C NSOUR NUMBER OF HEAT SOURCES (UP TO 50)

C P0 POWER DENSITY

C

C THE NEXT NSOUR LINES READ THE FOLLOWING INFORMATION FOR THE

C HEAT SOURCES (WITH ALL THE INFORMATION FOR EACH OF THE ELEMENTS

C ON A SINGLE LINE)

C

C WTSOUR( I) -WEIGHTING FACTOR OF I-TH SOURCE

C (POSITIVE FOR SOURCE, NEGATIVE FOR SINK)

C XSOUR(I) —X COORDINATE OF ORIGIN OF I-TH SOURCE

C LXSOUR (I) -LENGTH ALONG X DIRECTION OF I-TH SOURCE

C YSOUR(I) —Y COORDINATE OF ORIGIN OF I-TH SOURCE

C LYSOUR (I) -LENGTH ALONG Y DIRECTION OF I-TH SOURCE

C

C IF ITYPE=1, THE POINT FUNCTION CALCULATION CONTINUES WITH THE ABOVE

C SET OF (X,Y, Z) VALUES

C IF ITYPE=2, THE LINE AVERAGE CALCULATION READS THE FOLLOWING:

C NLINE THE NUMBER OF LINE SEGMENTS TO DO THE AVERAGE

C THE NEXT NLINE LINES THEN READ

C XLINE(J), LXLINE(J) —THE LOCATION AND LENGTH OF THE

C J-TH LINE ELEMENT

C IF ITYPE=3, THE AREA AVERAGE CALCULATION READS THE FOLLOWING:

C NAREA THE NUMBER OF AREA SEGMENTS TO DO THE AVERAGE

C THE NEXT NAREA LINES THEN READ

C XAREA (J) , LXAREA (J) , YAREA (J) , LYAREA (J) -THE LOCATION

C AND LENGTHS OF THE J-TH AREA ELMENT *********************************************************** *

65
APPENDIX A - TXYZ30 LISTING

C READ I TYPE (ANALYSIS TYPE, 1 FOR POINT, 2 FOR LINES, 3 FOR AREAS)

READ (10, *) I TYPE

C READ LX AND LY (THE X AND Y DIMENSIONS OF THE RECTANGULAR STRUCTURE)

READ (10,*) LX, LY

C READ L3 AND K3 (THICKNESS AND THERMAL CONDUCTIVITY OF THE TOP LAYER)

READ (10,*) L3, K3

C READ L2 AND K2 (THICKNESS AND THERMAL CONDUCTIVITY OF THE MIDDLE LAYER)

READ (10,*) L2, K2

C READ LI AND Kl (THICKNESS AND THERMAL CONDUCTIVITY OF THE BOTTOM LAYER)

READ (10, * ) LI, Kl

C READ NUP AND MUP (UPPER LIMIT OF THE SUMMATION OVER THE INDEX N (X-DIR)

C UPPER LIMIT OF THE SUMMATION OVER THE INDEX M (Y-DIR)

READ (10, *) NUP, MUP

C NUP AND MUP MUST BE LESS THAN OR EQUAL TO THE DIMENSIONALITY OF

C ARUZER AND ARFUNZ

C THE NEXT LINE ALLOWS VALUES OF NUP AND MUP UP TO THE PRESENT DIMENSION

C OF ARUZER AND ARFUNZ. TO GO BEYOND THIS VALUE, THE USER SHOULD EDIT

C THE DIMENSION STATEMENTS ACCORDINGLY AND THEN COMMENT OUT OR MODIFY

C THE NEXT LINE OF CODE

IF (NUP .GT. 500 .OR. MUP .GT. 500) GO TO 3999

READ (10, *) IEDGEX

GOTO (110, 115) IEDGEX

110 READ (10, *) ILX, XI, STEPX

DO 111 1=1, ILX

X(I)=X1+(I-1) *STEPX

111 CONTINUE

GOTO 119

115 READ (10,*) ILX

DO 116 1 = 1, ILX

READ (10, *)X (I)

116 CONTINUE

119 READ (10, *) IEDGEY

GOTO (120,125) IEDGEY

120 READ (10, *) ILY, Yl, STEPY

DO 121 1=1, ILY

Y (I) =Y1+ (1-1) *STEPY

121 CONTINUE

GOTO 12 9 125 READ (10,*) ILY

DO 126 1=1, ILY

READ (10, *) Y (I)

12 6 CONTINUE

129 READ (10, *) IEDGEZ

GOTO (130,135) IEDGEZ

130 READ (10, *) ILZ, Zl, STEPZ

Zl=-1 . 0*Z1

STEPZ=-1 . 0*STEPZ

DO 131 1=1, ILZ

Z (I) =Z1+ (1-1) *STEPZ

131 CONTINUE

GOTO 139

135 READ (10,*) ILZ

DO 136 1=1, ILZ

READ (10, *) Z (I) Z (I) = -1 . 0*Z (I)

13 6 CONTINUE

66
APPENDIX A - TXYZ30 LISTING

C READ THE NUMBER OF HEAT SOURCES AND THE POWER DENSITY

C NOTE-POWER DENSITY IS MULTIPLICATIVE FACTOR USUALLY SET EQUAL TO UNITY

139 READ (10, *)NSOUR,P0

C P0 IS THE POWER DENSITY, ASSUMED UNIFORM FOR ALL HEATERS

C NSOUR IS THE TOTAL NUMBER OF HEATING ELEMENTS ON THE SURFACE OF THE

C THE TOP LAYER (UP TO 50)

C THE NEXT LOOP READS IN THE COORDINATES OF THE ORIGIN OF THE

C HEATING ELEMENTS ALONG WITH THEIR LENGTHS AND WIDTHS

C THE WEIGHTING FACTOR IS ALSO ENTERED (THIS IS REAL, NONINTEGER)

DO 140 1=1, NSOUR

READ (10, * ) WTSOUR ( I ), XSOUR ( I ) , LXSOUR ( I ) , YSOUR ( I ) , LYSOUR(I)

C WTSOUR (I) IS THE WEIGHTING FACTOR FOR THE I-TH HEATER ELEMENT

C XSOUR (I) IS THE X COORDINATE OF THE ORIGIN OF I-TH HEATER ELEMENT

C LXSOUR (I) IS THE LENGTH OF THE I-TH HEATER ALONG THE X DIRECTION

C YSOUR (I) IS THE Y COORDINATE OF THE ORIGIN OF I-TH HEATER ELEMENT

C LYSOUR(I) IS THE LENGTH OF THE I-TH HEATER ALONG THE Y DIRECTION

140 CONTINUE

C THIS IF. . .THEN. . .ELSE IF CONSTRUCTION IS USED TO READ IN THE DATA FOR

C THE LINES OR AREAS TO BE CONSIDERED IN THE CALCULATION

C IF THE ANALYSIS IS FOR A POINT FUNCTION, THEN GO OUT OF THE IF. THEN. ELSE

IF (ITYPE.EQ.l) THEN

GOTO 17 0 C IF THE ANALYSIS IS FOR A LINE AVERAGE, THEN READ THE NUMBER OF LINE

C SEGMENTS AND THEIR LOCATION

ELSE IF (ITYPE.EQ.2) THEN

C READ THE NUMBER OF LINE SEGMENTS

READ (10, * ) NLINE

DO 150 J=l, NLINE

C READ THE ORIGIN AND LENGTH OF EACH LINE SEGMENT

C NOTE THAT THE AVERAGE IS ALONG THE X-DIRECTION FOR GIVEN Y VALUES

C TO DO THE AVERAGE ALONG THE Y-DIRECTION FOR GIVEN X, SIMPLY ROTATE

C THE STRUCTURE BY 90 DEGREES AND USE THE CORRESPONDING NEW X' S AND Y'

S

C XLINE(J) IS THE X COORDINATE OF THE ORIGIN OF J-TH LINE ELEMENT

C LXLINE(J) IS THE LENGTH OF THE J-TH LINE ALONG THE X DIRECTION

READ (10, * ) XLINE (J) , LXLINE (J)

150 CONTINUE

C IF THE ANALYSIS IS FOR A AREA AVERAGE, THEN READ THE NUMBER OF AREAS

C AND THEIR LOCATIONS

ELSE IF (ITYPE.EQ.3) THEN

C READ THE NUMBER OF AREAS

READ (10, * ) NAREA

DO 160 J=l, NAREA

C READ THE FOLLOWING

C XAREA(J) IS THE X COORDINATE OF THE ORIGIN OF J-TH AREA ELEMENT

C LXAREA(J) IS THE LENGTH OF THE J-TH AREA ALONG THE X DIRECTION

C YAREA(J) IS THE Y COORDINATE OF THE ORIGIN OF J-TH AREA ELEMENT

C LYAREA (J) IS THE LENGTH OF THE J-TH AREA ALONG THE Y DIRECTION

READ (10, * ) XAREA ( J) , LXAREA (J) , YAREA (J) , LYAREA (J)

160 CONTINUE

END IF

17 0 CONTINUE

P04LK = 4.0 * P0 / ( LX* LY* K3)

PILX = PI / LX

PILY = PI / LY

67
APPENDIX A - TXYZ30 LISTING

Q* *****************************************************************************

c

C END OF DATA INPUT SECTION

C BEGIN CALCULATION OF TEMPERATURE

C THE SUBROUTINES USED IN THE CALCULATION ARE:

C 1) UZERO(N,M) - CALCULATES THE FOURIER COSINE TRANSFORM OF THE

C FUNCTION, U(X,Y), THE POWER DENSITY FUNCTION FOR ALL OF THE

C HEAT SOURCES.

C 2) FUNZ(N,M, Z) - CALCULATES THE Z-DEPENDENT PORTION OF THE SUM

C REMEMBERING THAT THIS IS A FUNCTION OF THE SUMMATION

C INDICES (N, M)

Q* *****************************************************************************

C CALCULATE THE FOURIER COMPONENTS OF THE HEAT SOURCES, U (N, M) ^* *****************************************************************************

DO 300 MM=1,MUP

M = MM - 1 DO 250 NN=1,NUP

N = NN - 1 ARUZER (NN, MM) =UZERO (N, M)

2 50 CONTINUE

300 CONTINUE

q* *****************************************************************************

C END OF U(N,M) CALCULATION AND BEGINNING OF MAJOR LOOP FOR Z Q* *****************************************************************************

DO 6000 IZ=1, ILZ

Q* *****************************************************************************

C CALCULATE THE Z DEPENDENT PORTION, I.E., FUNZ(N,M, Z)

(2* **************************************************************************** * DO 400 MM=1,MUP

M = MM - 1 DO 350 NN=1,NUP

N = NN - 1 ARFUNZ (NN, MM) =FUNZ (N, M, Z (IZ) ) * ARUZER (NN, MM)

350 CONTINUE

400 CONTINUE

C THE FOLLOWING IF. THEN. ELSE CONSTRUCTION OPERATES ACCORDING TO THE

C TYPE OF ANALYSIS TO BE USED.

C FOR THE POINT FUNCTION ANALYSIS, THE 3000 LOOP IS USED

C FOR THE LINE AVERAGE ANALYSIS, THE 4000 LOOP IS USED

C FOR THE AREA AVERAGE ANALYSIS, THE 5000 LOOP IS USED

C THIS PORTION DOES THE POINT FUNCTION CALCULATION

IF (ITYPE.EQ.l) THEN

C BEGINNING OF POINT FUNCTION ANALYSIS

WRITE (12, 2)

WRITE (12, 40)

WRITE (12, 43)

DO 3000 IY=1, ILY

DO 3100 MM=1,MUP

M = MM - 1 COSYT (MM) =COS (FLOAT (M) *Y (IY) *PILY)

310 0 CONTINUE

68
APPENDIX A - TXYZ30 LISTING

DO 3000 IX=1, ILX

SUM=0 . 0 DO 3300 MM=1,MUP

M = MM - 1 DO 3200 NN=1,NUP

N = NN — 1 NDN=0

NDM=0

IF (N.EQ.O) NDN=1

IF (M.EQ.O) NDM=1

TOP = ARFUNZ (NN, MM) * COS (FLOAT (N) *X (IX) *PILX) * COSYT (MM)

BOTTOM= (NDN+1) * (NDM+1)

TSUM=TOP /BOTTOM

SUM=SUM+TSUM

32 0 0 CONTINUE

3300 CONTINUE

TEMP = P04LK * SUM

WRITE (12, 44)X(IX) , Y (IY) , Z (IZ) , TEMP

3000 CONTINUE

C END OF POINT FUNCTION CALCULATION

ELSE IF (ITYPE.EQ.2) THEN

C THIS PORTION DOES THE LINE AVERAGE CALCULATION

WRITE (12, 2)

WRITE (12, 41)

WRITE (12, 45)

DO 4000 IY=1, ILY

DO 4010 MM=1,MUP

M = MM - 1 COSYT (MM) =COS (FLOAT (M) *Y ( IY) *PILY)

4 010 CONTINUE

DO 4000 J=1,NLINE

SUM=0 .

DO 4200 MM=1,MUP

M = MM - 1 DO 4100 NN=1,NUP

N = NN - 1 NDN=0

NDM=0

IF (N.EQ.O) NDN=1

IF (M.EQ.O) NDM=1

IF(N.EQ.O) GO TO 4160

TERMX=SIN (FLOAT (N) *PI* (XLINE (J) +LXLINE (J) ) /LX)

1 -SIN (FLOAT (N) *PI*XLINE (J) /LX)

TERMX=TERMX*LX/ (FLOAT (N) *PI)

GO TO 4165

4160 TERMX=LXLINE (J)

4165 CONTINUE

TOP = ARFUNZ (NN, MM) * COSYT (MM) * TERMX

EOTTOM= (NDN+1) * (NDM+1)

TSUM=TOP /BOTTOM

SUM=SUM+TSUM

4100 CONTINUE

TEMP = P04LK * SUM / LXLINE(J)

42 00 CONTINUE

WRITE (12, 46) J, XLINE (J) , LXLINE (J) , Y (IY) , Z (IZ) , TEMP

4000 CONTINUE

69
APPENDIX A - TXYZ30 LISTING

C END OF LINE AVERAGE PORTION OF THE CODE

ELSE IF (ITYPE.EQ.3) THEN

C THIS PART PERFORMS THE AREA AVERAGE CALCULATION

WRITE (12, 2)

WRITE (12, 42)

WRITE (12, 47)

DO 5000 J=l , NAREA

SUM=0 .

DO 5100 MM=1,MUP

M = MM — 1 DO 52 00 NN=1,NUP

N = NN - 1 NDN=0

NDM=0

IF (N.EQ.0) NDN=1

IF (M.EQ.0) NDM=1

AREA=0 .

IF(N.EQ.O) GO TO 5160

TERMX = SIN (FLOAT (N) *PI* (XAREA (J) + LXAREA ( J) ) / LX) - SIN (FLOAT (N) *PI*XAREA (J) / LX)

TERMX=TERMX* LX/ (FLOAT (N) *PI

GO TO 5165

TERMX= LXAREA (J)

IF(M.EQ.O) GO TO 5164

TERMY = SIN (FLOAT (M) *PI* (YAREA (J) + LYAREA ( J ) ) / LY) - SIN (FLOAT (M) *P I* YAREA (J) / LY)

TERMY=TERMY* LY/ (FLOAT (M) *PI)

GO TO 5166

TERMY= LYAREA (J)

TERMI=TERMX* TERMY

AREA=TERMI

TOP = ARFUNZ (NN, MM) * AREA

BOTTOM= (NDN+1) * (NDM+1)

TSUM=TOP /BOTTOM

SUM=SUM+TSUM

CONTINUE

TEMP = P04LK * SUM / (LXAREA ( J) *LYAREA ( J)

CONTINUE

WRITE (12, 48) J, XAREA (J) , LXAREA (J) , YAREA (J) , LYAREA ( J) , Z (IZ) , TEMP

5000 CONTINUE

END IF

6000 CONTINUE

WRITE (12, 27) LX, LY

WRITE (12, 3)

WRITE(12,5) L3, L2, LI

WRITE (12, 4) K3, K2, Kl

WRITE (12, 6)NUP,MUP

WRITE (12, 7)NSOUR

WRITE(12,11) P0

WRITE (12, 8)

WRITE (12, 9)

DO 3888 I=l,NSOUR

WRITE(12, 10) I,WTSOUR(I) ,XSOUR(I) , YSOUR(I) , LXSOUR(I), LYSOUR(I)

3888 CONTINUE

1

5160

5165

1

5164

5166

5200

5100

70
APPENDIX A - TXYZ30 LISTING

IF (ITYPE.EQ.l) THEN

GOTO 7000

ELSE IF (ITYPE.EQ.2) THEN

WRITE (12, 32) NLINE

DO 6100 J=l, NLINE

WRITE (12, 33) J, XLINE (J) , LXLINE (J)

6100 CONTINUE

ELSE IF (ITYPE.EQ.3) THEN

WRITE (12, 34) NAREA

DO 6200 J=l, NAREA

WRITE (12, 35) J, XAREA(J) , LXAREA(J) , YAREA(J) , LYAREA(J)

6200 CONTINUE

END IF

GO TO 7 000

3999 WRITE (12, 88)

7000 STOP

END

Q* ************************************************ ********

C END OF THE MAIN PROGRAM

71
APPENDIX A - TXYZ30 LISTING

C FUNCTION UZERO LISTING

q* *****************************************************************************

C DESCRIPTION OF THE FUNCTION UZERO (N, M)

C THIS FUNCTION CALCULATES THE DOUBLE FOURIER COSINE TRANSFORM

C OF THE POWER DENSITY FUNCTION, U (X, Y) . THIS IS THE TRANSFORM

C FOR ALL OF THE HEAT SOURCES. THE ASSUMPTION IS MADE THAT THE

C POWER DENSITY IS UNIFORM AND EQUAL TO UNITY OVER THE SURFACE

C OF THE HEATING ELEMENTS. THAT IS,

C U(X, Y)=l (XSOUR(I)<=X<=XSOUR(I)+LXSOUR(I) AND

C YSOUR(I) <=Y<=YXOUR(I) +LYSOUR(I) ).

C U(X,Y)=0 OTHERWISE.

C UNDER THESE CONDITIONS, IT IS POSSIBLE TO ANALYTICALLY EVALUATE

C THE DOUBLE INTEGRAL FOR EACH HEATING ELEMENT. AS THE HEATING

C ELEMENTS ARE ASSUMED TO BE INDEPENDENT, THE CONTRIBUTION FROM

C EACH ELEMENT MAY BE ADDED TO OBTAIN THE U(N,M) FOR ALL.

C

C NONUNIFORM POWER DENSITIES MAY BE TAKEN CARE OF BY USING THE

C WEIGHTING FACTOR FOR EACH ELEMENT (READ IN THE MAIN PROGRAM)

C* ***************************************************************************** FUNCTION UZERO (N, M)

DIMENSION WTSOUR( 50), XSOUR(50), YSOUR(50)

REAL LXSOUR(50), LYSOUR(50), K3, K2, Kl, LX, LY, L3, L2, LI

COMMON K3, K2, Kl, LX, LY, L3, L2, LI

COMMON NSOUR, WTSOUR, XSOUR, YSOUR, LXSOUR, LYSOUR

PI=3. 14159265

UZERO=0 .

DO 500 1=1, NSOUR

IF(N.EQ.O) GO TO 100

TERMX = SIN (FLOAT (N) *P I* (XSOUR (I) + LXSOUR(I))/ LX)

1 - SIN (FLOAT (N) *PI *XSOUR ( I ) / LX)

TERMX=TERMX* LX/ (FLOAT (N) *PI)

GO TO 150

100 TERMX= LXSOUR (I)

150 IF(M.EQ.O) GO TO 200

TERMY = SIN (FLOAT (M) *PI* (YSOUR (I) + LYSOUR(I))/ LY)

1 - SIN (FLOAT (M) *PI*YSOUR (I) / LY)

TERMY=TERMY* LY/ (FLOAT (M) *PI)

GO TO 250

200 TERMY= LYSOUR (I)

250 TERMI=TERMX* TERMY

UZERO=UZERO+TERMI *WTSOUR ( I

500 CONTINUE

RETURN

END

72
APPENDIX A - TXYZ30 LISTING

q* **************************************************************************** * C FUNCTION FUNZ LISTING

C

C DESCRIPTION OF THE FUNCTION FUNZ (N, M, Z)

C THIS FUNCTION IS USED TO CALCULATE THE Z DEPENDENT PART OF THE

C FUNCTION USING THE S=0 VERSIONS OF EQUATIONS (15) -(17) OF KOKKAS

C THESE ARE USED IN CONJUNCTION WITH EQUATIONS (18) -(22) ALSO FOR

C S=0 (STEADY-STATE CONDITION) . THE SPECIFIC FORM OF FUNZ IS

C DETERMINED BY THE VALUE OF Z, I.E., IF Z FALLS IN THE TOP, MIDDLE,

C OR BOTTOM LAYER OF THE STRUCTURE. IN ADDITION, SPECIAL CARE IS

C TAKEN TO EVALUATE FUNZ FOR THE CASE WHERE GAMMA=0 AS THESE ARE

C SIMPLE, BUT THE COMPUTER DOES NOT KNOW HOW TO EVALUATE LIMITS.

C IN ADDITION, THE FORM OF THE THREE SOLUTIONS HAS BEEN CHANGED TO

C GET RID OF THE ARTIFICIAL OVERFLOW PROBLEMS COMING FROM THE

C HYPERBOLIC FUNCTIONS, COSH AND SINH, FOR THE CASES WHERE THE

C ARGUMENTS BECOME LARGE.

C

Q* ***************************************************************************** FUNCTION FUNZ(N,M, Z)

DIMENSION WTSOUR( 50) , XSOUR(50), YSOUR(50)

REAL LXSOUR(50), LYSOUR(50), K3, K2, Kl, LX, LY, L3, L2, LI

COMMON K3, K2, Kl, LX, LY, L3, L2, LI

COMMON NSOUR, WTSOUR, XSOUR, YSOUR, LXSOUR, LYSOUR

PI=3 . 14159265

GAMMA=SQRT ( (FLOAT (N) *PI/ LX) **2 + (FLOAT (M) *PI/ LY)**2)

VS=GAMMA* L3

VC=GAMMA* L2

VI=GAMMA* LI

VT=GAMMA* ( L3 + Z)

VM=GAMMA* ( L3+ L2+Z)

VB=GAMMA* ( L3+ L2+ Ll+Z)

BOTl=TANH (VS) *TANH (VI)

BOTl=BOTl+( Kl/ K2) *TANH (VS) *TANH (VC)

BOT2=( K2/ K3) *TANH (VI) *TANH (VC) + ( Kl/ K3)

GFUNC=1 . 0/ (BOT1+BOT2)

AZ=ABS (Z)

IF (AZ.GT. L3) GO TO 500

73
APPENDIX A - TXYZ30 LISTING ****************************************************

IF (GAMMA. EQ. 0.0) GO TO 100

TERMSl=TANH(VI)+( Kl/ K2)*TANH(VC)

TERMS2=( K2/ K3 ) *TANH (VT) * (TANH (VI ) *TANH (VC) + ( Kl/ K2)

TERMS=TERMS1+TERMS2

IF (Z.EQ.0.0) GO TO 90

IF (VS.GT.5.0.AND.VT.GT.5.0) GO TO 80

IF (VS.LT.5.0) GO TO 10

20 CONTINUE

IF (VT.LT.5.0) GO TO 30

C2=0 . 5*EXP (VT)

GO TO 40

30 C2=COSH(VT)

40 CONTINUE

FUNZ=GFUNC* TERMS * C

1 * C2 /GAMMA

RETURN

80 FUNZ=GFUNC* TERMS *EXP (GAMMA*Z) /GAMMA

RETURN

90 FUNZ=GFUNC* TERMS /GAMMA

RETURN

100 FUNZ=( L3+Z)+( K3/ K2)* L2+( K3/ Kl)* LI

RETURN

500 TOTAL= L3+ L2

IF (AZ .GT. TOTAL) GO TO 1500

10

CI' GO

CI

2 .0*EXP (-VS)

TO 2 0 1.0/COSH (VS)

74
APPENDIX A - TXYZ30 LISTING

Q* *****************************************************************************

C MIDDLE LAYER CALCULATION

C THIS IS THE MIDDLE LAYER CALCULATION WHICH IS DEFAULTED TO IF

C Z FALLS INTO THIS DOMAIN OF DEPTHS

Q* **************************************************************************** * IF (GAMMA. EQ. 0.0) GO TO 1000

TERMC=TANH (VI ) + ( Kl/ K2)*TANH(VM)

IF (VS.GT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0) GO TO 800

IF (VS.LT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0) GO TO 900

IF (VS.LT.5.0) GO TO 250

Cl=2 . 0*EXP (-VS)

GO TO 2 60

250 C1=1.0/COSH(VS)

2 60 CONTINUE

IF (VC.LT.5.0) GO TO 270

C2=2 . 0*EXP (-VC)

GO TO 280

270 C2=1.0/COSH(VC)

280 CONTINUE

IF (VM.LT.5.0) GO TO 320

C3=0.5*EXP (VM)

GO TO 330

320 C3=COSH(VM)

330 CONTINUE

FUNZ =GFUNC * TERMC* C

1 * C2 + C 3 / GAMMA

RETURN

800 FUNZ=GFUNC*TERMC*2 . 0*EXP (GAMMA*Z) /GAMMA

RETURN

900 FUNZ=GFUNC* TERMC *EXP (VS+GAMMA*Z) / (COSH (VS) * GAMMA)

RETURN

1000 FUNZ=( K3/ Kl)* Ll+( K3/ K2)*( L3+ L2+Z)

RETURN

75
APPENDIX A - TXYZ30 LISTING '************************************************

: BOTTOM LAYER CALCULATION

: THIS IS THE BOTTOM LAYER CALCULATION WHICH IS USED IF Z FALLS

J INTO THE BOTTOM LAYER

************************************

1500 IF (GAMMA. EQ . 0.0) GO TO 2000

IF (VS.GT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0.AND.

1 VB.GT.5.0) GO TO 1900

IF(VB. GT. 5.0. AND. VS. GT.5.0.AND.VC.LT. 5.0. AND. VI. LT. 5. 0)GO TO 2100

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO TO 2200

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.LT.5.0.AND.VI.GT.5.0)GO TO 2300

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO TO 2400

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.LT.5.0.AND.VI.GT.5.0)GO TO 2500

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0)GO TO 2600

IF (VS.LT.5.0) GO TO 1550

Cl=2 . 0*EXP (-VS)

GO TO 1560

1550 C1=1.0/COSH(VS)

1560 CONTINUE

IF (VC.LT.5.0) GO TO 157 0

C2=2 .0*EXP (-VC)

GO TO 1580

1570 C2=1.0/COSH(VC)

1580 CONTINUE

IF (VI.LT.5.0) GO TO 1590

C3=2 . 0*EXP (-VI)

GO TO 1600

1590 C3=1.0/COSH(VI)

1600 CONTINUE

C4=SINH (VB)

FUNZ =GFUNC * C

1 * C2 * C3 * C 4 / GAMMA

RETURN

1900 FUNZ=GFUNC * 4. 0*EXP (GAMMA* Z) /GAMMA

RETURN

2000 FUNZ=( K3/ Kl)*( L3+ L2+ Ll+Z)

RETURN

2100 FUNZ=GFUNC*EXP (VB-VS) / ( GAMMA*COSH (VC) *COSH(VI)

RETURN

2200 FUNZ=GFUNC*EXP (VB-VC) / (GAMMA* COSH (VS) *COSH (VI)

RETURN

2300 FUNZ=GFUNC*EXP (VB-VI) / (GAMMA*COSH (VS) *COSH(VC) ) RETURN

2400 FUNZ=GFUNC*2 . 0*EXP (VB-VS-VC) / (GAMMA*COSH (VI) ) RETURN

2500 FUNZ=GFUNC*2 . 0*EXP (VB-VS-VI) / (GAMMA* COSH (VC)

RETURN

2600 FUNZ=GFUNC*2 . 0*EXP (VB-VC-VI) / (GAMMA* COSH (VS)

RETURN

END

76
APPENDIX B - TML LISTING

Q* **************************************************************************** * C TML VERSION 1.0 - VERSION DATE 06/27/95

C

C THIS IS THE THERMAL MULTILAYER PROGRAM WHICH IS BASED UPON A THERMAL

C RECURSION RELATION WHICH IS SIMILAR TO THAT USED IN THE SOLUTION OF

C THE LAPLACE EQUATION FOR TWO-PROBE AND FOUR-PROBE RESISTANCE ANALYSIS.

C THIS PROGRAM CALCULATES THE STEADY-STATE SURFACE TEMPERATURE FOR A

C MULTILAYER RECTANGULAR STRUCTURE WITH AN ARBITRARY NUMBER OF LAYERS.

C THE CALCULATION IS PERFORMED FOR AN ARBITRARY NUMBER OF RECTANGULAR

C HEAT SOURCES/SINKS ON THE TOP SURFACE.

C THE TEMPERATURE MAY BE CALCULATED AS A POINT FUNCTION (X, Y) , A LINE

C AVERAGE, OR AN AREA AVERAGE ON THE TOP SURFACE.

C THE CALCULATION FOLLOWS FROM THE INPUT OF THE THICKNESSES AND THERMAL

C CONDUCTIVITIES OF ALL OF THE LAYERS IN THE STRUCTURE.

C - IT IS IMPORTANT TO EMPHASIZE THAT THE CALCULATION IS GENERAL FOR THE

C MULTILAYER STRUCTURE AND THE APPLICATION TO SEMICONDUCTOR STRUCTURES

C IS A SPECIAL CASE.

C

C THE STARTING POINT IS GIVEN IN EQUATIONS (13) -(23), WITH S=0

C (ZERO FREQUENCY, STEADY-STATE CONDITION), IN THE PAPER BY KOKKAS

.

C

C REFERENCES: THE ORIGINAL MATHEMATICAL ANALYSIS OF THE THREE-LAYER

C STRUCTURE WAS PERFORMED IN THE PAPER "THERMAL ANALYSIS

C OF MULTIPLE-LAYERED STRUCTURES" BY ACHILLES G. KOKKAS,

C IEEE TRANS. ELEC . DEV. VOL. ED-21, NO. 11, 674-681 (1974).

C THIS PAPER WAS DRAWN FROM HIS PHD THESIS: A. G. KOKKAS,

C "ANALYSIS AND DESIGN OF ELECTROTHERMAL INTEGRATED CIRCUITS,

C PH.D. THESIS, MIT, 1972.

C

C THE ORIGINAL FORTRAN IMPLEMENTATION OF THE STEADY STATE

C KOKKAS EQUATIONS IS CONTAINED IN THE TXYZ CODE AND

C WAS PRESENTED IN THE REPORT "SEMICONDUCTOR MEASUREMENT

C TECHNOLOGY: TXYZ: A PROGRAM FOR SEMICONDUCTOR IC THERMAL

C ANALYSIS" BY JOHN ALBERS, NBS SPECIAL PUBLICATION 400-7 6 C (APRIL 1984)

.

C

C VERSION 2.0 OF THE TXYZ CODE WAS PRESENTED IN THE REPORT

C "SEMICONDUCTOR MEASUREMENT TECHNOLOGY: VERSION 2 . 0 OF THE

C TXYZ THERMAL ANALYSIS PROGRAM: TXYZ20" BY JOHN ALBERS,

C NIST SPECIAL PUBLICATION 400-89 (JUNE 1992)

.

C

C VERSION 3.0 OF THE TXYZ CODE AND THE ACCOMPANYING THERMAL

C MULTILAYER CODES ARE DISCUSSED IN THIS REPORT

C "SEMICONDUCTOR MEASUREMENT TECHNOLOGY: HOTPAC : PROGRAMS

C FOR THERMAL ANALYSIS INCLUDING VERSION 3 . 0 OF THE TXYZ

C PROGRAM, TXYZ30, AND THE THERMAL MULTILAYER PROGRAM, TML"

C BY JOHN ALBERS, NIST SPECIAL PUBLICATION 400-96.

C

C THE REVIEW AND APPLICATION OF THE RECURSION RELATION

C TECHNIQUE FOR THE ANALYSIS OF THE LAPLACE EQUATION USED IN

C TWO-PROBE AND FOUR-PROBE RESISTANCE ARE CONTAINED IN THE

C REPORT "SEMICONDUCTOR MEASUREMENT TECHNOLOGY: A COLLECTION

C OF COMPUTER PROGRAMS FOR TWO-PROBE RESISTANCE (SPREADING

C RESISTANCE) AND FOUR-PROBE RESISTANCE CALCULATIONS, RESPAC

C BY JOHN ALBERS AND HARRY L. BERKOWITZ, NIST SPECIAL

C PUBLICATION 400-91, 1993.

C

77
APPENDIX B - TML LISTING

C THE THERMAL RECURSION RELATION TECHNIQUE HAS BEEN DISCUSSED

C IN THE PAPERS: "AN EXACT SOLUTION OF THE STEADY-STATE

C SURFACE TEMPERATURE FOR A GENERAL MULTILAYER STRUCTURE"

C BY JOHN ALBERS, PROCEEDINGS TENTH IEEE SEMI-THERM SYMPOSIUM,

C PP. 129-137, 1994 AND

C "AN EXACT RECURSION RELATION SOLUTION FOR THE STEADY-STATE

C SURFACE TEMPERATURE OF A GENERAL MULTILAYER STRUCTURE"

C BY JOHN ALBERS, IEEE TRANSACTIONS ON COMPONENTS, PACKAGING

C AND MANUFACTURING TECHNOLOGY - PART A, VOL. 18, NO. 1,

C PP. 31-38, 1995.

C

Q* **************************************************************************** * C IN THE PRESENT FORM, THE PROGRAM ALLOWS UP TO 500 TERMS TO BE INCLUDED

C IN BOTH THE N SUM (ALONG X) AND THE M SUM (ALONG Y) . TO GO BEYOND

C THIS NUMBER, SUBSTITUTE THE FOLLOWING TWO LINES WITH THE APPROPRIATE

C VALUES OF NX AND MY FOR THE FIRST TWO DIMENSION STATEMENTS - MAKE, SURE TO

C REMOVE THE COMMENTS FROM THE NEW LINES AND COMMENT OUT THE REPLACED LINES.

C. DIMENSION X(100), Y(100), Z(100), COSYT(MY)

C DIMENSION ARUZER (NX, MY) , ARFUNZ (NX, MY)

C ALSO REPLACE THE TESTING LINE ABOVE THE LINE LABELLED 110 WITH

C IF (NUP.GT.NX.OR.MUP.GT.MY) GO TO 3999

DIMENSION X (100) , Y(100), COSYT(500)

DIMENSION ARUZER(500, 500) , ARFUN0 (500, 500)

DIMENSION WTSOUR (50) , XSOUR(50), YSOUR(50)

REAL LXSOUR(50), LYSOUR(50), K(20), LX, LY, L(20)

DIMENSION XLINE (30) , XAREA (30) , YAREA (30)

REAL LXLINE (30) , LXAREA (30) , LYAREA (30)

COMMON LX, LY, NLAY, K, L

COMMON NSOUR, WTSOUR, XSOUR, YSOUR, LXSOUR, LYSOUR

PI=3. 14159265

C INPUT DATA IS READ FROM 10 AND OUTPUT GOES TO 12. THESE ARE WRITTEN IN

C LOWER CASE. MANY OPERATING SYSTEMS TAKE UPPER AND LOWER CASE AS

C EQUIVALENT. THE UNIX OPERATING SYSTEM VIEWS THE UPPER AND LOWER

C CASE-NAMED FILES AS DIFFERENT FILES.

C UNIX DEFAULTS TO THE LOWER CASE WHICH IS USED HERE.

open ( unit=10, file=' input .dat' , status = 'unknown')

open ( unit=12, file=' output .dat' , status = 'unknown' ) 2 FORMAT (IX, ' STEADY-STATE THERMAL MULTILAYER CALCULATION'/

1' USING THE THERMAL RECURSION RELATION IN KOKKAS EQUATIONS')

3 FORMAT (IX, ' THICKNESSES AND THERMAL CONDUCTIVITIES OF LAYERS')

4 FORMAT (IX, 'NUMBER OF LAYERS IN STRUCTURE^ , 12

5 FORMAT (IX, ' L' , 12, ' = ',F10.5,' K' , 12, ' = ' , F10 . 5)

6 FORMAT (IX, 'UPPER SUMMATION LIMITS ' , 2X, ' NUP=' , 15,

1 ' MUP=',I5)

7 FORMAT (IX, 'NUMBER OF HEAT SOURCES=' , 15)

8 FORMAT (IX, 'WEIGHTS, COORDINATES, LENGTHS, WIDTHS OF SOURCES')

9 FORMAT ( IX, ' HEAT SOURCE' , 4X, ' WTSOUR' , 8X, ' XSOUR' , 8X, ' YSOUR' , 7X,

1 ' LXSOUR' , 7X, ' LYSOUR'

10 FORMAT (4X, 13, 5X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5, 3X,F10.5)

11 FORMAT (IX, 'POWER DENSITY=' ,F1 1.6)

27 FORMAT ( IX, ' LX= ',F7.2,3X, ' LY= ',F7.2)

22 FORMAT (IX, F12 . 4, 2X, F12 . 4, 2X, F12 . 4, 2X, F12 .4)

31 FORMAT (IX, 617)

32 FORMAT (IX, 'NUMBER OF LINES=' , 13)

33 FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5)

34 FORMAT (IX, 'NUMBER OF AREAS=' , 13)

78
APPENDIX B - TML LISTING

35 FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5)

40 FORMAT (IX, 'POINT FUNCTION EVALUATION OF SURFACE T (X, Y) '

) 41 FORMAT (IX, ' LINE AVERAGE EVALUATION OF SURFACE TEMPERATURE')

42 FORMAT ( IX, ' AREA AVERAGE EVALUATION OF SURFACE TEMPERATURE')

43 FORMAT (6X, 'X' , 12X, ' Y' , 11X, ' T (X, Y) '

)

44 FORMAT(1X,F10.5, 3X,F10.5, 3X,F10.5)

45 FORMAT (IX, ' LINE #' , 5X, ' XLINE' , 7X, ' LXLINE' , 8X, ' Y' , 9X, ' AVE TEMP')

4 6 FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X,F10.5)

47 FORMAT (IX, 'AREA #' , 5X, ' XAREA' , 7X, ' LXAREA' , 8X, ' YAREA' , 7X,

1 ' LYAREA' , 3X, ' AVE TEMP'

)

48 FORMAT (IX, 13, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5, IX, F10 . 5)

51 FORMAT (IX, 14, 3X, 14)

52 FORMAT (IX, F10 . 5, 3X,F10.5)

53 FORMAT (IX, II)

54 FORMAT (IX, 14, 3X, F10 . 5, 3X, F10 . 5)

55 FORMAT (IX, F10 . 5)

56 FORMAT ( IX, 12 , 3X, F10 . 5

)

57 FORMAT(1X,F10.5, 3X, F10 . 5, 3X, F10 . 5, 3X, F10 . 5)

88 FORMAT (IX, ' YOUR UPPER LIMIT OF SUMMATION IS TOO LARGE. TRY AGAIN')

q* *****************************************************************************

C

C INPUT SECTION

C THE FOLLOWING VARIABLES ARE READ BY THE PROGRAM FROM FOR010

C

C I TYPE TYPE OF ANALYSIS

C =1 FOR POINT FUNCTION, T(XY, 0), EVALUATION

C =2 FOR LINE AVERAGE, <T(Y,0)>, EVALUATION

C =3 FOR AREA AVERAGE, <T(0)>, EVALUATION

C

C LX X DIMENSION OF THE MULTILAYER STRUCTURE

C LY Y DIMENSION OF THE MULTILAYER STRUCTURE

C NLAY NUMBER OF LAYERS IN THE STRUCTURE (UP TO 2 0 LAYERS

C BUT THE CODE CAN GO BEYOND THIS BY EDITING THE

C APPROPRIATE DIMENSION STATEMENTS IN THE MAIN PROGRAM

C AND SUBROUTINES/FUNCTIONS)

C

C THE NEXT NLAY LINES CONTAIN THE THICKNESSES AND THERMAL

C CONDUCTIVITIES OF THE N-LAYERS BEGINNING AT THE SURFACE AND GOING

C TO THE BOTTOM LAYER (CODE TAKES CARE OF NUMBERING FROM NLAY, ...,1)

C

C L(I), K (I) -THICKNESS AND THERMAL CONDUCTIVITY OF I-TH LAYER

C NUP UPPER LIMIT OF N SUM, X DIRECTION

C MUP UPPER LIMIT OF M SUM, Y DIRECTION

C

C IEDGEX INDEX FOR GENERATING THE VALUES OF X TO BE USED

C =1 TO READ DATA FOR FIXED INCREMENT X VALUES

C =2 TO READ IN ARRAY OF X VALUES OF NONFIXED INCREMENT)

C IF IEDGEX=1 THEN READ THE THREE VARIABLES (ON SAME LINE)

C ILX THE NUMBER OF X VALUES TO BE USED

C XI THE VALUE OF THE FIRST POINT IN X

C STEPX (THE INCREMENT IN X)

C IF IEDGEX=2 THEN READ THE VARIABLE AND ARRAY (ONE PER LINE)

C ILX THE NUMBER OF X VALUES TO BE USED

C X(I) THE ARRAY OF X VALUES (1=1, ILX)

C

79
APPENDIX B - TML LISTING

C IEDGEY INDEX FOR GENERATING THE VALUES OF Y TO BE USED

C =1 TO READ DATA FOR FIXED INCREMENT Y VALUES

C =2 TO READ IN ARRAY OF X VALUES OF NONFIXED INCREMENT)

C IF IEDGEY=1 THEN READ THE THREE VARIABLES (ON SAME LINE)

C ILY THE NUMBER OF Y VALUES TO BE USED

C Yl THE VALUE OF THE FIRST POINT IN Y

C STEPY (THE INCREMENT IN Y)

C IF IEDGEY=2 THEN READ THE VARIABLE AND ARRAY (ONE PER LINE)

C ILY THE NUMBER OF Y VALUES TO BE USED

C Y(I) THE ARRAY OF Y VALUES (1=1, ILY)

C

C NOTE: THERE IS NO INPUT FOR THE Z VALUES AS Z=0 THROUGHOUT THE

C CALCULATION

C

C IMPORTANT: THE POINT FUNCTION EVALUATION OF THE TEMPERATURE IS

C THE MOST ELEMENTAL VERSION. IN ORDER TO SIMPLIFY THE

C CONSTRUCTION OF THE DATA FILES FOR LINE AND AREA

C AVERAGES, THE PROGRAM EXPECTS TO SEE THE ABOVE IEDGEX,

C AND IEDGEY DATA. THIS IS READ EVEN IF IT IS NOT

C USED FOR THE AVERAGE VERSIONS. HOWEVER, THE LINE AND

C AREA INFORMATION MAY THEN BE SIMPLY APPENDED TO THE END

C OF THE DATA FILE IN ORDER TO RUN THESE VERSIONS. SEE THE

C SAMPLE I/O FILES FOR AN EXAMPLE OF THIS.

C

C NSOUR NUMBER OF HEAT SOURCES (UP TO 50)

C P0 POWER DENSITY

C

C THE NEXT NSOUR LINES READ THE FOLLOWING INFORMATION FOR THE

C HEAT SOURCES (WITH ALL THE INFORMATION FOR EACH OF THE ELEMENTS

C ON A SINGLE LINE)

C

C WTSOUR (I) -WEIGHTING FACTOR OF I-TH SOURCE

C (POSITIVE FOR SOURCE, NEGATIVE FOR SINK)

C XSOUR(I) —X COORDINATE OF ORIGIN OF I-TH SOURCE

C LXSOUR (I) -LENGTH ALONG X DIRECTION OF I-TH SOURCE

C YSOUR(I) —Y COORDINATE OF ORIGIN OF I-TH SOURCE

C LYSOUR (I) -LENGTH ALONG Y DIRECTION OF I-TH SOURCE

C

C IF ITYPE=1, THE POINT FUNCTION CALCULATION CONTINUES WITH THE ABOVE

C SET OF X,Y VALUES

C IF ITYPE=2, THE LINE AVERAGE CALCULATION READS THE FOLLOWING:

C NLINE THE NUMBER OF LINE SEGMENTS TO DO THE AVERAGE

C THE NEXT NLINE LINES THEN READ

C XLINE(J), LXLINE(J) —THE LOCATION AND LENGTH OF THE

C J-TH LINE ELEMENT

C IF ITYPE=3, THE AREA AVERAGE CALCULATION READS THE FOLLOWING:

C NAREA THE NUMBER OF AREA SEGMENTS TO DO THE AVERAGE

C THE NEXT NAREA LINES THEN READ

C XAREA (J) , LXAREA ( J) , YAREA (J) , LYAREA ( J) -THE LOCATION

C AND LENGTHS OF THE J-TH AREA ELMENT

80
APPENDIX B - TML LISTING

C READ ITYPE

READ (10,*) ITYPE

C READ LX AND LY (THE X AND Y DIMENSIONS OF THE RECTANGULAR STRUCTURE)

READ (10,*) LX, LY

C READ THE NUMBER OF LAYERS IN THE STRUCTURE

READ (10,*) NLAY

C READ LAYER THICKNESS AND THERMAL CONDUCTIVITY FROM SURFACE TO BOTTOM

DO 101 I=NLAY, 1,-1

READ (10, * ) L(I),K(I)

101 CONTINUE

C READ NUP AND MUP (UPPER LIMIT OF THE SUMMATION OVER THE INDEX N (X-DIR)

C UPPER LIMIT OF THE SUMMATION OVER THE INDEX M (Y-DIR)

)

C NUP AND MUP MUST BE LESS THAN OR EQUAL TO THE DIMENSIONALITY OF

READ (10, *) NUP, MUP

IF (NUP. GT. 500. OR. MUP. GT. 500) GO TO 3999

READ (10, *) IEDGEX

GOTO (110, 115) IEDGEX

110 READ (10, *) ILX, XI, STEPX

DO 111 1=1, ILX

X(I)=X1+ (1-1) *STEPX

111 CONTINUE

GOTO 119

115 READ (10,*) ILX

DO 116 1=1, ILX

READ (10, *) X (I)

116 CONTINUE

119 READ (10, *) IEDGEY

GOTO (120,125) IEDGEY

120 READ (10, *) ILY, Yl, STEPY

DO 121 1=1, ILY

Y(I)=Y1+(I-1) *STEPY

121 CONTINUE

GOTO 12 9 125 READ (10,*) ILY

DO 126 1=1, ILY

READ (10, *) Y(I)

12 6 CONTINUE

129 Z=0.0

C READ THE NUMBER OF HEAT SOURCES AND THE POWER DENSITY

C NOTE-POWER DENSITY IS MULTIPLICATIVE FACTOR USUALLY SET EQUAL TO UNITY

C P0 IS THE POWER DENSITY, ASSUMED UNIFORM FOR ALL HEATERS

C NSOUR IS THE TOTAL NUMBER OF HEATING ELEMENTS ON THE SURFACE OF THE

C THE TOP LAYER (UP TO 50)

139 READ (10, *)NSOUR, P0

C THE NEXT LOOP READS IN THE COORDINATES OF THE ORIGIN OF THE

C HEATING ELEMENTS ALONG WITH THEIR LENGTHS AND WIDTHS

C THE WEIGHTING FACTOR IS ALSO ENTERED (THIS IS REAL, NONINTEGER)

C WTSOUR(I) IS THE WEIGHTING FACTOR FOR THE I-TH HEATER ELEMENT

C XSOUR(I) IS THE X COORDINATE OF THE ORIGIN OF I-TH HEATER ELEMENT

C LXSOUR(I) IS THE LENGTH OF THE I-TH HEATER ALONG THE X DIRECTION

C YSOUR(I) IS THE Y COORDINATE OF THE ORIGIN OF I-TH HEATER ELEMENT

C LYSOUR(I) IS THE LENGTH OF THE I-TH HEATER ALONG THE Y DIRECTION

DO 140 1=1, NSOUR

READ (10, *)WTSOUR(I) , XSOUR(I) , LXSOUR ( I ) , YSOUR ( I ) , LYSOUR(I)

140 CONTINUE

81
APPENDIX B - TML LISTING

C THIS IF .. .THEN. . .ELSE IF CONSTRUCTION IS USED TO READ IN THE DATA FOR

C THE LINES OR AREAS TO BE CONSIDERED IN THE CALCULATION

C IF THE ANALYSIS IS FOR A POINT FUNCTION, THEN GO OUT OF THE IF. THEN. ELSE

IF (ITYPE.EQ.l) THEN

GOTO 17 0 C IF THE ANALYSIS IS FOR A LINE AVERAGE, THEN READ THE NUMBER OF LINE

C SEGMENTS AND THEIR LOCATION

ELSE IF ( I TYPE .EQ . 2 ) THEN

C READ THE NUMBER OF LINE SEGMENTS

READ (10, * ) NLINE

DO 150 J=l, NLINE

C READ THE ORIGIN AND LENGTH OF EACH LINE SEGMENT

C NOTE THAT THE AVERAGE IS ALONG THE X-DIRECTION FOR GIVEN Y VALUES

C TO DO THE AVERAGE ALONG THE Y-DIRECTION FOR GIVEN X, SIMPLY ROTATE

C THE STRUCTURE BY 90 DEGREES AND USE THE CORRESPONDING NEW X' S AND Y'

C XLINE(J) IS THE X COORDINATE OF THE ORIGIN OF J-TH LINE ELEMENT

C LXLINE(J) IS THE LENGTH OF THE J-TH LINE ALONG THE X DIRECTION

READ (10, * ) XLINE (J) , LXLINE (J)

150 CONTINUE

C IF THE ANALYSIS IS FOR AN AREA AVERAGE, THEN READ THE NUMBER OF AREAS

C AND THEIR LOCATIONS

ELSE IF (ITYPE.EQ.3) THEN

C READ THE NUMBER OF AREAS

READ (10, * ) NAREA

DO 160 J=l, NAREA

C READ THE FOLLOWING

C XAREA(J) IS THE X COORDINATE OF THE ORIGIN OF J-TH AREA ELEMENT

C LXAREA(J) IS THE LENGTH OF THE J-TH AREA ALONG THE X DIRECTION

C YAREA(J) IS THE Y COORDINATE OF THE ORIGIN OF J-TH AREA ELEMENT

C LYAREA(J) IS THE LENGTH OF THE J-TH AREA ALONG THE Y DIRECTION

READ (10, * ) XAREA (J) , LXAREA (J) , YAREA ( J) , LYAREA ( J)

160 CONTINUE

END IF

17 0 CONTINUE

£******************************************************************************

C END OF DATA INPUT SECTION

C BEGIN CALCULATION OF T (X, Y, Z)

C THE SUBROUTINES USED IN THE CALCULATION ARE:

C 1) UZERO(N,M) - CALCULATES THE FOURIER COSINE TRANSFORM OF THE

C FUNCTION, U(X,Y), THE POWER DENSITY FUNCTION FOR ALL OF THE

C HEAT SOURCES.

C 2) FUN0(N,M) - CALCULATES THE FOURIER COEFFICIENTS FOR THE LAYERED

C STRUCTURE USING THE THERMAL RECURSION RELATION

Q* *****************************************************************************

P04LK = 4.0 * P0 / ( LX* LY)

PILX = PI / LX

PILY = PI / LY

C CALCULATE THE FOURIER COMPONENTS OF THE HEAT SOURCES, U(N,M)

q* *****************************************************************************

DO 300 MM=1,MUP

M = MM - 1 DO 2 50 NN=1,NUP

N = NN - 1 ARUZER (NN, MM) =UZERO (N, M)

250 CONTINUE

300 CONTINUE

£*********************************+++++*+++*+*+++++++++++++*++++*+++++++*******

C END OF U(N,M) CALCULATION AND BEGINNING OF MAJOR LOOP FOR Z

82
APPENDIX B - TML LISTING

q* *********************************************************************

C CALCULATE THE FUNO(N,M) FROM THE RECURSION RELATION

Q* *********************************************************************

DO 400 MM=1,MUP

M = MM - 1 DO 350 NN=1,NUP

N = NN - 1 ARFUNO (NN, MM) =FUN0 (N, M) *ARUZER (NN, MM)

350 CONTINUE

400 CONTINUE

C THE FOLLOWING IF. THEN. ELSE CONSTRUCTION OPERATES ACCORDING TO THE

C TYPE OF ANALYSIS TO BE USED.

C FOR THE POINT FUNCTION ANALYSIS, THE 3000 LOOP IS USED

C FOR THE LINE AVERAGE ANALYSIS, THE 4000 LOOP IS USED

C FOR THE AREA AVERAGE ANALYSIS, THE 5000 LOOP IS USED

C THIS PORTION DOES THE POINT FUNCTION CALCULATION

IF (ITYPE.EQ.l) THEN

C BEGINNING OF POINT FUNCTION ANALYSIS

WRITE (12, 2)

WRITE (12, 40)

WRITE (12, 43)

DO 3000 IY=1, ILY

DO 3100 MM=1,MUP

M = MM - 1 COSYT (MM) =COS (FLOAT (M) *Y(IY) *PILY)

3100 CONTINUE

DO 3000 IX=1, ILX

SUM=0 .

DO 3300 MM=1,MUP

M = MM - 1 DO 3200 NN=1,NUP

N = NN - 1 NDN=0

NDM=0

IF (N.EQ.0) NDN=1

IF (M.EQ.0) NDM=1

TOP = ARFUNO (NN, MM) * COS (FLOAT (N) *X ( IX) *PILX) * COSYT (MM)

BOTTOM= (NDN+1) * (NDM+1)

TSUM=TOP /BOTTOM

SUM=SUM+TSUM

32 00 CONTINUE

3300 CONTINUE

TEMP = P04LK * SUM/K (NLAY)

WRITE (12, 44)X(IX),Y(IY), TEMP

3000 CONTINUE

C END OF POINT FUNCTION CALCULATION

ELSE IF (I TYPE . EQ . 2 ) THEN

C THIS PORTION DOES THE LINE AVERAGE CALCULATION

WRITE (12, 2)

WRITE (12, 41)

WRITE (12, 45)

83
APPENDIX B - TML LISTING

DO 4000 IY=1,ILY

DO 4010 MM=1,MUP

M = MM - 1 COSYT (MM) =COS (FLOAT (M) *Y(IY) *PILY)

4010 CONTINUE

DO 4000 J=1,NLINE

SUM=0 .

DO 4200 MM=1,MUP

M = MM - 1 DO 4100 NN=1,NUP

N = NN - 1 NDN=0

NDM=0

IF (N.EQ.0) NDN=1

IF (M.EQ.0) NDM=1

IF(N.EQ.O) GO TO 4160

TERMX=SIN (FLOAT (N) *PI* (XLINE (J) +LXLINE (J) ) /LX)

1 -SIN(FLOAT(N)*PI*XLINE(J)/LX)

TERMX=TERMX*LX/ (FLOAT (N) *PI)

GO TO 4165

4160 TERMX=LXLINE (J)

4165 CONTINUE

TOP = ARFUN0 (NN, MM) * COSYT (MM) * TERMX

BOTTOM= (NDN+1) * (NDM+1)

TSUM=TOP /BOTTOM

SUM=SUM+TSUM

4100 CONTINUE

TEMP = P04LK * SUM / (LXLINE (J) * K(NLAY))

42 00 CONTINUE

WRITE (12, 46) J, XLINE (J) , LXLINE (J) , Y(IY) , TEMP

4 000 CONTINUE

C END OF LINE AVERAGE PORTION OF THE CODE

ELSE IF (ITYPE.EQ.3) THEN

C THIS PART PERFORMS THE AREA AVERAGE CALCULATION

WRITE (12, 2)

WRITE (12, 42)

WRITE (12, 47)

DO 5000 J= 1 , NARE

SUM=0.0

DO 5100 MM=1,MUP

M = MM - 1 DO 5200 NN=1,NUP

N = NN - 1 NDN=0

NDM=0

IF (N.EQ.0) NDN=1

IF (M.EQ.0) NDM=1

AREA=0 .

IF (N.EQ.0) GO TO 5160

TERMX = SIN (FLOAT (N) *PI* (XAREA (J) + LXAREA ( J) ) / LX)

1 - SIN (FLOAT (N) *PI*XAREA (J) / LX)

TERMX=TERMX* LX/ (FLOAT (N) *PI)

GO TO 5165

84
APPENDIX B - TML LISTING

5160 TERMX= LXAREA(J)

5165 IF(M.EQ.O) GO TO 5164

TERMY = SIN (FLOAT (M) *PI* (YAREA (J) + LYAREA ( J) ) / LY)

1 - SIN (FLOAT (M) *P I* YAREA (J) / LY)

TERMY=TERMY* LY/ (FLOAT (M) *P I

)

GO TO 5166

5164 TERMY= LYAREA (J)

5166 TERMI=TERMX*TERMY

AREA=TERMI

TOP = ARFUNO (NN, MM) * AREA

BOTTOM= (NDN+1) * (NDM+1)

TSUM=TOP /BOTTOM

SUM=SUM+TSUM

5200 CONTINUE

TEMP = P04LK * SUM / (LXAREA (J) *LYAREA (J) *K (NLAY)

)

5100 CONTINUE

WRITE (12, 48) J, XAREA (J) , LXAREA (J) , YAREA (J) , LYAREA (J) , TEMP

5000 CONTINUE

END IF

WRITE (12,27) LX, LY

WRITE (12, 4) NLAY

WRITE (12, 3)

DO 3001 I=NLAY, 1,-1

WRITE (12, 5) I,L(I), I,K(I)

3001 CONTINUE

WRITE (12, 6)NUP,MUP

WRITE (12, 11) P0

WRITE (12, 7)NSOUR

WRITE (12, 8)

WRITE (12, 9)

DO 3888 I=l,NSOUR

WRITE (12, 10) I, WTSOUR(I) , XSOUR(l) , YSOUR(I) , LXSOUR(I), LYSOUR(I)

3888 CONTINUE

IF (ITYPE.EQ. 1) THEN

GOTO 7000

ELSE IF (ITYPE.EQ. 2) THEN

WRITE (12, 32)NLINE

DO 6100 J=1,NLINE

WRITE (12, 33) J, XLINE (J) , LXLINE (J)

6100 CONTINUE

ELSE IF (ITYPE.EQ. 3) THEN

WRITE (12, 34) NAREA

DO 6200 J=l, NAREA

WRITE (12, 35) J, XAREA (J) , LXAREA (J) , YAREA (J) , LYAREA (J)

6200 CONTINUE

END IF

GO TO 7000

3999 WRITE (12, 88)

7000 STOP

END ^*********************************^

; END OF THE MAIN PROGRAM

85
APPENDIX B - TML LISTING

: FUNCTION UZERO LISTING

DESCRIPTION OF THE FUNCTION UZERO (N,M)

THIS FUNCTION CALCULATES THE DOUBLE FOURIER COSINE TRANSFORM

OF THE POWER DENSITY FUNCTION, U(X,Y). THIS IS THE TRANSFORM

FOR ALL OF THE HEAT SOURCES. THE ASSUMPTION IS MADE THAT THE

POWER DENSITY IS UNIFORM AND EQUAL TO UNITY OVER THE SURFACE

OF THE HEATING ELEMENTS. THAT IS,

U(X,Y)=1 (XSOUR(I) <=X<=XSOUR(I) +LXSOUR(I) AND

YSOUR ( I ) <=Y<=YXOUR ( I ) +LYSOUR ( I ) )

U(X,Y)=0 OTHERWISE.

UNDER THESE CONDITIONS, IT IS POSSIBLE TO ANALYTICALLY EVALUATE

THE DOUBLE INTEGRAL FOR EACH HEATING ELEMENT. AS THE HEATING

ELEMENTS ARE ASSUMED TO BE INDEPENDENT, THE CONTRIBUTION FROM

EACH ELEMENT MAY BE ADDED TO OBTAIN THE U (N, M) FOR ALL.

C

C

C

C

C

c

c

c

c

c

cc

cc

c

c

100

150

200

250

500

NONUNIFORM POWER DENSITIES MAY BE TAKEN CARE OF BY USING THE

WEIGHTING FACTOR FOR EACH ELEMENT (READ IN THE MAIN PROGRAM)

***************+********************

FUNCTION UZERO (N,M)

DIMENSION WTSOUR( 50) , XSOUR(50), YSOUR(50)

REAL LXSOUR(50), LYSOUR(50), K(20), LX, LY, L(20)

COMMON LX, LY, NLAY, K, L

COMMON NSOUR, WTSOUR, XSOUR, YSOUR, LXSOUR, LYSOUR

PI=3 . 14159265

UZERO=0 .

DO 500 1=1, NSOUR

IF(N.EQ.O) GO TO 100

TERMX = SIN (FLOAT (N) *PI* (XSOUR (I) + LXSOUR(I))/ LX)

1 - SIN (FLOAT (N) *PI*XSOUR ( I ) / LX)

TERMX=TERMX * LX/ (FLOAT (N) *PI

GO TO 150

TERMX= LXSOUR (I)

IF(M.EQ.O) GO TO 200

TERMY = SIN (FLOAT (M) *PI* (YSOUR (I) + LYSOUR(I))/ LY)

1 - SIN (FLOAT (M) *P I* YSOUR (I) / LY)

TERMY=TERMY * LY/ (FLOAT (M) *PI)

GO TO 250

TERMY= LYSOUR (I)

TERMI=TERMX* TERMY

UZERO=UZERO+TERMI *WTSOUR ( I

CONTINUE

RETURN

END *********************************************

86
APPENDIX B - TML LISTING

C FUNCTION FUNO LISTING

C

C DESCRIPTION OF THE FUNCTION FUN0(N,M)

C

Q* ***************************************************************************** FUNCTION FUN0(N,M)

DIMENSION WTSOUR(50), XSOUR(50), YSOUR(50)

REAL LXSOUR(50), LYSOUR(50), K(20), LX, LY, L(20)

COMMON LX, LY, NLAY, K, L

COMMON NSOUR, WTSOUR, XSOUR, YSOUR, LXSOUR, LYSOUR

PI=3. 14159265

GAMMA=SQRT ( (FLOAT (N) *PI/ LX) **2 + (FLOAT (M) *PI/ LY)**2)

IF (GAMMA. EQ.O) GO TO 120

C FOR GAMMA. GT.O, START WITH BOTTOM LAYER

FUN0=TANH (GAMMA*L (1)

IF (NLAY.EQ.l) GO TO 110

C THE 100 LOOP DOES THE RECURSION RELATION FOR NONZERO GAMMA

DO 100 1=2, NLAY

TOP=FUN0*K(I) +K(I-1) *TANH (GAMMA*L (I)

BOT=K (I— 1) +K (I) *FUN0*TANH (GAMMA*L (I)

FUN0=TOP/BOT

100 CONTINUE

110 FUN0=FUN0 /GAMMA

RETURN

120 FUN0=L(1)

IF (NLAY.EQ.l) GO TO 130

C THE 125 LOOP DOES THE RECURSION RELATION FOR GAMMA EQUAL TO ZERO

DO 125 1=2, NLAY

FUN0= (FUN0*K(I)+K(I-1) *L(I) ) /K(I-l)

125 CONTINUE

130 RETURN

END

87



Technical Publications

Periodical

Journal of Research of the National Institute of Standards and Technology—Reports NIST research and development in those disciplines of the physical and engineering sciences in which the Institute is active. These include physics, chemistry, engineering, mathematics, and computer sciences. Papers cover a broad range of subjects, with major emphasis on measurement methodology and the basic technology

underlying standardization. Also included from time to time are survey articles on topics closely related to

the Institute's technical and scientific programs. Issued six times a year.

Nonperiodicals

Monographs—Major contributions to the technical literature on various subjects related to the

Institute's scientific and technical activities. Handbooks—Recommended codes of engineering and industrial practice (including safety codes) devel - oped in cooperation with interested industries, professional organizations, and regulatory bodies. Special Publications—Include proceedings of conferences sponsored by NIST, NIST annual reports, and

other special publications appropriate to this grouping such as wall charts, pocket cards, and bibliographies.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical

properties of materials, compiled from the world's literature and critically evaluated. Developed under a worldwide program coordinated by NIST under the authority of the National Standard Data Act (Public

Law 90-396). NOTE: The Journal of Physical and Chemical Reference Data (JPCRD) is published

bimonthly for NIST by the American Chemical Society (ACS) and the American Institute of Physics (AIP).

Subscriptions, reprints, and supplements are available from ACS, 1155 Sixteenth St., NW, Washington, DC

20056.

Building Science Series—Disseminates technical information developed at the Institute on building

materials, components, systems, and whole structures. The series presents research results, test methods, and

performance criteria related to the structural and environmental functions and the durability and safety

characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in their treatment of

a subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often serve as a vehicle for final reports of work performed at NIST under the sponsorship of

other government agencies.

Voluntary Product Standards—Developed under procedures published by the Department of Commerce

in Part 10, Title 15, of the Code of Federal Regulations. The standards establish nationally recognized

requirements for products, and provide all concerned interests with a basis for common understanding of

the characteristics of the products. NIST administers this program in support of the efforts of private-sector

standardizing organizations.

Order the following NIST publications—FIPS and NISTIRs—from the National Technical Information

Service, Springfield, VA 22 1 61. Federal Information Processing Standards Publications (FIPS PUB)—Publications in this series

collectively constitute the Federal Information Processing Standards Register. The Register serves as the

official source of information in the Federal Government regarding standards issued by NIST pursuant to the Federal Property and Administrative Services Act of 1949 as amended, Public Law 89-306 (79 Stat. 1127), and as implemented by Executive Order 11717 (38 FR 12315, dated May 11, 1973) and Part 6 of

Title 15 CFR (Code of Federal Regulations).

NIST Interagency Reports (NISTIR)—A special series of interim or final reports on work performed by

NIST for outside sponsors (both government and nongovernment). In general, initial distribution is handled

by the sponsor; public distribution is by the National Technical Information Service, Springfield, VA 22161,

in paper copy or microfiche form.
occ o ta u.E

v> ~°

i- So S

ca a oU IS s o o N u

cn e <u Q

ess

D

a |S

u

r c 5

usin 'Ca

aa ational I

aithersbi :ial B

lty

for £ BV 3 Z O O a.