NBS

PUBLICATIONS

,^lT OF

O NBS SPECIAL PUBLICATION 400-76

U.S. DEPARTMENT OF COMMERCE/National Bureau of Standards

Semiconductor Measurement Technology:

TXYZ: A Program

for Semiconductor

IC Thermal Analysis

QC-

. U57

/{GO-76

193'!
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards' was established by an act of Congress on March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and technology

and facilitate their effective application for public benefit. To this end, the Bureau conducts

research and provides: (1) a basis for the Nation's physical measurement system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau's technical work is per- formed by the National Measurement Laboratory, the National Engineering Laboratory, and

the Institute for Computer Sciences and Technology.

THE NATIONAL MEASUREMENT LABORATORY provides the national system of

physical and chemical and materials measurement; coordinates the system with measurement

systems of other nations and furnishes essential services leading to accurate and uniform

physical and chemical measurement throughout the Nation's scientific community, industry,

and commerce; conducts materials research leading to improved methods of measurement,

standards, and data on the properties of materials needed by industry, commerce, educational

institutions, and Government; provides advisory and research services to other Government

agencies; develops, produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

Absolute Physical Quantities^ — Radiation Research — Chemical Physics —

Analytical Chemistry — Materials Science

THE NATIONAL ENGINEERING LABORATORY provides technology and technical ser- vices to the public and private sectors to address national needs and to solve national problems; conducts research in engineering and applied science in support of these efforts; builds and maintains competence in the necessary disciplines required to carry out this research and technical service; develops engineering data and measurement capabilities;

provides engineering measurement traceability services; develops test methods and proposes

engineering standards and code changes; develops and proposes new engineering practices;

and develops and improves mechanisms to transfer results of its research to the ultimate user. The Laboratory consists of the following centers:

Applied Mathematics — Electronics and Electrical Engineering^ — Manufacturing

Engineering — Building Technology — Fire Research — Chemical Engineering^

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts

research and provides scientific and technical services to aid Federal agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and

economy in Government operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing the Federal Information Processing Standards Program, developing Federal ADP standards

guidelines, and managing Federal participation in ADP voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal agencies; and

provides the technical foundation for computer-related policies of the Federal Government.

The Institute consists of the following centers:

Programming Science and Technology — Computer Systems Engineering.

'Headquarters and Laboratories at Gaithersburg, MD, unless otherwise noted; mailing address Washington, DC 20234. ^Some divisions within the center are located at Boulder, CO 80303.
Semiconductor Measurement Technology:

TXYZ: A Program for Semiconductor

IC Thermal Analysis

John Albers

Semiconductor Devices and Circuits Division

Center for Electronics and Electrical Engineering

National Engineering Laboratory

National Bureau of Standards

Washington, DC 20234

Z

Q

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

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Issued April 1984
Library of Congress Catalog Card Number: 84-601026

National Bureau of Standards Special Publication 400-76

Natl. Bur. Stand. (U.S.), Spec. Publ. 400-76, 65 pages (Apr. 1984)

CODEN: XNBSAV

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1984

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, DC 20402
Table of Contents

Abstract 1

I. Introduction 2

n. Solution of Steady-State Heat Flow Equation: Single Rectangular Layer . 4

ni. Solution of Steady-State Heat Flow Equation: Three-Layer Structure ... 8

IV. Form of the Function U(n,ni) 18

V. Behavior of Fourier Coefficients: Small Values of the Argument 21

VI. Behavior of Fourier Coefficients: Large Values of the Argument 24

VII. Special Case of Power Source Covering Top Surface 30

VIII. Effect of Upper Summation Limits on Temperature Profile 32

DC. General Discussion of TXYZ Program . 34

X. Selected Program Examples 36

XI. References 38

XII. Appendix A: TXYZ Program Listing 50

iii
List of Figures

Figure 1. This figure presents the geometry of the three-layer structure

in which the steady-state temperature is calculated 39

Figure 2. This figure shows the n-dependence of the product of the power

density function and the cos terms in eq (92) 40

Figure 3. The top-layer Fourier coeflBcient, ri(n, m, 0), is shown as a function of n and m and at 2r = 0 for a three-layer structure where all three

layers have different thermal conductivities 41

Figure 4. The top-layer Fourier coefficient, 7-1(71, m, 0), is shown as a function of n and m and at 2: = 0 for a three-lay^j structure where all three

layers have the same thermal conductivity 42

Figure 5. This figure shows the dependence of the three Fourier coefficients

on n at the top of each of the layers 43

Figure 6. The behavior of the top-layer Fourier coefficient, ri[n,m ==

0, z = 0), along the n axis is presented in this figure 44

Figure 7. This figure shows the calculated temperature as a function of

the number of terms used in the sum (eq (98)) for the top of each of the

three layers 45

Figure 8. This figure presents the calculated temperature along the midfine of a 200 mil by 200 mil rectangular structure with a single, uniform

power density 1 mil by 1 mil heat source located at the center of the

surface of the top layer 46

Figure 9. By way of contrast with the previous figure, the calculated temperature along the midline of a 200 mil by 200 mil three-layer structure

is presented in this figure 47

iv
Semiconductor Measurement Technology:

TXYZ: A Program For Semiconductor IC Thermal Analysis

John Albers

Semiconductor Devices and Circuits Division

National Bureau of Standards

Washington, D.C. 20234

Abstract

A computer program, TXYZ, for the thermal analysis of semiconductor integrated

circuits is presented and its applications are discussed. The program makes use

of the closed form, analytic solution of the steady-state heat flow problem for a

rectangular three-layer structure with multiple heat sources on the top layer. The

temperature may be obtained for any point or set of points in the structure and

is useful in the determination of the steady-state thermal response of IC chips and

packages.

Key words: FORTRAN; Fourier analysis; integrated circuit; semiconductor; steadystate heat flow; thermal analysis; thermal resistance.

1
EsTTRODUCTION

Since the introduction of semiconductor integrated circuits, one of the most important sources of device failure is the lack of temperature control. Hence, an accurate

physical picture of the temperature distribution in the device package under the

power condition of actual operation is of utmost importance. The purpose of this

report is to describe the program, TXYZ, which has been developed for the thermal

analysis of integrated circuit packages. In particular, the basic physical model and

the mathematical analysis are described and an annotated listing of the program,

along with sample data, is presented.

The physical and mathematical model used here is taken in part from work previously

carried out by Kokkas [1]. In order to provide the reader with a self-contained

document for using TXYZ, much of the development presented by Kokkas has

been worked out and specialized to the steady-state heat flow problem. The discussion which follows presents an annotated description of Kokkas' analysis with

additional material added where necessary. In particular, the equations have been

analyzed in detail so as to investigate the convergence of the solutions used in

the numerical implementation. Specifically, the form of the solutions for small

and large values of the argument is shown to require special consideration to avoid

numerical overflow problems.

This report is naturally broken up into two parts. The first deals with the mathematical and numerical details of the construction of the program. The second

deals with the program and its use. For those readers who are interested in the use

of the program, the section entitled "GENERAL DISCUSSION OF THE TXYZ

2
PROGRAM" begins the portion of the report where the program and specific examples are discussed.

UNITS

The American semiconductor industry has traditionally used mixed English and

metric units, but presently there is a trend in the direction of the International

System (SI) units. For the purpose of conversion, it should be noted that 1 mil =

0.001 in. = 25.4 fim and that 1 cal/s = 4.184 W.

3
SOLUTION OF STEADY-STATE HEAT FLOW EQUATION:

SINGLE RECTANGULAR LAYER

Consider a material of uniform thermal conductivity {ki), in the form of a rectangular box of lateral dimensions L^, Ly and thickness Li . The problem is to determine

the temperature, T{x,y,z), inside the material. It is assumed that the temperature

satisfies the steady-state heat flow equation [2]

V^T(x,y,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 there is no heat flow out of the lateral boundaries of the material,

i.e.,

dx x=o,L, dy j,=o,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 with only the lateral boundary conditions being 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 0,Lx and 0, Ly. This is defined as [3]

4
ifx, fy, z)= / T{x, y, z) exp(-2;ri(a;/:, + yfy))dxdy, (3)

Jo'0 Jo

where f^, fy are the Fourier transform variables which are conjugate to the variables

2, y. The inverse Fourier transform is defined as /•f• +00oo /•r-too +00

/

T{fx, fy, z) exp(2;ri(a;A + yfy))dhdfy.

-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

/+00 p +OO

f

N

J

r(/:c,/j,,^)exp(2;rty/j,)2;r/i<-sin(27r/:ez) + icos(2Kfxx)>dfxdfy, (5)

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

would require that the cos term be removed. Next, consider the resulting expression

at the other lateral boundary, i.e., at z = Z/^ where it is supposed to be zero. The

only way that this could be the case is if the argument of the sin function is an

integer times tt, or that the Fourier transform variable is of the form

= ^-

The same argument applies to the y-dependent portion. Then, in the Fourier representation, the temperature (with the above lateral boundary conditions) may be written

as

/+00 /> +oo

I T{n,'m,z)cos(n7rx/Lx)cos{in7ry/Ly)dfxdfy. (7)

-oo J —oo

5
The Fourier transform equation is now written as

r(n, m, z)= / / T(x,y,z)cos[nKx/Lx)cos(fmry/Ly)dxdy. (8)

As the system is of finite size, it is convenient to write the integral in eq (7) as a sum over the Fourier cos terms which fit into the rectangular geometry. Further,

as the cos function is symmetric around the origin, the sums may be written over

only the non-negative 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

df — A— — ^ ^ _ -'L — J_

(g)

L^r. Lx I^X ^r.

Then, the Fourier representation of the temperature may be written as

rpf \ 4r(n, m, z) cos{nwx/Lx) cos{mny/Ly) l[x,y,z}= 2^ 2^ ——— . , (lOj

where 6nn' is the Kronecker delta and is equal to unity if n = n' and zero otherwise.

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

and

52

^ cos{n7rx/Lx) = —{riTr/Lx) cos{n7rx/Lx), (12)

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

E°°

4cos(w7r2;/La;)cos(m;rt//L«) r , ,^ , ,^ ^0 \ ,

(13)
As the sum is zero for arbitrary values of the variables x and y (and the cos terms

are, in general, nonzero), then a necessary and sufficient condition that eq (13) is

satisfied is that

|-(n7r/L,)2 - [rm^lLyf + ^^T{n,m,z) = 0. (14)

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

the variable, 7, is defined as

then the above equation may be written in the form

^r(n, m, z) - i^T(n, m, z) = 0. (16)

The solution of this equation is

r(n, m,z) = a cosh(7>2:) + ^ sinh(7>2r), [17]

where the coefficients ot and ^, which may be functions of 7, are determined from the

two -ar-dependent boundary conditions. The above equation is the general solution for

the ^-dependent Fourier expansion coefficients for a single rectangular layer. In the

discussion of the problem of a rectangular y-layer structure where all of the layers

have the same lateral dimensions, the solution in each of the layers can be expressed

in the form of the above equation where the coefficients are to be determined from

the two 2r-dependent boundary conditions appropriate to each of the layers. This

will be used in the next section where the three-layer problem will be discussed.

7
SOLUTION OF STEADY-STATE HEAT FLOW EQUATION:

THREE-LAYER RECTANGULAR STRUCTURE

The basic problem considered in this section is the calculation of the three-dimensional

temperature distribution in a three-layer structure which is assumed to be of the

geometric form presented in figure 1. The three layers are characterized in terms

of the thermal conductivities and thicknesses /c,-, L,(t = 1,2,3). One particular

example of a three-layer structure is an IC package. For this particular case, the

top layer is the semiconductor device, whereas the middle and bottom layers are

the die attach and substrate layers. Clearly, the thicknesses and thermal conductivities of these three layers are of considerable importance in the dissipation of heat

generated by the power sources on the surface of the semiconductor device. 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.

All three layers of the model are assumed to have the same lateral dimensions.

This assumption greatly simplifies the analysis in terms of Fourier series which are common to all three layers. In addition, it is assumed that there is no heat flow

out of the lateral boundaries of the structure due to either convection or radiation.

Again, this makes the analysis more manageable.

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

assumptions:

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

8
thicknesses may be different,

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

and heat flow in the structure takes place by conduction,

3) there are no heat losses due to interconnections to the chip (the top layer),

4) there is no thermal contact resistance between the various layers which are in

contact,

5) the heat sink, which is in contact with the bottom layer, is ideal and has a

temperature equal to ambient,

6) each layer is of uniform, isotropic, temperature-independent thermal conductivity, and

7) there is no input power density inside the structure, heat is generated only on

the top surface.

The starting equation for the analysis is the steady-state homogeneous heat flow

equation as discussed for the single-layer problem,

V'T(x,y,z) = 0. (18)

In the following discussion, the temperature in the first layer will be referred to as

Ti(x, y, z), while the temperature in the second and third layers will be denoted by

T2(x,y,z) and T2,(x,y,z), respectively. Built into the above assumptions is the fact

that the temperature and its normal derivative (proportional to the normal heat

flow) are continuous at the interfaces between the layers. The assumptions that heat

9
enters the structure only on the top surface through the heating elements and that

the bottom of the bottom layer is at the same temperature as the heat sink are also

available. These provide the six ^r-dependent boundary conditions on the system of

equations generated by using the solution of the homogeneous single-layer problem.

These boundary conditions may be expressed as follows. First, the assumption that

there is no heat flow out of the lateral boundaries is written as

dx x==o,L, dy y=o,Ly

Next, the assumption that heat enters only where the power is applied in the first

layer is expressed as

OZ z=0

where ki is the thermal conductivity of the top layer. Further, the assumption

that the temperature is continuous across the interfaces between the layers may be

written as

Tt{x,y,z)U^-L, = T2(x,y,z)U=.-L,, (21)

T2{x,y,z)\^—(Li+La) = Ts{x,y,z)\^—(Lx+l,)- (22)

The assumption that the heat flow is continuous across the interfaces between the

layers is expressed by the conditions that

dTi{x,y,z) _ dT2(x,y,z)

f^i ^ — » , {^^} OZ z Li OZ z=-Li

dT2{x,y,z) _ dTi(x,y,z)

K2 — K3 • dz z (L1+L2) dz z (Li+La)

10
Finally, the assumption that the temperature is continuous across the interface

between the third layer and the heat sink is written as

0 (25)

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

In the above equations, k,- is the thermal conductivity of the i-th layer and Lg =

Li + L2 + ^3 • It is also 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. For the

steady-state situation, the power density may be written as

where the function U[x, y) is the weighting function which describes the geometry

and uniformity (or nonuniformity) of the heat-generating components, and Pq is

the steady-state power density per unit area. As discussed previously, the ambient

temperature trivially satisfies the heat flow equation. Hence, it is convenient to

consider all temperatures relative to ambient. Therefore, when the power density is

set equal to zero, the temperature will be zero as it is relative to ambient. As the

three layers are of the same lateral dimensions and there is no heat flow out of the

lateral boundaries, the temperature in each of the three layers may be written in

the form of eq (10). Further, the ^-dependent Fourier expansion coefficients will be

of the form of eq (17) and may be written as

(26)

ri(n, m, z) = ai cosh(7xr) + sinh(72?), (27)

T2{n, m, z) — a2 cosh(7>2r) + ^2 sinh(7z). (28)

11
Ts{n, m, z) — 0:3 cosh(72r) + sinli(72r). (29)

There are six unknowns involved in these solutions and they are specifically the

set of expansion coefficients a{,^i{i = 1,2,3). These coefficients may be explicitly

evaluated by means of the boundary conditions on the temperature and its derivative evaluated at the interfaces. The Fourier coefficients may be used directly in

the above boundary conditions as the temperature is a sum over the Fourier expansion coefficients. By substituting the above equations into the appropriate boundary condition equations, it is possible to obtain a system of six equations in six

unknowns. For the present case, this system reduces to one with four equations

in four unknowns. The standard method for solving this system is to use Cramer's

rule [4] with the Laplace method [4] for the evaluation of the determinants involved.

However, instead of using this method explicitly here, the Fourier coefficients for

each of the layers will be presented and will be shown to satisfy the heat flow equation and the appropriate boundary conditions. In particular, the Fourier coefficients

in the three layers are given by eqs (15-17) in Kokkas' paper. Specializing these to

the steady-state situation, these are

72 (n, m, z)=a\d cosh(7(Li -i- L2 + z))

sinh(7(Li -h L2 z)) [, (31)

(n, m, z) = AIB cosh(7(Li -\- z)) -\- C sinh(7(Li + z)) (30)

(32)

where

(33)

12
B = D cosh(7Z/2) + E smh(7L2), (34)

C = sinh(7L2) 4- £;cosh(7L2)|, (35)

D = siiih(7L3), (36)

E = coshir^Ls), (37)

K>2

U{n,m)= I I U{x,y)cos{n7rx/Lx)cos[m7ry/Ly)dxdy, (38)

^0 Jo

and

0

is the double Fourier cos transform of the power density uniformity function. Now,

it will be shown that the above are the solutions of the equation (see eq (16))

Ti(n, m, z) - 7^ n(n, m, z) =

0, (39)

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

the case as

and

Q2

cosh(7(L + z)) = 7^ cosh(7(L + z)), (40)

q2

sinh(7(L + z]) = 7^ sinh(7(L + z)), (41)

dz^

where L is a constant and is equal to Li, Li + L2, or Lg in eqs (30-32). Hence, eqs

(30-32) satisfy the 2r-dependent differential equation. Next, it will be shown that

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

13
these is that

«i ' ^ = P(n, m) = U(n, m)Po . 42

az z=o

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

as

^7j_(n^_r7V^I _ ^^^^] Q _|_ (7 cosh(7Zyi)|

oz z=0 { J

/ci7 i^i?sinh(7Li) + C/cosh(7Li) J J = U{n,m)Po. (43)

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

the bottom layer boundary condition, i.e.,

Tz[n,m,z%— = 0. (44)

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

Tz(n, m, z%=-L, = Asinh(7(L^ - Lg)) = 0. (45)

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

Ti(n,m,z%—Lr= T2(n,m,z%=.-Li, (46)

T2(n, m, z)\^—(L1+L2) = Ts{n, m, z)\^—(L1+L2), (47)

14
dTi{n, m, z) dr2(n,m,z)

dz L=-L.

=

dz U-L/

dT2{n,m, 3r) drs(n,m, z)

^ dz z=-{Li+L2) dz z^-{Li+L2)

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

Ti{n,m,z)\^— = Tlj^ cosh(7(Li - Li)) + Csinh(7(Li -

= AB = A|z)cosh(7L2) + J^sinh(7L2)|- (50)

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

T2{n,m,z)\^—Li = cosh(7(Li + L2 - ^i)) + ^sinh(7(Li H-Z^ -Li))|

= aI^D cosh(7L2) + E sinh(7L2)|. (51)

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

condition. Next, consider the second temperature continuity equation.

r2(n, m, z)\^—(L1+L2) = rs{n, m, z)\^—(L1+L2), (52)

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

T2(n,m,z)\^ (L1+L2) = A|Dcosh(7(Li+L2-Li-L2))+£;sinh(7(Li+L2-^i-^))|

= AD = Asinh(7Z/3).

15

(53)
The right-hand side is

rs{n, m, z%—{L^^L^) = Asinh(7(L^ - Li - L2)) == Asinh(7L3). (54)

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

to be verified is the first heat flow continuity equation, i.e.,

dTi{n,m, z) dT2(n,m,z)

dz 1 r dz 1 t' ^^^^ OZ z=—Li OZ z=—Li

Evaluating the left-hand side leads to

= = AC27^|^sinh(7L2) + £;cosh(7L2)|- (56)

The right-hand side may be evaluated as

^^ iji) sinh(7L2) + E cosh(7L2)|. dT2(nm,z)

^ = /C27A^i)sinh(7L2) + ^cosh(7L2)^ (57)

uZ z=—Li

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

second and third layers, i.e.,

dz z (L1+L2) dz z=-{Li-j-L2)

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

dT2(n,m,z)

K,2

dz z=-{Li+L2)

= /c27a|d sinh(7(Li -\- L2 - Li - L2)) E cosh(7(Li + L2 - Li - L2))|

16
(59)

The right-hand side is

dn{n,m,z = Ki^Acosh['^(Lz — Li — L2) = /C37j4.cosh(7Z/3). (60) dz z=-{Li+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 U(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 7 very small and 7 very large may give rise

to overflow or underflow problems when the program is constructed.

17
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

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, t/,) and lengths along the x- and

y-directions given by (lxi,lyi). 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, ?7(a;, 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

^0 Jo

(62)

18
The integrals can be simply evaluated to give the result that

Ui{n, m) =

(n7r)(m;r)

X (63)

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

to give the cumulative heat source effect. In addition, as indicated previously, if

there are any nonuniform heat sources, their effect can be constructed by means of

overlapping uniform heat sources. 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 cos transform over the range of n, m required by the sum in eq (10). 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

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

that the function U(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

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 sin function for small values

(64)

19
of the argument. In particular

sin(a;) = z - — + •

• (65)

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

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 (63) shows that Ui{n,m) 0

Uke 1/nm as m,n -> oo. In addition, the cos terms in eq (10), i.e., cos[nKx/Lx)

cos{m7ry/ Ly), do not approach a definite limit for large values of the argument.

Consequently, the product U{n^7n)cos(nKx/Lx)cos{mwy/Ly) tends to zero slowly.

An example of the behavior of the z-dependent portion of this product function,

i.e., U{n) cos{n7rx /Lx) is presented in figure 2. The heat source used for this figure

is a 1 mil by 1 mil heat source centered on the surface of a 200 mil by 200 mil

structure. The product function is evaluated at the midpoint of the heat source.

Because of the symmetry used, the function shows damped oscillatory behavior. It

is clear from the figure that the product function does not fall off to zero sufficiently

fast and hence does not provide for rapid convergence of the calculated temperature

using the Fourier representation.

lim U{n) = lim —

re—0 re—0 nw

(66)

20
BEHAVIOR OF FOURIER COEFFICIENTS:

SMALL VALUES OF THE ARGUMENT

As lias 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 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 so as

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 (30-37) 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,m)Po/ k.i 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

(30-33)). Then, for small 7,

Ef^'^, (67)

K2

21
Df^^Ls, (68)

(69)

Bf^qLs-h —1L2, (70)

and, remembering that the factor U{n,m)Po/i<ii has been included explicitly in eq

(10),

A^-'^. (71)

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,

ri(n,m,z)^ - —Ll^ + —1L2 + —(7%^3 + —\i(Li-\-z)\

7 /C3 [ K2 K,i( K2} J —< L3 + —L2 + —(m + ^) r

t K2 J Then,

lim Ti(n, m, z) = (Li + 2r) + —L2 + —L3. (72)

7-*0 /C2 /C3

Next,

T2(n,m,z) i —I7L3 + —7(^1 +1^ +^^)1,

7 «3 I K2 J

or

lim r2(n, m, ^r) = —L3 + ~{Li + L2 + -sr). (73)

7-*0 ^3 K2

22
And finally,

Ts{n,m,z) - —\-f(Li + L2 + + z)\,

or

lim Ts{n, m,z) = —(Li + L2 + -^'3 + z). (74)

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.

23
BEHAVIOR OF FOURIER COEFFICIENTS:

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(a:), 8{x) = sinh(a;),

and t{x) = tanh(a;). Making use of this shorthand notation, the Fourier coefficients

may be written as

(75)

(n, m, z) = A\ Dc{^(Li + L2 + z)) + Es^^iL^^ + L2 + -sr)) [ (76)

Ts{n,m,z) = A8(q(L^ + 0). (77)

where.

(78)

B = I>c(7L2) + Es(^L2), (79)

24
C ^ — D8(^L2) + Ec{^L2) , (80)

(81)

(82)

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

U[n, in)Po/Ki has been deleted from eq (78) for convenience. This factor may simply

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

to all three layers. Now, the above equations (eqs (75-82)) will be rewritten by

making use of the definition of the hyperbolic tanh, i.e., t(x) = 3{x)/c(x). First,

consider the coefficient C.

C = —

c

(83)

Next, the coefficient B may be written as

B = Dc{'^L2) + ^5(7^2)

B = 8{r^Ls)c(^L2) + ^c(7L3)«(7L2)

(84)

25
Then,

B5(7Li) + Cc(7Li) =

= c(^LMnL2)c[iU){t[^L^)t[^^^^ (85)

^

AC2 K.I Kl )

Then, the coefficient A may be written as

(86)

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

n(7) = -, (87)

{«(7i,)<(7Li) + =f t(iZ„)%I^) + "^ti^LM^L,) + ^1

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

(88)

7c(7L3)c(7L2)c(7A)'

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

the Fourier coefficients as given by eqs (75-77) may be written as

Ti (n, m, z) =

n(7)<7(^i + z))

72 (n, m,z) =

and

73 (n, m, z) =

7c(7Li)c(7L2)c(7L3)

(91)

In eqs (89-91), 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

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

term PoU{n,m)/Ki has been written out explicitly and is no longer contained in

00 00 4U{n, m)r,(n, m, z) cos{nwx/Lx) cos{m7ry/ Ly)

i^nO + l)((^mO + l)L:^LyK,i

Ti{x,y,z) = Po J2 E (92)

«==0 m=0

27
the Fourier coefficients. It is important in obtaining the solution in a;, 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.

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.

In the subsequent discussion, 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 cos terms in eq (92) 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. Consequently, some discussion of the n and m behavior

of the Fourier coefficients is warranted. To show how the Fourier coefficients behave

as a function of n and m (or 7), these functions were studied in some detail. In figure

3, a three-dimensional plot of the n and m dependence of the top layer Fourier

coefficient is shown at z = 0 for a three-layer structure. In the figure, it can be seen

that the Fourier coefficients are peaked at the origin of the n, m plane and approach

zero as n and m become large. In addition, figure 4 contains the same type of plot

for the situation where all three layers have the same thermal conductivity, i.e., for

a thick one-layer structure.

28
The dependence on n of the three Fourier coefficients at the top of each layer is

contained in figure 5 for the structure used to generate figure 3. The curve marked

by A shows the behavior on the top surface of the structure, i. e., at z = 0. The

curves denoted by B and C present the results for the top surfaces of the second

and third layers, respectively. Figure 6 is a log-log plot form of figure 5a and is

especially important as it shows the slow convergence of the top surface Fourier

coefficient which will be discussed in detail in the section on the effects of the upper

limit on the calculated temperature.

29
SPECIAL CASE OF POWER SOURCE COVERING TOP SURFACE

As a special case of eq (92), 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 will be 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 (63),

Ui (n, m) =

(93)

For the particular case of uniform surface coverage, xi = yi

=

0, Ixi = Lx and

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

Ui{n, m) = sin(m7r)|. (94)

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 sin function gives rise to the result that

Ui(n,m) = LxLy6nQ6mQ- (95)

If this form of the U[n, m) function is substituted into the equation for the temperature in each of the three layers (eq (92)) and the form of the Fourier expansion

coefficients as 7 — 0 (eqs (72-74)) is used, it is readily shown that the temperatures

30
in the three layers may be written as

r.(.,!,,^) = Po{^ + ^ + ^} (95)

Ki K2 K,^ )

Ux,y,.)^Po{ '^+^ +

' + ^] (96)

Ts(x,y,z) = Fo| y (97)

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.

31
EFFECT OF UPPER SUMMATION LIMITS ON TEMPERATURE

Previously, the dependence of the Fourier coefficients on n and m has been discussed.

Figures 3 to 6 contain typical results of this behavior. The purpose of the present

section is to show how this behavior is mirrored in the calculated temperature. In

particular, eqs (89-91) will be used in eq (92) with a variable upper summation Hmit

in the latter equation. To simplify the analysis, a single heat source with y,- = 0

and Ij/i = Ly will be used. The width of this stripe heat source will be Ixi — 1

mil. Finally, the lateral dimensions of the three-layer structure will be taken as

Lx = Ly = 200 mil. This particular choice will require only the m = 0 term in

the sum to be retained while needing a large value of n terms in order to obtain

convergence of the sum. As there is complete coverage along the y-direction, the

7n-dependent portion of eq (95) may be used in eq (92) to obtain

sidered in this section. Curves A, B, and C in figure 7 show the temperature calculated at the center of the heat source as a function of the number of terms in

the sum, N, for the top of the first, second, and third layers, respectively. It is

clear from the curve in figure 7a that the calculation of the surface temperature (at

z = Qi) may require up to at least 350 terms in the sum while, from figures 7b and

7c, the temperature below the surface may need only 20 terms to be retained in the

sum. From the curve depicted in figure 7a, it may be argued that the number of

terms needed to adequately represent a feature size of Ax (typical of the lateral size

of a heat source) should be on the order of Lx/ Ax. For the present case, this ratio

(98)

where the dependence of T^(x, y, z) on N, the upper limit of summation, is con32
is about 200. Deeper into the structure, the heat flow has caused the approximate

size of the heat current to spread out and hence fewer terms are necessary. This is

borne out by the behavior of the temperature in the figure.

33
GENERAL DISCUSSION OF THE TXYZ PROGRAM

The annotated listing of the program is contained in the appendix of this report. In

the present section, several aspects of the program will be discussed. The purpose

of this discussion is to present the user with information concerning the implementation and/or modification of the program.

The first item is the format in which input data are to be entered. As listed

in the program, data are to be entered in a fixed-field format. However, many

FORTRAN compilers support free-field or list-directed read statements. Free-field

read statements are convenient as they allow the user to change input data without

having to be concerned with the tedious process of lining up of the data required

by fixed-field read statements. Hence, if the user's FORTRAN compiler allows for

free-field read statements, it is strongly suggested that the program be edited such

that all read statements are of this form.

The next item is concerned with the situation of the hyperbolic functions being

or not being built-in functions. If these functions are contained in the computer's

relocatable (or object code) mathematics library, then there are no changes required

in the code. If, however, these are not, then it is necessary to write separate function

subroutines for the hyperbolic functions. As the hyperbolic functions are easily

generated from the exponential function, this is a straightforward process.

Next, it is important to discuss the number of heat sources and the dimension

statements used in the program. In its present form, TXYZ will allow up to twenty

heat sources. This was deemed sufficient for most problems of interest. However,

should the user want to use more than this number of heat sources (especially

34
for the case of a number of nonuniform heat sources), it is necessary to change

the corresponding dimension statements. This is commented upon in the program

and the change is relatively straightforward. The only other comment concerning

dimension statements has to do with the number of points at which the temperature

is to be calculated and the maximum number of terms to be used in the calculation

of the temperature using eq (92). The program allows up to 501 as the maximum

number for both of these quantities. This was considered to be adequate for the

resolution of the spatial variation of the temperature and for the calculation of the

temperature (even for the smallest physically realizable heat source on a structure of

typical lateral dimensions). Hence, it is felt that 501 is a good maximum value which

gives sufficient detail without excessive use of computer memory space. Should

memory allocation be a problem, it may become necessary to reduce these numbers

in the dimension statements. It is important to carefully evaluate any reduction in

the maximum number of terms in the Fourier representation of the temperature

in light of the slow convergence of the surface temperature for small heat sources.

Finally, in regard to changing the dimension statements, it is important to change

the dimension statements in not only the main program but also in the function

subroutines as the dimensioned variables are declared in common and hence must

have the same dimensions everywhere (in the main program and the function

subroutines) to avoid any register difficulties.

A copy of this FORTRAN program may be obtained from the author by sending a

letter of request and a computer tape.

35
SELECTED PROGRAM EXAMPLES

In order to assist the user in the implementation of the TXYZ program, several

examples are presented in this section. In particular, the input data files are listed

and plots of the temperature profiles are presented in the accompanying figures.

The input data files have been annotated in order to facilitate familiarity with the

reading of these files. This is especially important because of the number of input

variables which set various "switches" inside the program. The corresponding data

files will be sent along with the FORTRAN listing of the TXYZ program so as to

simplify the implementation of the program and its subsequent use.

Two examples will be discussed in this section. The first is for the situation of a

single small heat source on the surface of the top layer. The second example is that

of a number of heat sources on the top layer. In particular, one of the heat sources

is nonuniform and its construction from uniform heat sources is illustrated in the

ccorresponding input data file.

The first example is that of a 1 mil by 1 mil uniform heat source on the surface

of a 150 mil by 200 mil three-layer structure. The specific structure is that of 15

mils of sihcon on 2 mils of die attach material on 30 mils of a substrate. This threelayer structure will be used in both this example and the one to follow. The input

data file for this first example is listed, with annotation, in Table 1. The calculated

temperature is presented in figure 8.

The second example is for the same three-layer structure as used in the first example.

However, in this case, there are several heat sources, with one of them being

nonuniform and built up from uniform heat sources. In particular, there are three

36
square uniform heat sources on the right portion of the surface and one uniform

heat source in the middle of the surface. The nonuniform heat source is a single

rectangular area. The surface temperature calculated along the line 0 < x < 200,

y = 101 is presented in figure 9. The input data file is presented in Table 2.

In addition to these examples, the reader's attention is drawn to the recent review

[5] where this program is used in the thermal evaluation of VLSI packages using

test chips.

ACKNOWLEDGMENTS

The author would like to thank Frank F. Oettinger for his continued interest during

the course of this work. His relentless use of the program brought many of the early

bugs to the surface.

The author would also like to thank Stephen E. Ross who changed the code into a

portable FORTRAN form. His assistance in getting the three-dimensional plots is

also appreciated.

37
REFERENCES

1. Kokkas, A. G., Thermal Analysis of Multiple-Layer Structures, IEEE Trans.

Electron Devices ED-21, 674-681 (1974).

2. Carslaw, H. S., and Jaeger, J. C, Conduction of Heat in Solids, Second Edition

(Oxford University Press, 1959).

3. Churchill, R. V., Fourier Series and Boundary Value ProblemSj Second Edition

(McGraw-Hill, 1963).

4. Ayres, F., Jr., Matrices (Schaum Publishing Co., 1950).

5. Oettinger, F. F., Thermal Evaluation of VLSI Packages Using Test Chips - A

Critical Review, Solid State Technology 27 (2), 169-179 (1984).

38
I

DISTRIBUTED

HEATERS

HEAT SINK

\ Ls SEMICONDUCTOR

W/ ^ 1 Lc t DIE ATTACH

Lj t SUBSTRATE

Figure 1. This figure presents the geometry of the three-layer structure in which

the steady-state temperature is calculated.

39
-0.006 I 1—1 I \ I I I I I I I \ I L__j I I ! L_J

0 100 200 300 400 500 600 700 800 900 1000

n

Figure 2. This figure shows the n-dependence of the product of the power density

function and the cos terms in eq (92), i.e., U{n) cos{n7rx/Lx), evaluated at the center

of a single small heat source which is in the form of a thin stripe along the ydirection. This specific form of the product is used as there is complete coverage

along the y-direction which means that only the m = 0 term contributes. The

general feature seen in this figure is that this function does not asymptote to zero

for n on the order of 500. This general kind of behavior has been found to be the

case for all situations considered.

40
ri (n.m.O)

Figure 3. The top-layer Fourier coefl&cient, ri(n, m, 0), is shown as a function of n

and m and at = 0 for a three-layer structure where all three layers have different

thermal conductivities. The specific structure is for 15 mils of silicon (/ci = 0.00267)

over 2 mils of die attach (/C2 = 0.00064) over 30 mils of substrate material (/C3 =

0.00999). The peak value of ri is at the origin of the n,m plane.

41
r 1 (n.m.O)

Figure 4. The top-layer Fourier coefficient, Ti{n,m,0), is shown as a function of n

and m and at 2: = 0 for a three-layer structure where all three layers have the same

thermal conductivity, i.e., a thick one-layer structure. The specific example is for 47

mils of silicon (/c = 0.00267). As with the previous figure, ri is peaked at the origin

of the n, m plane. Both this and the previous figure are meant to convey qualitative

features of the Fourier coefficient. More detailed information will be presented in

the next two figures where specific attention will be focused upon just the n axis

(with m = 0) behavior.

42
360 f

'—'—I

> ' I —'

I I 1 1

1

1

1—-I—I

I I 1 1

340 - -

10° 10^ 10^ 10^

n + 1

Figure 5. This jBgure shows the dependence of the three Fourier coefficients on n at

the top of each of the layers. The curve denoted by A is Ti{n,m = 0,z = 0) (for

the top of the first layer). B shows the Fourier coefficient, T2{n,m = 0,z = —Li)

(for the top of the second layer). Finally, the curve denoted by C represents the

Fourier coefficient, Ts{n,m = 0,z = —Li — L2), i.e., at the top of the third layer.

This figure and the next represent a cross section along the n axis of figures of the

type presented in figures 3 and 4.

43
10° 10^ 10^ 10^

n+ 1

Figure 6. The behavior of the top layer Fourier coefficient, Ti{n,m = 0,z =

0),

along the n axis is presented in this figure. These results (curve A of fig. 5) are presented in a log-log plot to show the slow convergence z = 0. When curves B

and C of figure 5 are plotted in the same manner, they fall off to zero much more

rapidly.

44
12 - B (i = 2, z = -Li)

10

8

C (i = 3, z = -Li-L2)

4 ' ' '

I 'I I I 'I I ' 'I I I I I ' ' I ' ' ' I 'I I ' ' ' I ' ' ' 0 50 100 150 200 250 300 350 400 450 500

N

Figure 7. This figure shows the calculated temperature as a function of the number

of terms used in the sum (eq (98)) for the top of each of the three layers. Curve A is Ti^(z = 0) while curves B and C represent T^{z = -Li) and T^{z = -Li - L2),

respectively. This figure clearly shows how the behavior of the Fourier coefficients

for each of the layers is mirrored in the corresponding calculated temperatures.

45
20r—r—

Oo

LUGC

Z)

\-

<

OC

LU

0.

LU

1—I—I—I—I

nr\—I—I—1—I—I—I—I—I—I—I—

I I r 92 94 96 98 100 102 104

DISTANCE (mils)

106 108 110 112

Figure 8. This figure presents the calculated temperature along the midline of a 200

mil by 200 mil rectangular structure with a single, uniform power density 1 mil by

1 mil heat source located at the center of the surface of the top layer. Note that

only the temperature in the region of the heat source is plotted, i.e., from 90 mils

to 110 mils. Also, note that the temperature falls off rapidly away from the heat

source. The specific structure is that of 15 mils of silicon on 2 mils of die attach

material on 30 mils of substrate material. The annotated input data file used to

obtain this temperature data is to be found in Table I.

46
20 40 60 80 100 120 140 160 180 200

DISTANCE (mils)

Figure 9. By way of contrast with the previous figure, the calculated temperature

along the midline of a 200 mil by 200 mil three-layer structure is presented in this

figure. However, in this case, two nonuniform heat sources and one uniform heat

source are used. The three-layer structure is the same as used in the previous figure.

The construction of the nonuniform heat sources from a number of uniform heat

sources is contained in the annotated input data file contained in Table 11.

47
TABLE I INPUT DATA FOR SINGLE HEAT SOURCE

INPUT DATA DESCRIPTION OF INPUT DATA

150

0.00267

0.00064

0.00999

500

90 0.1

200

15

230

500

1201

0

75.5

0

0

1 0.093423

99 1 75

0

0

0 0

0

X AND Y DIMENSIONS OF RECTANGULAR STRUCTURE

THICKNESS AND THERMAL CONDUCTIVITY OF TOP LAYER

THICKNESS AND THERMAL CONDUCTIVITY OF MIDDLE LAYER

THICKNESS AND THERMAL CONDUCTIVITY OF BOTTOM LAYER

UPPER SUMMATION LIMITS FOR N AND M SUMMATIONS

lEDGEX (=0 FOR SINGLE POINT, =1 FOR NUMBER OF POINTS)

NUMBER OF POINTS ALONG X, FIRST POINT, STEP INCREMENT

lEDGEY (=0 FOR SINGLE POINT, =1 FOR NUMBER OF POINTS)

Y POINT FOR CALCULATION

lEDGEZ (=0 FOR SINGLE POINT, =1 FOR NUMBER OF POINTS)

Z POINT FOR CALCULATION

NUMBER OF HEAT SOURCES AND POWER DENSITY

X COORDINATE, LENGTH ALONG X AXIS, Y COORDINATE, LENGTH

ALONG Y AXIS

48
TABLE II

INPUT DATA FOR MULTIPLE, NONUNIFORM HEAT SOURCES

INPUT DATA DESCRIPTION OF INPUT DATA

200 200

15 0.00267

2 0.00064

30 0.00999

500 500

1200 0 1 0

101

0

0

7 0.001

10 65 10 180

30

50

100

160

45

25

2

10 180

10 180

100 2

20 10 20

160 20 90 20

160 20 170 20

0

0

0 0

0

X AND Y DIMENSIONS OF RECTANGULAR STRUCTURE

THICKNESS AND THERMAL CONDUCTIVITY OF TOP LAYER

THICKNESS AND THERMAL CONDUCTIVITY OF MIDDLE LAYER

THICKNESS AND THERMAL CONDUCTIVITY OF BOTTOM LAYER

UPPER SUMMATION LIMITS FOR N AND M SUMMATIONS

lEDGEX (=0 FOR SINGLE POINT, =1 FOR NUMBER OF POINTS)

NUMBER OF POINTS ALONG X, FIRST POINT, STEP INCREMENT

lEDGEY (=0 FOR SINGLE POINT, =1 FOR NUMBER OF POINTS)

Y POINT FOR CALCULATION

lEDGEZ (=0 FOR SINGLE POINT, =1 FOR NUMBER OF POINTS)

Z POINT FOR CALCULATION

NUMBER OF HEAT SOURCES AND POWER DENSITY

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER #1

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER #2

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER #3

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER #4

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER US

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER #6

X COORDINATE, LENGTH ALONG X AXIS,

ALONG Y AXIS FOR HEATER #7

COORDINATE, LENGTH

COORDINATE, LENGTH

COORDINATE, LENGTH

COORDINATE, LENGTH

COORDINATE, LENGTH

COORDINATE, LENGTH

COORDINATE, LENGTH

49
APPENDIX A~TXYZ PROGRAM LISTING

MAIN PROGRAM LISTING

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

C INSTRUCTIONS FOR INPUT DATA FILE SETUP

C DATA TO BE READ BY THIS PROGRAM ARE SET UP AS FOLLOWS:

C

C RLX (X DIMENSION OF 3 LAYER STRUCTURE)

C RLY (Y DIMENSION OF 3 LAYER STRUCTURE)

C RLS (THICKNESS OF TOP LAYER)

C RKS (THERMAL CONDUCTIVITY OF TOP LAYER)

C RLC (THICKNESS OF MIDDLE LAYER)

C RKC (THERMAL CONDUCTIVITY OF MIDDLE LAYER)

C RLI (THICKNESS OF BOTTOM LAYER)

C RKI (THERMAL CONDUCTIVITY OF BOTTOM LAYER)

C NUP (UPPER LIMIT OF N SUM, X DIRECTION)

C MUP (UPPER LIMIT OF M SUM, Y DIRECTION)

C

C lEDGEX (=0 FOR DEFAULT, =1 TO ALTER RANGE OF X VALUES)

C IF IEDGEX=1 THEN THE FOLLOWING SEGMENT OF INPUT IS READ

C ILX (THE NUMBER OF INCREMENTS IN X TO BE USED)

C (ILX+1=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 ELSE IEDGEX=0 READ XI (THE CONSTANT X VALUE)

C

C lEFGEY (=0 FOR DEFAULT, =1 TO ALTER RANGE OF Y VALUES)

C IF IEDGEY=1 THEN THE FOLLOWING SEGMENT OF INPUT IS READ

C ILY (THE NUMBER OF INCREMENTS IN Y TO BE USED)

C (ILY+1=THE NUMBER OF Y VALUES TO BE USED)

C Y1 (THE VALUE OF THE FIRST POINT IN Y)

C STEPY (THE INCREMENT IN Y)

C ELSE IEDGEX=0 READ XI (THE CONSTANT X VALUE)

C

C lEDGEZ (=0 FOR DEFAULT, =1 TO ALTER RANGE OF Z VALUES)

C IF IEDGEZ=1 THEN THE FOLLOWING SEGMENT OF INPUT IS READ

C ILZ (THEN NUMBER OF INCREMENTS IN Z TO BE USED)

C (ILZ+1=THE NUMBER OF Z VALUES TO BE USED)

C Z1 (THE VALUE OF THE FIRST POINT IN Z)

C STEPZ (THE INCREMENT IN Z)

C ELSE IEDGEX=0 READ XI (THE CONSTANT X VALUE)

C

C NOTE: AS THE CALCULATION TAKES THE Z VARIABLE TO BE ZERO OR

C NEGATIVE, ENTER Z1 AND STEPZ AS POSITIVE QUANTITIES.

C THE PROGRAM CHANGES THE SIGN OF Z1 AND STEPZ TO MAKE

C THE CALCULATION FOR NEGATIVE Z.

CC NSOUR (NUMBER OF HEAT SOURCES, UP TO 20)

C PO (POWER DENSITY)

C XSOUR(I) (X COORDINATE OF ORIGIN OF 1ST SOURCE)

50
APPENDIX A~TXYZ PROGRAM LISTING

C YSOUR(I) (Y COORDINATE OF ORIGIN OF 1ST SOURCE)

C RLXSOUR(I) (LENGTH ALONG X DIRECTION OF 1ST SOURCE)

C RLYSOUR(I) (LENGTH ALONG Y DIRECTION OF 1ST SOURCE)

C REMAINING HEAT SOURCES WITH SAME INPUT STRUCTURE AS

C

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

c DIMENSION X(501),Y(501),Z(501),COSYT(501)

DIMENSION ARUZER(501,501 ),ARFUNZ(501,501

)

DIMENSION XSOUR(20),YSOUR(20),RLXSOR(20),RLYSOR(20)

COMMON RKS,RKC,RKI,RLX,RLY,RLS,RLC,RLI

COMMON NSOUR,XSOUR,YSOUR,RLXSOR,RLYSOR

C

c *****************************************************

c

C DOCUMENTATION AND BACKGROUND

C THIS PROGRAM CALCULATES THE SURFACE TEMPERATURE T(X,Y,Z) DUE TO DC

C POWER INPUTS ONLY. THE TEMPERATURE IS THE TEMPERATURE RELATIVE

C TO THE AMBIENT. THE SPECIFIC EQUATIONS USED ARE GIVEN IN

C EQUATIONS (13)-(23), WITH S=0 (STEADY-STATE CONDITION),

C IN THE PAPER BY KOKKAS (REF: "THERMAL ANALYSIS OF MULTIPLEC LAYERED STRUCTURES" BY ACHILLES G. KOKKAS, IEEE TRANS. ELEC. DEV.

C VOL. ED-21, NO. 11, 674-681 (1974)).

C VARIABLES USED: THE VARIABLES LISTED AS REAL IN THE ABOVE

C ARE THE FOLLOWING-KS,KC, AND KI ARE THE THERMAL CONDUCTIVITIES

C OF THE SEMICONDUCTOR, CONDUCTOR, AND INSULATOR, RESPECTIVELY.

C LX AND LY ARE THE LATERAL DIMENSION OF THE CHIP WHILE

C LS,LC, AND LI ARE THE THICKNESSES OF THE SEMICONDUCTOR, THE

C CONDUCTOR, AND THE INSULATOR, RESPECTIVELY.

C IMPORTANT NOTE: WHILE THE VARIABLES HAVE THE NOTATION WHICH

C SEEMS TO IMPLY A SEMICONDUCTOR, A CONDUCTOR, AND AN INSULATOR,

C THESE REFER TO THE WAY IN WHICH KOKKAS FORMULATED THE PROBLEM.

C THE THING TO KEEP IN MIND IS THAT THE TOP LAYER HAS (LS,KS),

C THE MIDDLE LAYER HAS (LC,KC), AND THE BOTTOM LAYER HAS (LI,KI) . C IT IS NOT NECESSARY THAT THEY BE WHAT THEY SEEM TO BE, THEY ARE

C DETERMINED BY THE RESPECTIVE THICKNESSES, L, AND CONDUCTIVITIES, K

C

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

c

1 FORMAT(IHI)

2 FORMATdX, 'STEADY-STATE THERMAL ANALYSIS CALCULATION USING EQS

1 (13)-(23) OF KOKKAS'/)

3 FORMATdX, 'THERMAL CONDUCTIVITIES AND LAYER THICKNESSES')

4 FORMATdX, 'KS= ',F10.8,' KC= ',F10.8,' KI= ',F10.8)

5 FORMATdX, 'LS= ',F10.5,' LC= ',F10.5,' LI= ',F10.5)

6 F0RMAT(/1X, 'UPPER SUMMATION LIMITS ',2X,' NUP=',I5,

( 1 ' MUP=',I5/)

7 F0RMAT(//1 X, 'NUMBER OF HEAT SOURCES= ' ,15)

8 F0RMAT(/1X, 'COORDINATES, LENGTHS, AND WIDTHS OF HEAT SOURCES'/)

9 FORMATdX, 'HEAT SOURCE ' ,7X, ' XSOUR ' ,1 OX, ' YSOUR ' ,9X, ' LXSOUR ' ,

1 9X,'LYS0UR'/)

10 F0RMAT(7X,I3,5X,F10.5,5X,F10.5,5X,F10.5,5X,F10.5)

51
APPENDIX A~TXYZ PROGRAM LISTING

11 FORMATdX, 'POWER DENSITY=' ,F1 2.6)

12 F0RMAT(/1X, 'CALCULATING ',13,' X POINTS WITH A FIRST POINT OF '

1 ,F5.1,' AND A STEP SIZE OF ',F5.1)

13 F0RMAT(/1X,'THE CONSTANT X COORDINATE IS ',F5.1)

14 F0RMAT(/1X, 'CALCULATING ',13,' Y POINTS WITH A FIRST POINT OF '

1 ,F5.1,' AND A STEP SIZE OF ',F5.1)

15 F0RMAT(/1X,'THE CONSTANT Y COORDINATE IS ',F5.1)

16 F0RMAT(/1X, 'CALCULATING ',13,' Z POINTS WITH A FIRST POINT OF '

1 ,F5.1,' AND A STEP SIZE OF ',F5.1)

17 F0RMAT(/1X,'THE CONSTANT Z COORDINATE IS ',F5.1)

27 F0RMAT(1X,'LX= ',F7.2,3X,' LY= ',F7.2)

22 F0RMAT(1X,3F8.2,F10.4)

31 F0RMAT(1X,6I7)

51 F0RMAT(I4,3X,I4)

52 F0RMAT(F10.5,3X,F10.5)

53 FORMAT(II)

54 FORMAT(I4,3X,F10.5,3X,F10.5)

55 FORMATCFIO.S)

56 F0RMAT(I2,3X,F10.5)

57 F0RMAT(F10.5,3X,F10.5,3X,F10.5,3X,F10.5)

CC ********************************************** *********^

C DATA INPUT SECTION (DATA READ FROM F0R010)

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

c IXFLG = 0

lYFLG = 0

IZFLG = 0

ILX=0

ILY=0

ILZ=0

STEPX=1.0

STEPY=1.0

STEPZ=-1 .0

WRITE(11,1

)

WRITE(11,2)

WRITE(11,3)

C RLX AND RLY ARE THE X AND Y DIMENSIONS OF THE RECTANGULAR STRUCTURE

READ(10,52)RLX,RLY

WRITE(11,27)RLX,RLY

C RLS IS THE THICKNESS AND RKS IS THE THERMAL CONDUCTIVITY OF THE TOP LAYER

READ(10,52)RLS,RKS

C RLC IS THE THICKNESS AND RKC IS THE THERMAL CONDUCTIVITY THE MIDDLE LAYER

READ(10,52)RLC,RKC

C RLI IS THE THICKNESS AND RKI IS THE THERMAL CONDUCTIVITY THE BOTTOM LAYER

READ(10,52)RLI,RKI

WRITE(11,5)RLS,RLC,RLI

WRITE (11, 4) RKS, RKC, RKI

C NUP IS THE UPPER LIMIT OF THE SUMMATION OVER THE INDEX N (X-DIR)

C MUP IS THE UPPER LIMIT OF THE SUMMATION OVER THE INDEX M (Y-DIR)

READ(10,51)NUP,MUP

WRITE(11,6)NUP,MUP

52
APPENDIX A~TXYZ PROGRAM LISTING

NUP = NUP + 1 MUP = MUP + 1 IF (NUP. GT. 501 .OR. MUP. GT. 501) GO TO 3999

888 FORMATCIX,' YOUR UPPER LIMIT OF SUMMATION IS TOO LARGE. TRY AGAIN')

READ(10,53)IEDGEX

C SET IEDGEX=0 TO GET CONSTANT POINT FOR X

C SET IEDGEX=1 TO SET THE LIMITS OF THE CALCULATION ALONG X IF(IEDGEX.EQ.O) GO TO 115

READ (10,54) ILX,X1,STEPX

WRITE(11,12)ILX,X1,STEPX

GO TO 125

115 READ(10,55)X1

WRITE(11,13)X1

IXFLG = 1 125 READ(10,53)IEDGEY

C SET IEDGEY=0 TO GET CONSTANT POINT FOR Y

C SET IEDGEY=1 TO SET THE LIMITS OF THE CALCULATION ALONG Y IF(IEDGEY.EQ.O) GO TO 135

READ (10,54) ILY,Y1,STEPY

WRITE (11,1 4) ILY,Y1,STEPY

GO TO 150

135 READ(10,55)Y1

WRITE(11,15)Y1

lYFLG = 1 150 READ(10,53)IEDGEZ

C SET IEDGEZ=0 TO GET CONSTANT POINT FOR Z

C SET IEDGEZ=1 TO SET THE LIMITS OF THE CALCULATION ALONG Z IF (lEDGEZ.EQ.O) GO TO 160

READ (10,54)ILZ,Z1,STEPZ

WRITEd 1,16)ILZ,Z1 ,STEPZ

Z1 =-1.0*Zl

STEPZ=-1 .0*STEPZ

GOTO 175

160 READ(10,55)Z1

WRITE(11,17)Z1

Z1 = -1.0*Z1

IZFLG = 1 ILX = ILX + 1 ILY = ILY + 1 . ILZ = ILZ + 1 175 READ(10,56)NS0UR,P0

C PO 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 20 ARE POSSIBLE)

WRITE(11,7)NS0UR

WRITE(11,11) PO

WRITE(11,8)

WRITE(11,9)

C THE NEXT LOOP READS IN THE COORDINATES OF THE ORIGIN OF THE

C HEATING ELEMENTS ALONG WITH THEIR LENGTHS AND WIDTHS

DO 100 1=1, NSOUR

READ(10,57)XS0UR(I),RLXS0R(I),YS0UR(I),RLYS0R(I)

53
APPENDIX A~TXYZ PROGRAM LISTING

C XSOUR(I) IS THE X COORDINATE OF THE ORIGIN OF I-TH HEATER ELEMENT

C RLXSOR(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 RLYSOR(I) IS THE LENGTH OF THE I-TH HEATER ALONG THE Y DIRECTION

WRITE(11,10)I,XS0UR(I),YS0UR(I),RLXS0R(I),RLYS0R(I)

100 CONTINUE

WRITE(11,1 ) C

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

C

C END OF DATA INPUT SECTION

C END OF INPUT SECTION. THE THERMAL CONDUCTIVITIES OF THE

C SEMICONDUCTOR, CONDUCTOR, AND INSULATOR, THE THICKNESSES OF THE

C SEMICONDUCTOR, CONDUCTOR, AND INSULATOR, THE X AND Y DIMENSIONS

C OF THE CHIP, THE NUMBER OF HEATING SOURCES AND THEIR X,Y AND

C LENGTH AND WIDTH HAVE BEEN ENTERED.

C

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

c 180 PI=3. 14159265

P04LK = 4.0 * PO / (RLX*RLY*RKS)

PILX = PI / RLX

PILY = PI / RLY

C

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

C CALCULATE THE X(0:INTX), Y(0:INTY), AND Z(0:INTZ) ARRAYS

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

c X(1)=X1

DO 200 1=2, ILX

X(I)=X(I-1)+STEPX

200 CONTINUE

Y(1)=Y1

DO 220 1=2, ILY

Y(I)=Y(I-1)+STEPY

220 CONTINUE

Z(1)=Z1

DO 240 I=2,ILZ

Z(I)=Z(I-1)+STEPZ

240 CONTINUE

CC ***************************************************************

C

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) 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).

C 54
APPENDIX A~TXYZ PROGRAM LISTING

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

c

cc *****************************************************************

C CALCULATE THE FOURIER COMPONENTS OF THE HEAT SOURCES, U(N,M)

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

c DO 300 MM=1,MUP

M = MM - 1 DO 250 NN=1,NUP

N = NN - 1 ARUZER(NN,MM)=UZERO(N,M)

250 CONTINUE

300 CONTINUE

C

c *****************************************************************

c

C END OF U(N,M) CALCULATION AND BEGINNING OF MAJOR LOOP FOR Z

Cc *****************************************************************

c DO 3000 IZ=1,ILZ

c

c *****************************************************************

C CALCULATE THE Z DEPENDENT POTION, I.E., FUNZ(N,M,Z)

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

c 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 DO 3000 IY=1,ILY

DO 700 MM=1,MUP

M = MM - 1 COSYT(MM)=COS(FLOAT(M)*Y(IY)*PILY)

700 CONTINUE

C DO 3000 IX=1,ILX

SUM=0.0

DO 1900 MM=1,MUP

M = MM - 1 DO 1700 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) *P ILX) * COSYT(MM)

55
APPENDIX A- -TXYZ PROGRAM LISTING

B0TT0M=(NDN+1 )*(NDM+1)

TSUM=TOP/BOTTOM

SUM=SUM+TSUM

1700 CONTINUE

1900 CONTINUE

TEMP = P04LK * SUM

WRITE (20,22)X(IX),Y(IY),Z(IZ),TEMP

3000 CONTINUE

GO TO 4000

3999 WRITE(6,888)

4000 STOP

END

C **************************************************************^

C END OF THE MAIN PROGRAM

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

56
APPENDIX A~TXYZ PROGRAM LISTING

FUNCTION UZERO LISTING

FUNCTION UZERO(N,M)

DIMENSION XSOUR(20),YSOUR(20),RLXSOR(20),RLYSOR(20)

COMMON RKS,RKC,RKI,RLX,RLY,RLS,RLC,RLI

COMMON NSOUR,XSOUR,YSOUR,RLXSOR,RLYSOR

C

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

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)=1 (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 ***************************************************************

C PI=3.U159265

UZERO=0.0

DO 500 I=1,NS0UR

IF(N.EQ.O) GO TO 100

TERMX = SIN(FLOAT(N)*PI*(XSOUR(I) + RLXSOR (I) ) /RLX)

1 - SIN(FLOAT(N)*PI*XSOUR(I)/RLX)

TERMX=TERMX*RLX/ (FLOAT(N)*PI)

GO TO 150

100 TERMX=RLXSOR(I)

150 IF(M.EQ.O) GO TO 200

TERMY = SIN(FLOAT(M)*PI*(YSOUR(I) + RLYSOR (I) ) /RLY)

1 - SIN(FLOAT(M)*PI*YSOUR(I)/RLY)

TERMY=TERMY*RLY/(FLOAT(M)*PI)

GO TO 250

200 TERMY=RLYSOR(I)

250 TERMI=TERMX*TERMY

UZERO=UZERO+TERMI

500 CONTINUE

RETURN

END

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

57
APPENDIX A~TXYZ PROGRAM LISTING

FUNCTION FUNZ LISTING

FUNCTION FUNZ(N,M,Z)

DIMENSION XS0UR(20),YS0UR(20),RLXS0R(20),RLYS0R(20)

COMMON RKS,RKC,RKI,RLX,RLY,RLS,RLC,RLI

COMMON NSOUR,XSOUR,YSOUR,RLXSOR,RLYSOR

C

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

CC 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 1, I.E., IF Z FALLS IN THE TOP, MIDDLE

C OR BOTTOM LAYER OF THE STRUCTURE. IN ADDITION, SPECIAL CARE IS

C TAKEN AS 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

c *****************************************************

c PI=3. 14159265

GAMMA=SQRT((FL0AT(N)*PI/RLX)**2 + ( FLOAT(M) *PI/RLY) **2)

VS=GAMMA*RLS

VC=GAMMA*RLC

VI=GAMMA*RLI

VT=GAMMA*(RLS+Z)

VM=GAMMA*(RLS+RLC+Z)

VB=GAMMA*(RLS+RLC+RLI+Z)

B0T1=TANH(VS)*TANH(VI)

B0T1=B0T1+(RKI/RKC)*TANH(VS)*TANH(VC)

B0T2=(RKC/RKS)*TANH(VI)*TANH(VC)+(RKI/RKS)

GFUNC=1 .0/(B0T1+B0T2)

AZ=ABS(Z)

IF (AZ.GT.RLS) GO TO 500

58
APPENDIX A--TXYZ PROGRAM LISTING

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

C TOP LAYER CALCULATION

C THIS PORTION IS THE TOP LAYER CALCULATION WHICH IS DEFAULTED

C TO IF Z FALLS INTO THE TOP LAYER

C IF (GAMMA. EQ. 0.0) GO TO 100

TERMS1=TANH(VI)+(RKI/RKC)*TANH(VC)

TERMS2=(RKC/RKS)*TANH(VT)*(TANH(VI)*TANH(VC)+(RKI/RKC)

)

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

C1=2.0*EXP(-VS)

GO TO 20

10 C1=1 .0/COSH(VS)

20 CONTINUE

IF (VT.LT.5.0) GO TO 30

C2=0.5*EXP(VT)

GO TO 40

30 C2=C0SH(VT)

40 CONTINUE

FUNZ=GFUNC*TERMS*C1*C2/GAMMA

RETURN

80 FUNZ=GFUNC*TERMS*EXP(GAMMA*Z) /GAMMA

RETURN

90 FUNZ=GFUNC*TERMS/GAMMA

RETURN

100 FUNZ=(RLS+Z)+(RKS/RKC)*RLC+(RKS/RKI)*RLI

RETURN

500 TOTAL=RLS+RLC

IF (AZ.GT. TOTAL) GO TO 1500

59
APPENDIX A—TXYZ PROGRAM LISTING

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

C MIDDLE LAYER CALCULATION

C THIS IS THE MIDDLE LAYER CALCULATION WHICH IS DEFAULTED TO IF

C Z FALLS INTO THIS DOMAIN OF DEPTHS

C IF (GAMMA. EQ. 0.0) GO TO 1000

TERMC=TANH(VI)+(RKI/RKC)*TANH(VM)

IF (VS.GT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0) 60 TO 800

IF (VS.LT.5.0) GO TO 250

C1=2.0*EXP(-VS)

60 TO 260

250 C1=1.0/C0SH(VS)

260 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=C0SH(VM)

330 CONTINUE

FUNZ=GFUNC*TERMC*C1*C2*C3/GAMMA

RETURN

800 FUNZ=6FUNC*TERMC*2.0*EXP(6AMMA*Z) /GAMMA

RETURN

1000 FUNZ=(RKS/RKI)*RLI+(RKS/RKC)*(RLS+RLC+Z)

RETURN

60
APPENDIX A—TXYZ PROGRAM LISTING

************ ******************************************

BOTTOM LAYER CALCULATION

THIS IS THE BOTTOM LAYER CALCULATION WHICH IS USED IF Z FALLS

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)G0 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)G0 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)G0 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

C1=2.0*EXP(-VS)

GO TO 1560

1550 C1=1 .0/COSH(VS)

1560 CONTINUE

IF (VC.LT.5.0) GO TO 1570

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*C1*C2*C3*C4/GAMMA

RETURN

1900 FUNZ=GFUNC*4.0*EXP(GAMMA*Z) /GAMMA

RETURN

2000 FUNZ=(RKS/RKI)*(RLS+RLC+RLI+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) /(6AMMA*C0SH(VC)

RETURN

2600 FUNZ=GFUNC*2.0*EXP(VB-VC-VI) /(GAMMA*C0SH(VS))

RETURN

END

61
NBS-n4A iREv. 2-8C) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBS SP 400-76

2. Performing Orjzui. Report No 3. Publication Date ^ril 1984

4. TITLE AND SUBTITLE

Semioonduotor Measurement Technotogy : TXYZ: A Program for Semiconductor IC

Thermal Analysis

5. AUTHOR(S)

J. Albers

6. PERFORMING ORGANIZATION f/f )o/nt or other than NBS. see /n struct/on s; NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

Final

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

Same as item 6.

10. SUPPLEMENTARY NOTES

Library of Congress Catalog Card Number: 84-601026

I I Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here)

A computer program, TXYZ, for the thermal analysis of semiconductor integrated

circuits is presented and its applications are discussed. The program makes use of

the closed form, analytic solution of the steady-state heat flow problem for a rectangular three-layer structure with multiple heat sources on the top layer. The

temperature may be obtained for any point or set of points in the structure and is

useful in the determination of the stea5y-state thermal response of IC chips and

packages

12. KEY WORDS (S/x to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

FORTRAK; Fourier analysis; integrated circuit; semiconductor; steady-state heat

flow; thermal analysis; themal resistance.

13. AVAILABILITY

fxl Unlimited

1 1 For Official Distribution. Do Not Release to NTIS

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

20402.

Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

65

15. Price

USCOMM-DC 6043-P80
Announcement of Semiconductor Measurement Technology

List of Publications 72 - 1962-1982

Semiconductor Materials and Processes Division

National Bureau of Standards

Bldg. 225, Room A305

Washington, DC 20234

Please send a copy of your latest "Semiconductor Measurement Technology, List

"of Publications 72."

Name

Company

Address

City State Zip Code

NBS TECHNICAL PUBLICATIONS

PERIODICALS

JOURNAL OF RESEARCH—The Journal of Research of the National Bureau of Standards reports NBS research and develop- ment in those disciplines of the physical and engineering sciences in which the Bureau 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 lime are survey articles on topics

closely related to the Bureau's technical and scientific programs. As a special service to subscribers each issue contains complete

citations to all recent Bureau publications in both NBS and non- NBS media, issued six times a year. Annual subscription: domestic

$18; foreign $22.50. Single copy, $5.50 domestic; $6.90 foreign.

NONPERIODICALS

Monographs—Major contributions to the technical literature on

various subjects related to the Bureau's scientific and technical activities. Handbooks—Recommended codes of engineering and industrial practice (including safely codes) developed in cooperation with in- terested industries, professional organizations, and regulatory bodies

Special Publications—Include proceedings of conferences sponsored by NBS, NBS annual reports, and other special publications

appropriate to this grouping such as wall charts, pocket cards, and

bibliographies. Applied Mathematics Series—Mathematical tables, manuals, and

studies of special interest to physicists, engineers, chemists,

biologists, mathematicians, computer programmers, and others engaged in scientific and technical work.

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 NBS under

the authority of the National Standard Data Act (Public Law

90-396).

NOTE; The principal publication outlet for the foregoing data is the Journal of Physical and Chemical Reference Data (JPCRD)

published quarterly for NBS by the American Chemical Society

(ACS) and the American institute of Physics (AiP). Subscriptions,

reprints, and supplements available from ACS, 1 155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical irformalion developed at the Bureau 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 NBS 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. NBS administers this program as a supplement to the activities of the private sector standardizing

organizations. Consumer Information Series—Practical information, based on NBS research and experience, covering areas of interest to the con- sumer. Easily understandable language and illustrations provide

useful background knowledge for shopping in today's tech- nological marketplace. Order ihe above NBS publications from: Superintendent of Documents. Government Printing Office, IVashington, DC 20402. Order the following NBS publications—FIPS and NBSIR's—from

the National Technical Information Service . Springfield, VA 22161

.

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 Govern- ment regarding standards issued by NBS 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 1 1717 (38 FR 12315, dated May II, 1973) and Part 6

of Title 15 CFR (Code of Federal Regulations).

NBS Interagency Reports (NBSIR)—A special series of interim or

final reports on work performed by NBS for outside sponsors

(both government and non-government). In general, initial dis- tribution is handled by the sponsor; public distribution is by the National Technical Information Service , Springfield, VA 22161,

in paper copy or microfiche form.
U.S. Department of Commerce

National Bureau of Standards

Washington, D.C. 20234

Official Business

Penalty for Private Use S300