PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES totum on or baton date duo.

 

DATE DUE DATE DUE DATE DUE

__"___l
_:
ll '
4C:
= II____

JLJE IDE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

—
-

 

MSU ls 5n Affirmative ActioNEqual Opportunity Institution

SYMBOLIC, ALGEBRAIC, AND NUMERIC SOLUTIONS TO HEAT CONDUCTION

PROBLEMS USING GREEN'S FUNCTIONS

BY

Paul Henry Zang

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

1987

O

’7
4-"

@763

W”)

°Copyright by Paul Henry Zang 1987
All Rights Reserved

SYHJ

QM:

Paul Hen:
Doctor 01

 

 

 

This is to certify that the
dissertation entitled

SYMBOLIC, ALGEBRAIC, AND NUMERIC SOLUTIONS TO HEAT CONDUCTION
PROBLEMS USING GREEN'S FUNCTIONS

presented by

Paul Henry Zang

3- \ ‘2.- 81

Date

   
  

II
Doctor of Philosophy andidate

has been accepted towards fulfillment
of the requirements for

Doctor of Philosophy Degree in Mechanical Engineering

Michigan State University

by

qr/JeLwV/‘I- A ”7&0 // /W‘/

Jame V. Beck Date
P essor

D ssertation Advisor

Department of Mechanical Engineering

@flwflm¢( é '1 #5 7
‘i‘fgsi‘liifi Chaim ”“8

Department of Mechanical Engineering

 

 

SYMBL

 

..
Hf:

 

ABSTRACT

SYMBOLIC, ALGEBRAIC, AND NUMERIC SOLUTIONS TO HEAT CONDUCTION

PROBLEMS USING GREEN'S FUNCTIONS

by
Paul Henry Zang

Symbolic calculations that involve tedious, error-prone evaluation
have plagued the scientist and engineer for many centuries. A new tool
called computer algebra can be used to evaluate complex mathematical
operations, such as differentiation and integration, which can be
repetitious when applied to partial differential equations. Symbolic
results from computer algebra systems can offer insight to problems
which numerical results lack.

Symbolic manipulation of expressions and operations are used
extensively in this thesis. A new technique, using computer algebra,
for the symbolic solution of heat diffusion-type problems is examined.
The new technique involves the Green's function approach and uses
Laplace transforms and separation of variables for determining
appropriate Green's functions. Data bases of Green's functions and
integrals are used to speed up the calculation time of solutions.

A systematic and orderly procedure for developing Green's functions
which are computationally efficient for small dimensionless times is
examined. The small time Green's functions are used in a partitioning

scheme which accelerates the evaluation time of the symbolic solutions.

tempera
problem.
for var:

energy 1

examples
literatu

effiCien.

 

 

Two computer programs are presented that symbolically calculates
temperature distributions for a limited number of heat transfer
problems. The one dimensional program called CANSS generated solutions
for various types of boundary conditions, initial conditions, and volume
energy heat sources. The two dimensional program called CANSSZD
generates temperature distributions for boundary conditions of the
zeroth, first and second kinds. The temperature distributions of
examples presented in this thesis match solutions found in the
literature and are partitioned in time to increase evaluation

efficiency.

ACKNOWLEDGMENTS

The author wishes to express his sincere appreciation and gratitude
to his major thesis adviser, Dr. James V. Beck, for his insight, con-
tributions and painstaking editing of the rough manuscript. His
participation in the dissertation process provided the author with a
true appreciation of what it means to be a research scientist.

The author also wishes to thank the members of the guidance com-
mittee. Professors Steven Shaw and Craig Somerton gave the proper
perspective to the dissertation process. Dr. Richard Phillips deserves
special recognition for his cool and collected thoughts and eternal
patience. Thanks are due also to Dr. Kevin Cole for his valuable in-
puts.

Partial financial support for this dissertation was provided by a
grant from the National Science Foundation, Grant No. MBA 81-21499, the
Department of Mechanical Engineering, and the Department of Engineering
Research. The author extends his appreciation to the agency and his
gratitude to the departments.

Many thanks are due to my parents, Jack and Wanda Zang, for their
never ending confidence in me, my wife’s parents, Jerry and Joanne
McCarthy, for their patience and trust, and all my friends, particularly
Maureen Clements and Brian Agar.

Last, to my wife, Mary Ellen, the author dedicates this disserta-
tion for her unflagging confidence and optimism in me. Through her,

this dissertation was made possible.

LIST OF I

LIST OF P

LIST OF S

r-ql

~- Gas:

 

A)

(AJ

 

 

 

TABLE OF CONTENTS

LIST OF TABLES ............................................ ix
LIST OF FIGURES ........................................... x
LIST OF SYMBOLS ........................................... xiii
1. INTRODUCTION .......................................... 1
1.1 Previous work ................................ 5
1.2 Thesis Objectives ............................ 9
1.3 Synopsis of Thesis ........................... 11
2. GREEN'S FUNCTION FORMULATION .......................... 13
2.1 Introduction ................................. 13

2.2 Mathematical Derivation of Heat Diffusion in

a Body ....................................... 16
2.2.1 Heat Conduction Equation and Boundary
Conditions .............................. 16

2.2.2 The Auxiliary Green's Function Equation . 25

2.3 The Green's Function Approach ................ 28
2.3.1 Mathematical Derivation of the Green's

Function Approach ....................... 28

2.3.2 Determination of the Green's Function ... 36

2.3.3 Products of Green's Functions ........... 39

2.4 Formalism of the Green's Function Approach ... 44

2.5 Summary ...................................... 48

3. SMALL TIME GREEN'S FUNCTIONS OBTAINED USING LAPLACE
TRANSFORMS ............................................ 50

3.2 Mathematical Development of the Small Time
Green's Function ............................... 57

3.3 Green's Functions for Some Semi-infinite Cases
in One Dimension ............................... 73

3.4 Small Time Green's Functions for Finite Bodies . 82
3.4.1 Slab Insulated On Both Boundaries (X22).. 92

3.4.2 Slab With a Carslaw Left Boundary
Condition and Right Boundary Insulated

(X42) ................................... 96

3.5 Summary ........................................ 99
TRANSIENT ONE DIMENSIONAL CANSS PROGRAM ............... 101
4.1 Introduction ................................... 101
4.2 One Dimensional CANSS Program .................. 102
4.3 One Dimensional CANSS Examples ................. 109

4.3.1 Semi-infinite Example Problem (X2OB1T0).. 110

4.3.2 Finite Slab With a Transendental
Initial Condition (X21B11T6) ............. 114

4.3.3 Finite Slab With a Boundary Condition
a Function of the Square Root of Time

(X22B3OTO) .............................. 119

4.4 Some Integrals Used in One Dimensional
Problems ....................................... 125
4.4.1 The Dawson Integral ..................... 125

4.4.2 An Exponential Integral in One
Dimensional Problems .................... 129

4.5 Time Region Partitioning for One Dimensional

Problems ....................................... 131

4.6 One Dimensional CANSS Flowchart/Example ........ 134
4.7 Summary ........................................ 136
TRANSIENT TWO DIMENSIONAL CANSSZD PROGRAM ............. 138
5.1 Introduction ................................... 138
5.2 Transient Two Dimensional CANSSZD Program ...... 140
5.2.1 CANSSZD Program ......................... 140

5.2.2 Two Dimensional CANSSZD Flowchart/
Example ................................. 144

vii

 

 

APPEjm IX
Ufimm

LIST OF 5

 

 

5.3 Two Dimensional CANSSZD Examples ............... 146

5.3.1 A Two Dimensional Plate with Heating .... 146
5.3.3 A Partially Heated Plate ................ 151
5.4 Two Dimensional Problems ....................... 153
5.4.1 An Integral in Time Region One .......... 155
5.4.2 An Integral in Time Region Two .......... 156

5.5 Time Partitioning in Two and Three Dimensions .. 163

5.6 Summary ........................................ 169
6. SUMMARY AND CONCLUSIONS ............................... 170
6.1 Summary ........................................ 170
6.2 Conclusions .................................... 173
6.3 Recommendations ................................ 178
APPENDIX A. A SHORT TABLE OF LAPLACE TRANSFORMS .......... 182
APPENDIX B. SOME USEFUL INTEGRALS ........................ 184
APPENDIX C. CANSS AND CANSSZD PROGRAM EXAMPLES ........... 200
LIST OF REFERENCES .......................................... 211

viii

 

“‘1'“

. r-v2;

3.5...“133 "+4-

 

Iable 3.

Table 3.

Table 3.

Table 3.

Table 3.:

Table 3_

Table 3.

(7)

6

 

 

Table

Table

Table

Table

Table

Table

Table

Table

Table

Table

LIST OF TABLES

Small Time for Finite and Infinite Body Source at

x' - 0, Point of Interest at x - 0.25 ........... 53

Exponential as a Function of Time ............... 53
-q(x+x')

Coefficients of e ........................ 68
-q(2L-x-x')

Coefficients of e ..................... 68
-q(2L+x-x') -q(2L-x+x')

Coefficients of e and e .... 69

-qx
Inverse Laplace Transforms of A(-) e .......... 72

X42 Case Cl - 0.10
Small and Large Time Green’s Function for Various
Number of Terms ................................. 98

Small Time Green's Function for a(t-r)/L 5 0.025. 100

Key for Heat Conduction Data Base (Excerpted from
Tzeng and Beck [1985]) .......................... 105

Dimensionless Change in Temperature from Equation

(5.3.10) When the Partition Time is Set to 0.1 and
the x Coordinate is Set to L/2 .................. 154

ix

Figure
Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

2.2a

LIST OF FIGURES

Structure of an Expert System ................... 4
Thermal Energy Wave With Exponential Decay ...... 15
Semi-infinite Slab With a Temperature Boundary
Condition (X10) ................................. 21
Semi-infinite Slab With a Heat Flux Boundary
Condition (X20) ................................. 21
Semi-infinite Slab With a Convective Boundary
Condition (X30) ................................. 21
Finite Slab With a Non-convective Thin Film and a
Heat Flux Condition (X42) ....................... 24

Finite Slab With a Convective Thin Film and a Heat

Flux Condition (X52) ............................ 24

A Semi-infinite Body and an Infinite Body ....... 26

A Description of Normals at the Surface of a Finite
Body ............................................ 32

Reflections of Sources and Sinks in a Finite Body 38

Distinct Cases of Green's Function for One

Dimensional Problems ............................ 46
Small Time Green's Functions for Various Lengths

of Finite Body Versus Time ...................... 52
Reflections of Sources and Sinks in a Finite Body
and the Locations of Additional Reflection Terms. 66

Green's Function Versus Time for a Semi-infinite
Body With Boundary Conditions of the Zeroth, First,

and Second Kind ................................. 75
Green's Function Versus Time for a Semi-infinite
Body With a Boundary Condition of the Third Kind
and Various Values of the Parameters ............ 77
Green's Function Versus Time for a Semi-infinite
Body With a Boundary Condition of the Fourth Kind
and Various Values of the Parameters ............ 78

Green's Function Versus Dimensionless Position for

an Infinite Body ................................ 79

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

.6b

.6c

.7a

.7b

.7a

.8a

.8b

.8c

.9a

.9b

.9c

.10a

.10b

.11

.1a

.1b

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
First Kind ...................................... 80

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Second Kind ..................................... 81

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Third Kind, B1 - 0.1 ............................ 83

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Third Kind, B1 - 1.0 ............................ 84

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Third Kind, Bl - 10.0 ........................... 85

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Fourth Kind, Cl - 0.1 ........................... 86

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Fourth Kind, C1 - 1.0 ........................... 87

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Fourth Kind, C1 - 10.0 .......................... 88

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Fifth Kind, C1 - 0.1, Bl - 0.1 .................. 89

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Fifth Kind, Cl - 0.1, Bl - 1.0 .................. 90

Green's Function Versus Dimensionless Position for
an Finite Body With a Boundary Condition of the
Fifth Kind, Cl - 0.1, Bl - 10.0 ................. 91

Finite Body Insulated on Two Sides .............. 93

Finite Body Insulated on One Side and With a Non-
convective Thin Film on the Opposite Side ....... 93

Number of Terms for Convergence of the Green's
Function Versus Dimensionless Time for an Insulated
Body (X22). Accuracy - 0.00001 ................ 95

Semi-infinite Body With Constant Heat Flux on the
Surface (X20B1T0) ............................... 111
Semi-infinite Body With Constant Heat flux on the
Surface and a Transcendental Initial Temperature

(X22811T6) ...................................... 111

xi

 

Figur

Figure

Figure

&'
r3
H
(D
U!

I.“
"a
'1
(D
LII

H
U?
r.
*1
(D
Ln

 

 

 

 

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

.1c

.4a

.4b

.4c

.1a

.1b

.6a

.6b

Semi-infinite Body With a Heat Flux Condition as a
Function of Time (X22B30T0) ..................... 111

Dimensionless Temperature Versus Dimensionless Time
for a Semi-infinite Slab With Constant Heat Flux 113

Integral Error Function as a Function of n ...... 115

Dimensionless Temperature Versus Dimensionless Time
for a Finite Body in Example Two for the Boundary
Condition ....................................... 118
Dimensionless Temperature Versus Dimensionless Time

for a Finite Body in Example Two for the Initial
Condition ....................................... 120
Dimensionless Temperature Versus Dimensionless Time

for a Finite Body in Example Two for the Initial

and Boundary Condition Added Together ........... 121
Two Types of Approximations for Small Time

Integrals ....................................... 124
Dawson Integral ................................. 128
Convolution of the Green's Function and a Function
of Time ......................................... 133
Thin Plate With Heat Flux and Temperature

Conditions ...................................... 147
Partially Heated Thin Plate ..................... 147

Integral in Equation (5.4.4) Versus Dimensionless
Time When X - 0 ................................. 157
Integral in Equation (5.4.5) Versus Dimensionless
Time When Y - 0 ................................. 158
The cgshe Function Versus Dimensionless Time When

C1 - a and Various Values of C2 ................. 161

The coshe Function Versus Dimensionless Time When
C1 - x and Various Values of C2 ................. 161

Time Partitioning for an Aspect Ratio Greater Than
One for a Two Dimensional Body .................. 165

Time Partitioning for an Aspect Ratio Less Than
One for a Two Dimensional Body .................. 165

Time Partition Regions for a Three Dimensional Body

4.
> 1, r
zx

Where r+ > 1, and r:y >1 .............. 168

xii

S and s'

 

 

 

 

 

 

 

 

Q5.°£xifit

 

m
5
O.
U)

xuz»a<n:nxo:z:zt~7<cap4::c>~:c:c>u:>»x 8 < c n'm n.0'o a B WW*i*:TGIFhO'O-O O‘m
’~<:
N

}<
N

LIST OF SYMBOLS

Parameter in Equation (3.2.37)
Thin Film Thickness

Specific Heat

Parameter in Equation (3.2.37)
Exponential Function

Forcing Function

Volume Energy Heat Source
Heat Transfer Coefficient
Index

Index

Thermal Conductivity

Index

Index

Laplace Transform Parameter
Heat Flux

Spatial Coordinate

Number of Boundaries

Time

Instantaneous Source Solution
Volume Element

Laplace Solution

Length Variables

Function for Equation (3.2.29)
Function for Equation (3.2.29)
Constants

Constants

Initial Condition

Green's Function

Heavyside Function

Index for Numbering System
Index for Numbering System
Kernel Function

Length of Finite Slab

Matrix in Equation (3.2.20)
Matrix in Equation (3.2.20)
Heat Source Strength

Region of Space

Surface Boundary

Temperature

Transformation Variable
Numbering System Parameters

GREEK SYMBOLS

ym'S-TDQ

Thermal Diffusivity
Eigenvalue

Start-up Time

Small Number

Dummy Variable

xiii

Jpfligffrov 0.001nlon

 

Superscr.

Subscript

 

 

 

 

 

 

Square Root Sign
Density

Constant

Dummy Variable
Dummy Time Variable
Constant

Delta Function
Laplacian Operator
Infinity

Partial Derivative Sign
Large Time Index
Gamma Function

Pi (3.1415....)

aflhmaqmmamdbk

Superscripts

* Dimensionless
+ Dimensionless

subscripts

L Normalized With Respect to L

xiv

CHAPTER 1

INTRODUCTION

The first introduction of a child to mathematics is frequently
symbolic. The child learns that symbols have meaning - such as, two
crossed lines mean the numbers are added and two horizontal lines means
what follows is the sum total. Calculus brings the onset of the numeri-
cal approach and the acceptance of symbols becomes lost in the quest for
approximate solutions.

Computers, long used to reduce the repetition and increase the
accuracy of numerical calculation, can now do the same for symbolic
analysis. Computer algebra systems are used to define new mathematical
concepts and increase the speed of repetitive symbolic analysis pre-
viously performed by hand. Symbolic analysis can give insight to the
structure of the physics of problems that purely numerical solutions
miss.

Computer algebra programs work best on algorithmic representations
of systems that involve complex mathematical operations, such as in-
tegration or differentiation, but are straightforward. Symbolic
calculations of this type are said to be "fierce". At another extreme,
computer algebra programs also work well on systems that represent a
broad class of problems and can be expressed using an algebraic struc-
ture. The calculations for this case are not "fierce", but tedious and

repetitive.‘

Symbolic manipulation, or computer algebra, was brought into the
public domain by a program called MACSYMA [The MATHLAB Group, 1983] in
the early 70's. Since that time, many computer algebra systems have
been brought to the market. The early versions of computer algebra
programs were expensive, memory intensive, and dependent on specific
computer hardware. Today, the computer algebra software is more user
friendly, less machine dependent, and less expensive.

The introduction of symbolic software compatible with micro-
computers, personal computers, and, more recently, hand held
programmable calculators encourages the use of computer algebra and will
ultimately lead to educational programs in symbolic manipulation.
Computer algebra can be used as an educational tool in many branches of
science including calculus, physics, chemistry, and engineering. It is
now possible for most universities to offer computer algebra to students
as a learning aid.

This thesis presents a study of the application of computer algebra
techniques to a field of engineering and can be described as gomputer
aided symbolic engineering (CASE). Symbolic computer programs use
mathematical concepts to describe, analyze, and evaluate the mathemati-
cal operations that formerly required pencil-and-paper analysis.
Symbolic methods improve accuracy by evaluating solution in closed form.
The accuracy of the solutions obtained using computer algebra techniques
is dependent on the accuracy of the parameters or variables input to the
problem thus reducing or eliminating human error caused by evaluating
repetitive algebraic processes. The pattern revealed by a group of
variables or parameters will cover an infinite number of cases.

Computer algebra can be used in the areas of fluid dynamics, to solve
large full symbolic matrices, finite element methods, to generate ac-

curate trial functions, and machine dynamics to name a few. The CASE

field will continue to grow as more uses are discovered for computer

algebra.

Hayes and Michie [1984] state an expert system applies a structured
set of rules to a data base to evaluate input and calculate output, and
is a rudimentary form of artificial intelligence (AI). The expert
system discussed in this thesis employes a rule based procedure
(inference engine) along with known facts and assertions (knowledge
base) and is used to treat problems in heat diffusion. Figure 1.1
represents the structure of a expert system.

An expert system called computer algebraic, numeric, and symbolic
system (CANSS) applies symbolic manipulation to problems in mathematical
physics. The CANSS program developed in this thesis has the capacity to
solve heat diffusion problems not found in the literature. The CANSS
program represents a new technique in CASE for obtaining temperature
distributions for linear, multi-dimensional, transient heat diffusion
problems by the application of symbolic manipulation computer software.
As a tool for heat transfer engineers in the 80's and 90's, symbolic
analysis can be compared to the introduction of finite element tech-
niques to transient heat transfer in the early 70's.

The CANSS program uses a Green's function approach and a symbolic
manipulation program called SMP [1983] to generate analytical tempera-
ture distributions for problems that involve various linear boundary
conditions, initial conditions and volume energy heat sources. The
technique can calculate temperature distributions for some nonlinear
heat diffusion problems by using a transformation found in Ozisik [1980,
pg. 440] using the Kirchhoff transformation. Some nonlinear diffusion
problems can be treated as linear for special combinations of the physi-

cal properties.

Eoum>m uumoxm so we ousuoauum H.H wusmgm

mxm3mz<

llilliillillll

 

 

mzHUzm
mozmmmmzH

 

 

mm<m
NUQMAZOZM

 

 

 

 

 

I

monemmac

xmamrm Hmmmxm

Transient temperature distributions for finite, infinite, and semi-
infinite single or multi-dimensional bodies are typically calculated
using finite difference or finite element methods. These methods gener-
ate information about the temperature for a finite set of points or
nodes and times. Care must be taken with the distribution size of the
points or mesh and the time step to insure stability of the solution.
Solutions obtained by numerical techniques must be re-calculated for the
entire domain of space and time when one or more input parameters
change. Simple substitution of new parameters into the symbolic form of
the solution result in new solutions without extensive re-calculation.

Computer algebra methods offer accurate symbolic results that can
act as test cases for the calibration of purely numerical techniques
such as finite element or finite difference. Specifically, symbolic
solutions to diffusion problems can complement the numerical solutions
that are normally used in heat transfer analysis. Symbolic solutions to
basic geometries can be included as trial functions in the finite ele-
ment and finite difference methods. Symbolic temperature distributions
can be also used as a starting point for numerical procedures or for
procedures where a poor initial guess will lead to solutions that do not

converge or oscillate.

1.1 Previous Herr

The Green's function approach to the solution of partial differen-
tial equations of heat diffusion is well documented. Morse and Feshbach
[1953], and more recently, Ozisik [1980] and Beck [1984a] provide a
structure for determining temperature distributions using the Green's

function approach.

Walters [1949] uses the Green's function approach to solve tran-
sient heat conduction and vibration problems analytically and, more
recently, Hassanein and Kulcinski [1984] use the Green's function ap-
proach to examine the rapid heating of fusion reactor walls. Hassanein
and Kulcinski reported comparisons of the Green's function approach to
the finite difference method. They report the Green's function ap-
proach requires more analytical calculations than the finite element
method but the time step may be larger and the calculations are more
straightforward. The analytical calculations in their study could be
calculated by a computer algebra system to decrease the complexity of
the problem.

The Green's function approach has been used to generate influence,
or kernel, functions for the unsteady surface element method developed
by Keltner and Beck [1981]. The unsteady surface element method splits
up the boundary and uses the influence functions to obtain solutions for
both linear and nonlinear boundary conditions. Cole [1986] used the
Green's function approach to determine some influence functions for
conjugate heat transfer problems in the unsteady surface element method.

The Green's function is geometry specific - a different Green's
function is needed for each geometry and each set of boundary condi-
tions. Many references and texts such as Butkovski [1982], Stakgold
[1979] and Beck [1984b] include lists of the Green's functions. Beck
and Litkouhi [1985] proposed a numbering system, which is used exten-
sively in this thesis, to generate a data base for the Green's
functions.

Numerically inefficient solutions to heat and mass diffusion
problems have been the bane of the heat transfer engineer for many
years. Miller and Gordon [1937] report the solutions obtained using the

traditional methods of Fourier series, while mathematically acceptable,

are inefficient for small regions of time because of the slow rate of
convergence. Aizen, et. a1., [1971] show additional methods of increas-
ing the speed of convergence for certain problems in heat diffusion.

The Green's function approach is well suited for symbolic calcula-
tion because it is algorithmic and uses operations not typically found
in numerical analysis. The capability for integration and differentia-
tion of the computer algebra systems make them the unique vehicle for
calculating temperature distributions using the Green's function method.
Integration and differentiation are the key to the Green's function
approach. The algorithmic structure of the approach is used with a
computer algebra system to generate new symbolic or numeric solutions.

In a recent paper by HaJi-Sheikh and Lakshminarayanan [1986],
symbolic analysis is applied to the solution of diffusion type problems
through the use of the Galerkin method. The temperature solutions
reported in this work are efficient for large dimensionless times and
for complex geOmetrical shapes.

The integration routines in SMP are based on an algorithm first
proposed by Risch [1969]. Risch bases his algorithm on the text by Ritt
[1949] which describes the integration of expressions in finite terms.
The Risch algorithm states that the integral must be represented as an
elementary function and the integral solution must also be expressible
in elementary functions. An elementary function is a function composed
of polynomials, exponentials, and logarithms using only rational and
algebraic operations. The Risch algorithm is still in its infancy, but
the work of Cherry [1986] and Knowles [1986] continue to expand the
functions to which the algorithm is applied. Cherry extended the types
of elementary functions to include some logarithmic integrals and

Knowles extended Cherry's work to include some exponential integrals and

a broader class of logarithmic integrals. Ng [1977] describes proce-
dures based on the Risch algorithm as pattern recognition strategies.

Roach and Steinberg [1984] use symbolic manipulation in the area of
computational fluid dynamics. They report the ability of speeding up
code development time and the prospect of virtually error-free testing
of constitutive equations and difference forms.

Rand [1984] describes the application of a computer algebra program
for solving ordinary differential equations, finding eigensolutions to
eigenvalue problems, and solving some examples of boundary value
problems. Rand reports the application of computer algebra for finding
approximate solutions to differential equations that contain small
parameter by using perturbation methods. Mathematicians and engineers
have been interested in finding closed form expressions for summations.
Moenck [1977] describes a method using symbolic manipulation to express
the sum of a rational function as a rational function part and a
transcendental part involving derivitives of the gamma function.

Char, et. al., [1986] describes the application of computer algebra
to undergraduate mathematics curriculum. They have found it feasible to
offer courses in computer algebra to large groups of students. The
initial findings suggest the introduction of a computer algebra system
to undergraduates has met with limited success. They suggest more
powerful facilities for integration, a smooth interface to numerical

procedures, and a more user friendly interface.

1.2 Thesis Objectives

The first objective of this thesis is to use the computer algebra
system SMP [1983] to develop two computer algebra programs that calcu-
late symbolic temperature distributions, one for one-dimensional bodies
and the second program for two-dimensional bodies. A unified method of
solution based on Green's functions is developed. The symbolic solu-
tions generated by the programs are to be computationally efficient for
the whole range of dimensionless times.

The second objective is the generation of symbolic temperature
distributions for some heat diffusion problems using the programs men-
tioned above. Experience needs to be obtained with the capabilities and
limitations of such computer programs. The solutions to basic heat
transfer problems can be used as kernel functions in numerical tech-
niques such as finite element, finite difference, and boundary element
methods. The programs in this thesis are based on the Green's function
approach to heat transfer which implies a linear problem model.

Litkouhi [1982] shows the distributions developed by the linear Green's
function approach can be used in a numerical surface element method to
obtain distributions to more physically complex or nonlinear problems.
Some temperature distributions that are not available are also gener-
ated.

The third objective of this thesis is to investigate a new class of
applications of computer algebra systems to mathematical problems in
engineering. The programs and examples in this thesis are chosen from
heat transfer, but the computer algebra environment is not restricted to

the diffusion equation.

10

The fourth objective of this thesis is to study the need for a
structure in the field of interest. Structure in a field of interest
means the availability of a formalism for the calculation of solutions
and a data base of pertinent information relating to the field. The
structure is important because cases with different parameters can be
treated with the same formalism. For example, the formalism for solu-
tion of a conduction heat transfer problem could be the application of
Fourier series analysis and the data base would include solutions to
fundamental ordinary differential equations and some basic integrals.

The fifth objective is an examination of the knowledge bases neces-
sary for efficient computation of solutions. Hayes-Roth, et. al.,
[1983] state that the knowledge base aid the efficiency of an expert
system by eliminating "blind alleys", eliminating repetitive calcula-
tions, and applying specific information about the problem.

A sixth objective is to use the new concepts of time partitioning
by Beck and Keltner [1984] to obtain symbolic solutions to heat transfer
problems that converge rapidly for the whole range of dimensionless
times. Efficient symbolic solutions can give valuable insight to the
physics of the solutions. An investigator can determine how the solu-
tion to a specific system would react to a change of input parameters
without re-calculating the solution over the entire domain of time and
space.

The final objective is to survey areas associated with mathematical
physics in engineering for the mathematical structure or formalism that

can effectively use the elements of computer algebra.

 

 

f?

11

1.3 Synopsis of the Thesis

This thesis is divided into six chapters. Chapter 1 is the intro-
duction and gives the motivation for studying heat diffusion using a
Green's function approach and computer algebra. Previous work and the
scope of the thesis are examined in Chapter 1.

Chapter 2 introduces the mathematical development and structure of
the Green' function approach to the solution of heat diffusion problems
that appear in this thesis. The partial differential equations and the
_associated boundary and initial conditions are examined. A procedure
for obtaining the Green's functions and a procedure for generating the
products of one dimensional Green's functions to obtain the Green's
functions for two and three dimensional cases in rectangular coordinates
are presented.

Chapter 3 describes a method for obtaining approximate expressions
for Green's functions in one dimension using the technique of Laplace
transforms. The expressions generated by this method converge quickly
for small times. Examples of Green's functions are given which examine
the effect of various boundary conditions for semi-infinite bodies.
Three example problems are examined for finite, one dimensional bodies
using various boundary conditions.

Chapter 4 describes a computer algebra program, called CANSS, that
calculates symbolic temperature distributions for semi-infinite and
finite bodies in one dimension. Various types of boundary conditions,
initial conditions, and volume energy heat sources may be applied to the

bodies using the CANSS algorithm. Three example problems are examined

12

along with some important integrals that relate to one dimensional
problems. A flowchart/example is presented.

Chapter 5 describes a symbolic computer algorithm, called CANSSZD,
that generates symbolic temperature distributions for two dimensional
plates. Two example problems are examined along with some special
integrals that occur during the calculation. A discussion of time
partitioning in two and three dimensions concludes the chapter.

Chapter 6 presents the conclusions and summary. This chapter also

offer suggestions for the extension of this thesis.

CHAPTER 2

erunq'S FUNCTION FORMULATION
2.1 Introduction

This chapter traces the development of the Green's function ap-
proach to the solution of linear partial differential equations for heat
conduction. This approach has become an accepted solution technique due
to the recent work by Greenberg [1971], Ozisik [1980] and Beck [1984b],
and the previous work of Morse and Feshbach [1953]. These authors show
that the Green's function approach to the solution to heat conduction
problems offers both simplicity and structure and gives the engineer and
scientist an alternate approach for solving diffusion type problems.

The unifying structure implicit in the Green's function approach
provides an ideal testing vehicle for demonstrating the use of symbolic
computation.

The Green's function approach is simple because of the straightfor-
ward manner in which the boundary condition, initial condition and heat
generation terms of the solution are generated. The approach has struc-
ture due to the use of a single Green's function for the complete
solution.

Classical heat conduction theory states that the heat transfer

rate, q, is proportional to the temperature gradient in the medium, or,

~fll
q ari (2.1.1)

13

 

 

 

 

 

 

 

 

 

where
mediun
conduc

heat f

on or
consti

sion 1.

that i]

The mo<

"a5 fix
and She
CauSES
velocit
“WP-S t;
{1958} '
to ,3ch
Perature
Th:
Ve
Small, ec

small r85

14

where ri is a coordinate in the direction of the heat flow in the
medium. Inserting a constant of proportionality called the thermal
conductivity, k, the constitutive equation that describes the rate of

heat flow by conduction in the ri direction as,

q _ -k g¥ . (2 1.2)
1

One consequence of the above equation is that a thermal disturbance
on or in the medium is propagated everywhere instantaneously. The
constitutive equation is very accurate even though instantaneous diffu-
sion is unrealistic. Maxwell [1867] and, more recently, Vick and Ozisik
[1984] modified the constitutive equation by including a relaxation term
that includes some start-up time, d, for the initiation of heat flow.

The modified constitutive equation,

as__ 3.1
q + ¢ at k 3‘1 (2.1.3)

was first suggested by Maxwell to account for finite diffusion velocity
and shows that the spontaneous release of a finite pulse of energy
causes a thermal wave front to move through the body at a finite
velocity. The wave of thermal energy dissipates exponentially as it
moves through the medium, see Figure 2.1. Other researchers (Vernotte
[1958] and Chester [1963]) have used the modified constitutive equation
to account for thermal waves in their experiments with helium at tem-
peratures close to absolute zero.

This thesis will not concern itself with the property of the ther-
mal wave that moves through the medium. When the relaxation time o is
small, equation (2.1.3) reduces to equation (2.1.2) except for extremely

small real times or temperatures near absolute zero.

15

.haooa Howucocodxm Suds o>o3 mmuocm Hmauozh

sun

a.~ ouawwm

 

 

ucoum o>e3
Haauoch mcu>oz

l

 

 

Ammonm
« Heauofi.

 

 

 

 

 

 

Ti
tion a;
equatic
the hon
the tim
techniq
functio
and fort

summariz

2.2 H

2.2.1 5:

equation
Petature
The part1

dimensiOn

16

The solutions to heat conduction problems using the Green's func-
tion approach is presented in Section 2.2. The partial differential
equations for multi-dimensional heat conduction are developed along with
the boundary and initial conditions. In Section 2.3, the concepts of
the time-reversed auxiliary Green's function are presented along with a
technique to determine the Green's function and multiplying Green's
functions to obtain multi-dimensional Green's functions. The structure
and formalism of the approach are discussed in Section 2.4. Section 2.5

summarizes and concludes this chapter.

2.2 Mathematical Derivation of'Heat Diffusion in a Body

2.2.1 d 0 da 0 di 0 s

The purpose of this section is to discuss the partial differential
equation and its boundary and initial conditions used to generate tem-
perature distributions in infinite, semi-infinite and finite bodies.
The partial differential equation that describes the transient, multi-

dimensional, linear heat conduction in cartesian coordinates is,

2
since). :21 16.1:
v 'r + k - 13 arj + 71' -a at (2.2.1)

2
The symbol V is the Laplacian operator defined as,

2
V - L- (2.2.2)

29
6r

 

 

 

the s
is th
assum
and k
solid
fusiv
440] :
nonlit

exampj

and p

17

the symbol g(z,t) is the internal heat generated in the solid body and r
is the rectangular coordinate (i.e. x, y, or z). The boundaries are
assumed to be parallel to the rectangular coordinates. The symbols a
and k are the thermal diffusivity and conductivity, respectively, of the
solid body. The diffusion equation is nonlinear when the thermal dif-
fusivity or p Cp are functions of the temperature. Ozisik [1980, pg.
440] shows that for certain boundary conditions and restrictions, the
nonlinear equation can be transformed into a linear equation. For

example, if the thermal conductivity can be expressed as,

k - 1 + fllT (2.2.3)
and p C1) is expressed as,

p Cp - l + flzT, (2.2.4)

the diffusion equation is nonlinear. If the constants 81 and 82 are
equal, the diffusion equation can be transformed into a linear equation
by the Kirchhoff transformation and the thermal diffusivity becoming a
constant.

The term 8 3%“ represents energy carried by a convective flow in

j

the rj direction and the term 7T could represent generation that is
proportional to the local temperature. The terms 8 and 1 are constant.
In heat transfer analysis, the 7T term may also represent the effect of
side losses for a fin.

The nonhomogeneous boundary conditions that are applied to the

surfaces are written as,

QT 11 _
k1 ani + (pcb)1 at + hiT fi(ri’t)’ (2.2.5)

18

where the integer i represents the surface to which the boundary condi-
tion is applied and n1 is the outward pointing normal to the i-th
surface. The symbol h1 is the convection coefficient and fi(r1,t) is
the nonhomogeneous forcing function associated with the i-th surface.
The k1 symbol will represent the thermal diffusivity of the body, except
when a temperature condition is imposed on the surface, k1 is set equal
to zero.

An additional term, (pcb)1 gg', has been added to the traditional
boundary equation that takes into account a thin film occurring at the
i-th surface. It could represent a thin film of gold on a glass or
silicon substrate or (see Carslaw and Jaeger [1959, pg. 128]) a slab in
contact with a well stirred fluid. A laminar sublayer will form across
a boundary of a slab placed in a flow of fluid. If the flow velocity of
the sublayer can be considered constant, the sublayer can be considered
a thin film.

A thin film term acts to hold or store some energy. It is assumed
that the thickness of the thin film is small enough so that the tempera-
ture at the surface of the film is the same as the temperature at the
surface of the body. The term (pcb)1 represents the storage capacity of
the thin film at the i-th surface.

The initial condition necessary to complete the description for the

solution of equation (2.2.1) is,

T(r.0) - F(r). (2.2.6)

The heat conduction equation (2.2.1) in the rectangular coordinate
system with convection in a direction parallel to the rectangular coor-
_ dinate directions can be reduced to a more desirable form by defining a

new variable (see Ozisik [1980, pg. 75]), W(;,t), where

 

 

 

TiLt)

Substf

VW+

where

rapre

{Epre

aPPli

duCti

COndi

PV‘
24;):

“here

19

51 51 52 fie
T(r.t) - Wir.t) expi-vt] exvlg; (x -'5' t)] exeig; (y - 5‘ t)] -
fie fie
exp[§; (2 -‘E- t)] (2.2.7)

substituting equation (2.2.7) into equation (2.2.1) yields,

2 we; is?!
v w + k - a at (2.2 8)
where,
51 1 52 52
8'(r.t) - g(r.t) expl-vti cuppa; (x - -2- t)] exp[-§; (y - T t)] .
B: s
exp[-§; (z - 5‘ t)] (2.2.9)

represents the heat generation in the solid body.

Equation (2.2.8) is easier to solve than equation (2.2.1), and it
represents a broad class of conduction problems that include convective
diffusion, heat generation and generation proportional to the local
temperature. The transformation given by equation (2.2.7) must be
applied to the boundary and initial condition as well as the heat con-
duction equation. Applying the transformation to the general boundary
condition, equation (2.2.5), yields,

2! as . _ .
k1 dn1+ (pcb)1 at + hiW fi(ri,t), (2.2.9)

where,

2
B. fl
. .41 ..i

20

the new convection coefficient, and,

51 51 B2 52
fi - fi(ri’t) exp[-1t] exp[-§; (r1 - 5' t)] exp[-§; (r2 - §—'t)] .
fl: 5:
exp[-§; (r3 - 3‘ t)] (2.2.11)

represents the new forcing function.
The boundary conditions for a diffusion problem has been stated
previously as,

1T 11 _
k1 ani + (pcb)1 at + hiT fi(ri’t)° (2.2.12)

Five distinct classes of boundary conditions can be obtained from this
equation.
The first kind is when the temperature is prescribed condition at a

boundary, see Figure 2.2a,
T - fi(ri,t). (2.2.13)

Equation (2.2.13) can be obtained from equation (2.2.12) by letting ki =
0, (pcb)1 - 0, and h1 - l. The term fi(ri’t) is the prescribed tempera-
ture history at the boundary. If the temperature at the boundary is

zero, equation (2.2.13) on the boundary becomes the homogeneous boundary

condition,
T - 0. (2.2.14)

The temperature boundary condition is also called a Dirichlet condition.

 

 

T=T
o

21

. Figure 2.2a

Semi-infinite Slab
With a Temperature Boundary
Condition.(X10). .

 

 

 

 

Figure 2.21: ‘3 "' 90 L ---""x

Semi-infinite Slab
With a Heat Flux Boundary
Condition (X20). .1

 

 

 

 

 

Figure 2.2c

Semi-infinite Slab
With a.Convective Boundary
Condition (X30).

 

 

22

The second kind of boundary condition (also called a Neumann

condition) is one in which the heat flow is prescribed, see Figure 2.2b,

fl ..
k1 ani fi(r1,t), (2.2.15)
‘1

where f1(r1,t) is a prescribed heat flux. This equation is obtained
from equation (2.2.12) by letting h1 and (pcb)i equal zero.
If the heat flux at the boundary is zero, the boundary condition

becomes homogeneous and is referred to as being insulated,

31 _ 0. (2.2.16)
ni r

i

A Robin condition or convective boundary condition occurs when
there is a linear combination of the temperature and the normal deriva-
tive of the temperature -- a boundary condition of the third kind; see
Figure 2.2c. The value of the thin film coefficient in equation

(2.2.12) is set to zero and the boundary condition equation becomes,

 

 

ii -

k1 ani + h1 T fi(r1,t) (2.2.17)

r1
or,
L %I + 31 T - f"(r ,c) - L, (2.2.18)

n1 1
'1
hiL

where f"(ri,t) was previously defined, Bi - _k— is the Biot number for

the solid body at the i-th surface and L may be the thickness of the

 

 

 

 

 

body.
thermal
homogen
A .
exists
the sur
thin fi
gradiet

the hon

W
0.10)
ales

 

 

‘ihere
thin f

equal

23

body. The Biot number is the convective coefficient divided by the
thermal conductivity. If f”(r1,t) is zero, the boundary conditions is
homogeneous.

A boundary condition of the fourth kind, or Carslaw condition,
exists when a "thin" film with no convective heat loss is prescribed at
the surface of the thin film's outer boundary; see Figure 2.3a. The
thin film is assumed to have high conductivity therefore no temperature

gradient exists through the film. The value of h is set to zero and

 

 

i
the boundary condition equation (2.2.12) becomes,
21 a: _
k1 an, + (pcb)1 a: fi(r1,t) (2.2.19)
1‘1
or,
b
21' +‘_'£_’_1u_£1‘_21'_‘l (2220)
an k at k ’ ' ‘
i i i -
1"’1

where k1 is the thermal conductivity of the solid body and (pcb)i is the
thin film coefficient. This equation is homogeneous if fi<ri’t) is set
equal to zero.

The boundary condition of the fifth kind, or the Jaeger condition,
includes all of the terms in the boundary condition equation described
by equation (2.2.12), see Figure 2.3b. The term k1 is set to the ther-
mal conductivity of the solid body, the term hi is the convection
coefficient and (pcb)1 is the thin film coefficient. It is homogeneous
if f1(r1,t) is set equal to zero.

Another type of boundary condition that is considered will be

called a boundary condition of the zeroth kind. This boundary condition

24

i 11

ND
q = qo "—"'
1——-—
a

\\X\\\\\ 0'
(T

 

 

 

\\\\\\\\\\

L

 

.(

 

X
II
0
X
It
I"

Figure 2.3a Finite Slab With a Non-convective Thin Film and a
Heat Flux Condition (X42).

i

’2
.\\\\\\§\\Xi v

f

V

 

\\\\\\\\\\

 

 

 

X
II
C
X
r

Figure 2.3b Finite Slab With a Convective Thin Film and a Heat
Flux Condition (X52).

 

 

has no '
infioit

referre

sion e

2.2.2

func-
equa

Peta

the

The

fro

SCH;

25

has no physical boundary, such as in the case of an infinite or semi-
infinite body; see Figure 2.4. The zeroth condition has also been
referred to as the natural condition.

The initial condition for the heat conduction equation,
T(r.0) - F(r). (2.2.21)

This concludes the description of the transient linear heat diffu-

sion equation that will be used in this thesis.

2.2.2 ' o u on

The purpose of this section is to generate the auxiliary Green's
function equation. This equation, when combined with the heat diffusion
equation, gives a formalism to the solution of heat diffusion for tem-
perature distribution.

The development of the auxiliary Green's function equation follows.
the work of Morse and Feshbach [1959], Ozisik [1980], and Beck [1984].
The auxiliary Green's function equation for heat diffusion is obtained
from the heat conduction equation model with an instantaneous heat pulse
source of unit strength. The boundary conditions for this auxiliary
problem are homogeneous and the solution to this auxiliary problem is

called the Green's function, G. The auxiliary Green's function equation

is,

a V26(;,t|;',r) + a; - 105(1- - 7) - 3% (2.2.22)

in the region R and t > r,

26

q=qo

 

 

(a) Semi-infinite Body (X20)

 

l
i
I
l
l

 

3*).

(b) Infinite Body (X00)

Figure 2.4 A Semi-infinite Body and an Infinite Body.

 

 

sub j e

....
H.
3’ rs

and s

of the
at so:
SI’HIbol

reSpe:

27

subject to the boundary conditions,

29. flfi _
k1 ani + (pcb)1 at + hIG 0 (2.2.23)

on the i-th surface,
and subject to the causality relationship,

C(r.t|r'.f) - 0 (2.2.24)

for t < r.

The physical interpretation of the Green's function is the response
of the system to a unit impulse of heat that occurs at some time 1 and
at some position.;'. The solution to the auxiliary problem is given the
symbol G(;,t|;',r) where I and t are the point and time of interest
respectively, 1' and r are the position and time of the impulse.

Equation (2.2.24) is called the causality relationship because it
expresses the relationship between the impulse, which occurs at time 1,
and the effect of the impulse, which can occur only after time 1. This
means that for times t < r or, mathematically, for times -t > -r , there
is no effect and the Green's function is zero.

The heat diffusion equation and the auxiliary Green's function
equation are parabolic in time due to the appearance of the first
derivative with respect to time. This means, for example, that a solu-
tion to the heat conduction equation, T(r,t), is not the same as the
solution T(r,-t). The diffusion equation and the auxiliary Green's
function equation are asymmetric in time; they can distinguish between
the past and future.

The Green's function also satisfies a reciprocity condition,

The fum
time -t
true an

reverse

2
avoc-

Where
auxili

Ufa:

2.3

C011:
plac
aux:

ture

28

G(:.t|r’.f) - G<r'.-rlr.-t). (2.2.25)

The function G(;',-r|;,-t) is the effect of a source at location 3 and
time -t at a point 1' and at time -r. The causality condition is still
true and the reciprocal function G(;',-r|;,-t) satisfies the time

reversed auxiliary Green' function equation,

2
a voc + u; - y) 5(c - r) - - 3%, (2.2.26)
2
where V0 is the Laplacian operator for the ;' coordinates. This
auxiliary Green's function now describes the development of the effect

of a source placed at position 1 and at time -t.

2.3 The Green's Function.Approach
2.3.1 c va 0 0 th G een' unct o A roach

The purpose of this section is to combine the heat diffusion and
the auxiliary Green's function to obtain a formalism for the solution of
heat diffusion problems for the temperature distribution. A general
expression for the solution to the heat conduction equation [0zisik,
1980] can be generated by combining the heat conduction equation and the
time-reversed auxiliary Green's function equation. Multiplying the heat
conduction equation (equation (2.2.8) with r replaced by r' and t re-
placed by r) by the Green's function, multiplying the time-reversed
auxiliary Green's function equation (equation (2.2.26)) by the tempera-

ture and subtracting gives,

 

This

and

Wher

Equa

,1

’
,.
..
(D

the

29

(szT-TV§G)+uL1'z‘L)'G-%8(;-;')6(t-r)T-;'a'('g—Tl.

a 87
(2.3.1)

This equation is integrated over the total region, R, with respect to r'
and r. The term r goes from 0 to t+¢, where e is arbitrarily small value

of time which will be made to approach zero. This yields,

t+e t+e
2 2
I I a ( G VOT - T VoG ) dv' dr + I I’% G g(;',r) dv' dr
r-O R 1-0 R
t+e
- T(r,t) - I’[ G T ] dv'. (2.3.2)
R 7-0

where dv' is an volume element in the region R. Rearranging the above

equation for the temperature distribution gives,

t+e t+e
T(r,t) - - I [ G T ] dv' + I Iii G g(;',r) dv' dr
R r-O r-O R
t+e
2 2
+ I I a ( G VOT - T VOG ) dv' d1. (2.3.3)
r-O R

The only term on the left hand side is the temperature distribution of
the body at location 1 and at time t. This equation is now examined
term by term.

The first term on the right hand side of equation (2.3.3) is,

t+e

-Jl[c'r] dv'. (2.3.4)

R r-O

 

EvaluA‘

functi
before

is the

and i

tribu

This

S\

30

t-i-e
Evaluating - [c 'r] - - (o - c T|f_o) - + c TIM). The Green's
r-O

function evaluated at r-t+e is zero because the effect cannot begin
before the impulse. When 1 - 0, the temperature distribution, T(;,0),

is the initial temperature distribution, F(;). This term becomes,

Il-IMLdrm)mr>wu (21»
R

and is the effect of the initial distribution on the temperature dis-

tribution and is designated 11.

The second term on the right hand side of equation (2.3.3) is,

t+¢

I3 - I I § G(1,tl;',r) g(;',r) dv' dr. (2.3.6)
r-O R

This term is the effect of a distributed heat source, g(;,t), on the

temperature distribution and is designated I3.

The third term on the right hand side of equation (2.3.3) is,

t+e
2 2
r-O R

Green's theorem can be used to change the volume integral to a

surface integral so that,

t+e

2 2
I I a ( G VOT - T VOG ) dv' dr
r-O R

 

 

where ‘

the boi

the nu

bounda:
Green':
(2.2.11

temper,

S’If‘a’

Integra

Yields

31

 

t+e s

.. fl _ 2.9.
I E Ila ( C an, T a“: ) dSi dr, (2.3.8)
r-O i-l S1 r -ri r -r1

.0.

where the term a
“1

the boundary surface 81’ see Figure 2.5, where i-l,2,3,...,s, and s is

is differentiation along an outward drawn normal to

the number of boundaries.

The integrand of this integral can be expressed in terms of the
boundary conditions of the heat conduction equation and the auxiliary
Green's function equation. Multiplying the boundary condition equation
(2.2.12) by the Green's function, multiplying equation (2.2.23) by the

temperature and subtracting yields,

 

G21 _ 25L _5211‘3 fl 2c:_ 21
(an Tan ) k 6* k (Taf +667”
1 r'-r 1 r'-r 1 i
1 1
or, (2.3.9)
(Gm: 429. )_5£1'_“’G+ (“human
an1 ani R1 R1 67
r'-r1 r'-ri

(2.3.10)

Integrating over the surface S and the time, 1, for constant thermal

i
diffusivity and boundary conditions of the second through fifth kinds,

yields,

t+e s

11
I a (G an

r-O i-l S.
1

wrfi )dS.dr
an 1

I i I
r -r. r -r.
1 1

 

i

Pi

32

 

aiii-hr

flI.+ h I
8x \ ’////__-ax

 

 

/\_/
P.

Figure 2.5 A Description of Normals at the Surface of a Finite
Body.

 

 

Th.
effect
will be
aquatic
thin f;
design

F

15 2e:

33

t+e s f (r' 7)
- a I E I '1'E1——' G(;,t|;',r) dSi dr
r-O i-l Si 1
s (p c b)
+ a} I "T—i c<;,t|;',0) F(r') dSi. (2.3.11)
1
1-1 31

The first term on the right hand side of equation (2.3.11) is the
effect of the boundary conditions on the temperature distribution and
will be designated as 14' The second term on the right hand side of
equation (2.3.11) is the effect of the thermal storage capacity of the
thin film on the temperature distribution of the solid body and is
designated 12.

For a boundary condition of the first kind, since G at the boundary

is zero, equation (2.3.8) becomes,

t+e s'
- - . 2L
I4 0 I E I fj(rj,r) an ‘ de dr. (2.3.12)
r-O j-l Sj j r'-rj

where s' is the number of boundary conditions of the first kinds.
Drawing together the four terms yields an expression for the tem-

perature distribution as,

T(r,t) - I1 + 12 + 13 + 14

- I C(r.tlr'.0) F(r') dv'

R
s
(pcb)
+ a E I -——E-——1 G(r,t|r' 0) F(r') d8
1

 

bc

*9)

Ci

(1)

rt

34

t+€
+ I [fiandrn>yyn>w'm
r-O R

t+e s

f (r',r)
+ a I’ E I -1—E1——- G(;,t|;',r) dsi d1
r-O 1-1 s 1

i
(for boundary conditions of the second through fifth kinds)

t+e s'
. fig.
- a E I fj(rj,r) anj de (11'. (2.3.13)
r-O j-l SJ r'-rj

(for boundary condition of the first kind only)

The four terms on the right hand side of equation (2.3.13) describe
a formalism to be used to solve for the temperature distribution in a
body for boundary conditions of the second through fifth kind. The
first two terms represent the effect on the temperature distribution
caused by a nonzero initial condition. The third term represents the
effect caused by a volume energy source and the last two terms repre-
sents the effect caused by nonhomogeneous boundary conditions.

If a boundary condition of the zeroth kind occurs at a surface i=j,
the last two terms in equation (2.3.13) are omitted for that surface.

The temperature distribution in the one dimensional cartesian

coordinate x for a one dimensional slab is,

L
T(x,t) - I G(x,t|x',0) F(x') dx'

x'-0

2
p c b)
+ a E __—k__—i [G(x,t|x',0) F(x')]

(
1-1 1 x"x1

35

t L
+ I I § G(x,t|x',r) g(x',r) dx' dr
r-O x'-0
s f (x'
a! E -1—-1——- C(x, tIx' ,r)x d7 .
r-O i-l

(for boundary conditions of the second through fifth kinds)

t s'
- a {I §f1(xj,,r) 39- dr, (2.3.14)
r-O 3-1anj 3 -xj

(for boundary condition of the first kind only)
where L is the length of the slab in the x direction.
The temperature distribution in the two dimensional cartesian

coordinates x and y for a two dimensional plate is,

T(X.y.t) - I I C(X.y.tIX'.y'.0) F(x'.y') dx' dy' +

 

x' y'
(p c kb)
6.3—i [C(X'thlx':}".0) F(x',y')] +
i-l x'-xi
Y 'Yi
t
I J. I§G(X,y,tlx',y1’7) g(X',Y',T) dA d? +
7'0 x' y'
t s fi(Xi,yi,r)
a I E ki C(x’y’tlx"y"')x'-xid' ,
1.0 1-1 Y'-yi

(for boundary conditions of the second through fifth kinds)

I

t S
- a I E ”1(x ,y. 7) %§— - - dr, (2.3.15)
,_0 1-1 J X' xj y' yj

(for boundary condition of the first kind only)

36

23.2 W

The objective of this section is to discuss some methods to deter-
mine the correct Green's function based on the boundary conditions.
There are many procedures for the determination of the Green's function.
One procedure, which will be discussed in the next chapter, is by the
method of Laplace transforms.. A second method, which produces results
similar to the Laplace transform method, is the method of images. Both
of these techniques produce Green's functions that are computationally
efficient at small dimensionless time. The third method, which produces
Green's functions that are computationally efficient at large dimension-
less times, is obtained by using the traditional method of separation of
variables.

The method of images requires the temperature in a finite or semi-
infinite body caused by a supply of heat at certain points, and the
removal of heat at other points. A supply of heat at a point is called
a source while the removal of heat at a point is called a sink. In a
one dimensional infinite body, the temperature distribution due to an

instantaneous heat pulse is,

T(x.x'lc.r) - Q [a a a (.-.)1-1/2 expt-(x-x'>2/<a «2 (t-rm.

(2.3.16)

where Q is the strength of the source, x' and r are the location and
starting time of the source, and x and t are the location and time of
interest. Multiplying this temperature distribution by pc and integrat-

ing over the entire infinite body yields,

37

I ch dx' - Qpc I [a a a (c-r)]’1/2 exp[-(x-x')2/(4 a (t-r))] dx'

- Qpc. (2.3.17)

which shows the total amount of heat liberated in the infinite solid is
Qpc.

When the infinite body is bounded by planes, the bounding planes
can be considered to act as heat mirrors which reflect the sources and
sinks. The Green's function of a bounded body is simply the distribu-
tion of the original source plus the effects of these reflections, see
Figure 2.6. The method in Chapter 3 will give a technique for the
placement of the sources, sinks and correction terms for the positions
shown in Figure 2.6.

In his book on partial differential equations, Sommerfeld [1967]
describes the method of images for heat diffusion in a slab as a room
with parallel mirrors. This simulates a finite body of length L. A
light placed in the room will be reflected by both mirrors not once, but
in infinite repetition. The reflections of the light source are the
heat poles with a period 2nL where L is the distance between the mirrors
and n is an integer index. The summation of the poles yield the Green's
function for small dimensionless time. The reflections can be con-
Sidered to form background correction factors to the effect of the
source.

The method of images is restricted to systems that have boundaries
comPosedof straight lines. Other shapes may be considered but the
method may be applied only approximately.

A second approach that is frequently used to solve heat conduction
Prohlems is by using the separation of variables. This technique gives

Green's functions that are efficient for large times. The homogeneous

38

////////////// “L
T "T

Ix-x'i

l .

T" _,
. //.{,. c 5///////// "

 

 

 

.(2L+x+x')

F1Sure 2.6 Reflections of Sources and Sinks in a Finite Body.

a heat C
, technf

tion,

Since

and

C
L U]

01).

A!

39

heat conduction equation is solved using the separation of variables
technique, with homogeneous boundary conditions and an initial condi-

tion, F(z), and the solution is expressed in the form,

T(X3t) - J. 1((191' 3:) F(I') dV' 0 (203-18)
R

Since the Green's function solution for this case is,

T(I.t) - I G(I.I'It.0) F(I') dv'. (2.3.19)
R

and it follows that,

G(r.r'|t.0) - K(r.r'.t). (2.3.20)
Replacing t with t-r in K(;,;',t) gives the general Green's function for
solution to nonhomogeneous problems, G(;,;'|t,r).

It is important to note that to solve the nonhomogeneous problem
for temperature, the homogeneous problem for the Green's function is all

that is needed to be considered.

2.3.3 zzgggggg of Green'; Eugggigns

The purpose of this section is to show the combinations of Green's
functions that are mathematically possible for rectangular coordinates.

.A.unique feature of the Green's function approach to the solution
°f transient heat diffusion problems is the capability of multiplying

one dimensional Green's functions to obtain Green's functions for two

and th}

 

the zei

the te

C(m.

T0 pr:
into -
and t

added

auxil

Ihe

40

and three dimensions in cartesian coordinates for boundary conditions of
the zeroth through third kinds. The development of this idea follows
the techniques of Ozisik [1980] and Beck and Yen[l984b].

In cartesian coordinates it is desired to prove,

C(X.y.2.tIX'.y'.z'.r) -

Gl(x,tlx',r) - G2(y,t|y',r) - G3(z,tlz',r). (2.3.21)

To prove this, the right hand side of equation (2.3.21) is substituted

into the auxiliary Green's function equation and the boundary conditions

and then compared to the equation of the individual coordinate equations
added together.

Substituting the right hand side of equation (2.3.21) into the

auxiliary Green's function equation (equation (2.2.26)) yields,

2 2 2
a c, a c, a G3
a G29: 2 + 0103 2 + C162 2

3x ay 82
361 6G2 6G3
- G26a 5;— + 6163 -3:_ + G162 -§:- -
6(x-x') 6(y-y') 6(z-z') 6(t-r). (2.3.22)

The one dimensional equations that describe G1, 62, and G3 are,

 

 

2
a 61 as
l 1 .
<1 [ 6x2 - a a, ] - 6(x-x ) 6(t-r) (2.3.23)

 

 

2
2 5(y y ) 6(t 1) (2.3.24)
By a 61

 

 

Multi

and a

inf .

41

 

 

2
6 G, 6G,
a [ az, - i af ] - 6(2-2') 6(t-r) (2.3.25)

Multiplying these three equations by G2G3, GIGS, and GIG2 respectively,

and adding gives,

2 2 2
a G, a G2 a G3
a G26: 2 + 6163 2 + C162 2
8x ay 62
ac, 3G2 8G3
- GzG, '5:' + GIGa -5; + G1G2 ‘3; -

c,c3 8(x-x') 6(t-r) + G163 6(y-y') 6(t-r) + 0,0, 6(z-z') 6(t-r).

(2.3.26)

The terms on the left hand side of equations (2.3.22) and (2.3.26) match
therefore it must be shown that the right hand terms of these equations

are equal. This means, for example, that,
0263 6(x-x') 8(t-r) - % 6(x-x') 6(y-y') 5(2-2') 6(t-r) (2.3.27)

Integrating equation (2.3.24) with respect to r from minus to plus

infinity ,

 

 

m 2 m
a G2 ac2
In [ 3y: - i 81 ] dr - 6(y-y') J. 6(t-r) d1 , (2.3.28)

Which yields ,

92 - 6(y-y') H(t-r). (2.3.29)

 

When

Tvfhen

func

Su‘r

42

When tflf, the left and right hand side of equation (2.3.27) are zero.
When t-r, G2 and G, act like unit step functions, H(t-r), multiplied by

functions independent on time. Therefore,
2
62G, - 6(y-y') 6(z-z') H (t-r). (2.3.30)

Substituting equation (2.3.30) in equation (2.3.27), integrating over

time yields,
I H (t-r) 6(t-r) dt - I 3 6(t-r) dt
11’ 1
_ 3 - 3 _ (2.3.31)
t-f

 

The same procedure is used for the remaining terms on the right hand

side of equation (2.3.27) and results in,
3['% 6(x-x') 6(y-y') 6(z-z') 6(t-r) } -
6(x-x') 6(y-y') 6(z-z') 6(t-r), (2.3.32)

which completes the proof for the differential equation.
The boundary conditions must also be satisfied for the product
relationship to hold. The general boundary condition for boundary

conditions of the first, second and third kind is,

LG. _
ki ani + h1 G 0. (2.3.33)

43

Substituting equation (2.3.21) in this expression gives,

ac1

or, dividing through by G2G3 gives,

ac1

k1 5;— + h1 G1 - 0, (2.3.35)
i

which is the same expression as equation (2.3.33). This proves that the
product relationship will hold for boundary conditions of the first,
second and third kind. 3

The general equation for boundary conditions of the fourth and

fifth kind is,

as. as _
k1 ani + (pcb)1 at + h1 c 0. (2.3.36)

Substituting equation (2.3.21) in this equation gives,

6G1
62G, 32:— i G1G2G3h1
8G1 6G2 6G3
i (pcb)1 G263 _az + GIG3 '3; + GIG2 '5? - 0 (2.3.37)

which cannot be reduced to the form of equation (2.3.36). Consequently,
the product relationship of the one dimensional Green's functions is not
valid for two or three dimensional diffusion problems when there are

boundary conditions of the fourth or fifth kinds.

44

2.4 Formalism of the Green's Emotion Approach

The description of the temperature distribution based on the
Green's function approach, equation (2.3.13), provides a logical and
mathematical structure to the solution of heat diffusion problems. This
formalism is used to begin a structured data base of solutions based
solely on the geometry and the boundary condition types of the specified
problems.

The dictionary defines formalism as a rigorous adherence to recog-
nized forms. A mathematical formalism uses some basic "building blocks"
or functions in a structure to obtain analytical solutions. In the
Green's function approach, the basic "building blocks" are contained in
the data bases. Both data bases in the Green's function approach, the
Green's function data base and the integral data base, are finite and
analytical. This will restrict the capabilities of the method to a
finite number of cases, but the solutions will be analytical. The data
bases do no need to be analytical in general. For example, in the
integral data base, an integral that does not have a closed form could
be expressed in numerical form, but the solution based on the numerical
form will not be analytical.

The Green's function approach to heat diffusion problems relies on
the availability of Green's functions. A single Green's function is used
in each term on the right hand side of equation (2.3.13). Independent
of the initial condition and volume energy source, the Green's function
is determined from the boundary conditions. This suggests that a look-
up table or data base of Green's functions is important to the
formalism. Furthermore, a complete table of one dimensional Green's

functions for the cartesian coordinate system is available.

 

si

bc

In

4S

Beck and Litkouhi [1984] have suggested a numbering system for
single and multi-dimensional Greens functions based on the types of
boundary conditions that occur on the surfaces of the body. Figure 2.7
shows the distinct number of cases for the Green's function for the one
dimensional rectangular coordinate system. For example, a one dimen-
sional slab with a boundary condition of the first kind at the left
surface and insulated on the right surface is designated as X12. The
numbering scheme allows the Green's functions to be manually and com-
puter catalogued and the effort necessary for locating the functions
will be the effort necessary to establish the catalogue number for the
specified problem. Figure 2.7 shows the distinct number of cases for
the Green's function for the rectangular coordinate system.

The numbering system proposed is important since the temperature
distributions obtained by the Green's function approach can be
catalogued and stored and need not be re-evaluated. The numbering
system uses the types of boundary conditions and determines a Green's
function that represents a plane, line, or point source that occurs at
some point x' in the slab (or on the boundary), and at some time r. Due
to the linearity of the problems, this function is then multiplied by
the forcing function and integrated over the boundary, area, and/or time
of interest.

Since all of the Green's functions for the six linear boundary
conditions have been catalogued, it is easy to calculate new solutions
by implementing the Green's function formalism. The convergence of the
new solutions is greatly improved over the solutions based solely on the
method of separation of variables or the method of Laplace transforms.

For small times, the Laplace transform method is most efficient.

For dimensionless times greater than 0.05, more terms are necessary in

46

maoanoum HmCOchwEHo mco pom coauocsm m.:mmuo mo mmmmo uocfiumfin n.N muawwm

>mHm20mo WCz_LIII“N >mhufiouo WtzhuzT_§um

 

 

  
        

max sex mmx max sex
sex mex max esx
mmx smx
mmx smx
eex

mh<km >Q<whm Oz

 

cm“
3%
Gmx
omx
2x

 

 

 

Dex

 

 

E5280 u:z_....z_‘\

 

 

the La
ing to
the nu

tract:

dimen:
solut
probl
kinds

of th

due 1
will
ary .
vari
vOlu

With

in t
tier

redL

type

47

the Laplace transform series for convergence. At this point, by switch-
ing to the separation of variables function, convergence is retained and
the number of terms evaluated in the series remains small or at least
tractable.

The ability to use the product of Green's functions for multi-
dimensional temperature distributions is advantageous since the same
solution kernel may be programmed for one, two, or three dimensional
problems for boundary conditions of the zeroth, first, second, and third
kinds. The Green's function approach does not allow boundary conditions
of the fourth and fifth kinds for problems in two and three dimensions.

The great usefulness of the Green's function procedure is startling
due to the vast number of cases that are involved. The CANSS program
will treat heat conduction problems that deal with nonhomogeneous bound-
ary conditions of the zeroth through fifth kind, time and space
variation in the boundary conditions and initial condition, constant
volume energy sources, terms associated with fins, and terms associated
with flow.

The symbolic formalism of the Green's function approach permits the
Lase of either the small time Green's function or the large time function

LJn the solution to the stated problem. Time partitioning of the solu-
tion, as suggested by Beck and Keltner [1985] may be used effectively to
reduce the computational load.

Time partitioning for the solution of linear, transient diffusion
type equations uses a linear combination of two separate, but equivalent
EKpressions. One expression is obtained using a Laplace transform
technique and converges quickly when time is small; the other expression
can be obtained using the method of separation of variables and will
converge rapidly when the dimensionless time is large. Therefore, two

expressions are available for use, each having a region of dimensionless

48

time for which it is best suited. The Green's function expressions can

be used to generate a solution that converges for any dimensionless time

by combining both the small and large dimensionless time expressions .
An example of the combined Green's function for a one dimensional

slab insulated on both boundaries is,

G(xL,xi,t*) -

i [4 1r 314/2 E (exp[-(2m+xL-x£)2/4t*] + exp[-(2m+xL+x£)2/4t*])

111—um
(2.4.1)
* *
for t small, and for large t is,
*
C(xL,x£,t ) - Constant
O
+ l 2 Eexp[ n21r2t*] cos[n ] cos[nzx'] (2 4 2)
L J; «XL L . .
n-l
where xL and x1" are normalized with respect to L and t* - 91%)- . Note
L

that the exponential term converges rapidly to zero for small dimension-
less times in the small time expression and rapidly to zero for large

dimensionless times in the large time expression.

2.5 Su-ary

A formalism for solving simple, multi-dimensional, linear, tran-

Sient heat diffusion problems is presented in this chapter. This

49

formalism consists of a multi-dimensional equation and a group of kernal
functions called Green's functions. The formalism allows the applica-
tion of six distinct types of boundary condition, and the facility for
handling nonhomogeneous initial conditions and volume energy sources.

The unique structure of the Green's function formalism and the
availability of the Green's functions themselves are ideal candidates
for symbolic solution. A simple transformation allows the formalism to
handle terms that include convective diffusion and temperature gener-
ation proportional to the local temperature.

Two methods are discussed to determine the Green's functions. A
third method will be described in the following chapter. These methods
provide a data-base form for cataloging the Green's functions.
Multiplying one dimensional Green's functions for certain combinations

of boundary conditions is shown.

Chapter 3

SHALL rm “'8 FUNCTIONS os'ranum USING IAPLACE TRANSFORMS

3 . 1 Introduction

Carslaw and Jaeger [1959] and Horse and Feshbach [1953]
demonstrate the use of the small time Green's functions, and Ozisik
[1980] summarizes the importance of these functions in their classical
heat transfer and physics texts. The small time functions permits
efficient investigation of transient thermal activity at very small
times. For small times, the heat conducting bodies are thermally semi-
infinite; that is, temperature changes are contained only in the region
of the body near the heated surface. The small time Green's functions
are particularly important in problems where the forcing functions
resemble instantaneous sources, such as in the areas of robotics,

e lectronics, and measuring energy deposition by pulsed lasers.

A limited group of small time Green's functions can be found in
the previously cited references, but no care was taken to organize these
functions. This chapter extends the work of Beck[l984] , Ozisik[1980],
and Horse and Feshbach[ 1954] by (a) organizing and systematizing the
small time solutions for boundary conditions of the first, second, and
third kind, and (b) generating additional small time functions for
fourth and fifth types of boundary conditions.

The boundary condition of the fourth kind involves a surface

film of finite heat capacity and the fifth kind involves a surface film

50

 

 

but

film

per:
man
of
of
dii
mo1
ha‘
pr
to

DL‘

r-i

(“I

51

but also has convective heat transfer at the outside surface of the
film.

The method for obtaining small time Green's functions incor-
porates a symbolic manipulation program called SMP [1983] for the
manipulation and evaluation of large algebraic expressions. The ability
of SMP to perform complicated algebraic manipulations and the necessity
of dealing with advanced mathematical constructs such as integration,
differentiation, factoring, and expansions, to name a few, is the prime
motivation for using SMP. The procedures written in the SMP language
have a generality not found in the traditional scientific numerical
programs (i.e. BASIC, FORTRAN, PASCAL) due to the ability of the program
to manipulate symbolic expressions, as well as the capability for
numerical evaluation.

Table 3.1 and Figure 3.1 show a comparison of the Green's func-
tion for a semi-infinite body versus a finite body at small times. The
location of the source is at coordinate x', the location of the point of
interest is the coordinate x, and the thickness of the finite body is L.
The symbol t is the time of interest and the symbol r is the time when
time source begins. Each body has a heat flux condition at the left

b<>undary and the finite bodies are insulated on the right side. The
diffusivity of the medium is constant and equal to one.

Figure 3.1 shows the semi-infinite body solution is a good
approximation to the finite body solution when the time is small and the
thickness of the body is large. It is not a good approximation when the
time becomes large or the thickness of the body is small. As the thick-
ness of the finite body increases, the Green's function for the finite
bOdies converge to the solution of the semi-infinite body, as expected.
The Green's function for the larger width bodies can be approximated at

Small times by the semi-infinite Green's function. The dependence of

52

.oBHH mamuo> hpom ouacam a mo

 

 

 

mnumceq msowum> How macauocah m.:oouu mafia Hausa H.m ouswwm
U
0N.0 0 v.0 N _..0 00.0 #00 00.0
L p h p b b L P h lb h P b b b h L l- ]. O.—.
j '
CNX .
a ..N P
T / I
r. T.v._.
.. T
.. um...
III00.H n A NNX
.. a
ll [moP
mh.0 l A NNN
W ....0.N
l 0m.0 I A NNX I
b D 1.] - D I P - I b b] {-1 L D P - r 1.] P

 

 

 

00000000000000000000

53

Table 3.1 Small Time for Finite and Infinite Body
Source at x' -

X20

L - o

HHHHHHHHHHHHHHHHHHHH

time

.001
.005
.010
.025
.050
.100
.250
.500
.000
.000

mr-‘00000000

.18261
.82649
.93495
.90875
.84596
.77522
.70583
.64080
.58090
.52604
.47584
.42983
.38757
.34862
.31262
.27924
.24820
.21924
.19216
.16676

L - 0.5

bananas:bananas:haeaeaeaharehfihahdrdrdhi

0 .
X22

.18261
.83002
.96495
.99278
.99851
.99969
.99994
.99999
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000
.00000

P‘P‘P‘h‘h‘h‘h‘h‘P‘P‘P‘P‘P‘F‘F‘P‘P‘P‘F‘h‘

X22

L - 0.75

.18261
.82649
.93496
.90891
.84698
.77865
.71391
.65606
.60579
.56278
.52626
.49542
.46943
.44758
.42922
.41380
.40086
.38999
.38088
.37323

Table 3.2 Exponential as a Function of Time

0000000000

LT expansion

e.1/4c

.00000
.00000
.00000
.00005
.00674
.08209
.36788
.60653
.77880
.95123

0000000000

e-4/4t

.00000
.00000
.00000
.00000
.00000
.00005
.01832
.13534
.36788
.81873

0000000000

SOV expansion

0000000000

-4n t
e

.96129
.82087
.67383
.37271
.13891
.01930
.00005
.00000
.00000
.00000

Point of Interest at x - 0.25

P‘F‘F‘k‘h‘h‘h‘h‘P‘P‘P‘P‘F‘P‘F‘P‘h‘h‘h‘h‘

X22

L - 1.0

.18261
.82649
.93495
.90875
.84596
.77523
.70587
.64094
.58129
.52689
.47747
.43264
.39199
.35516
.32178
.29154
.26414
.23932
.21683
.19645

54

the thickness scale may be removed from consideration of the solution by

*
defining a non-dimensional time parameter, t , where,

t* 9Ԥ=11 . . (3.1.1)

Small and large time solutions can be combined to obtain solu-
tions that are computationally efficient for the total time region. The
solutions are split into time partitions for which the resulting solu-
tions may be evaluated. At early times, the transient solutions are
more efficiently represented by a summation of exponential functions
derived from the Laplace transform (LT) or the method of images. The
exponential terms generated by the Laplace transform technique are
functions of -Cn/t* where CIn is a function of m2 and increases in value
as the index m increases, and t* is the dimensionless time. As the
dimensionless time becomes small, the exponential terms rapidly ap-
proaches zero. For example, the first two typical exponential terms of
the summation derived by the Laplace transform technique for a one

dimensional slab, insulated on both boundaries, are shown in Table 3.2

-<1/4c*> and e'<“/4t*>.

under the LT heading, namely e Notice that for

c*< 0.05, the term e'(1/4‘*) is less than 0.0068 and for c*< 0.01 there
is no contribution out to the fifth decimal place. This means that as
the dimensionless time gets small, only the first few terms need be kept
for an accurate approximation. The additional terms may be dropped

2
without diminishing the accuracy of the computation. When m-2 (m =4),

4 a *
one additional term, e-( / t ), is also shown in Table 3.2. Five
decimal place accuracy is obtained for the sum of these two terms by

retaining only the first term when the dimensionless time is less than

0.1.

55

The solutions derived using separation of xariables (SOV) are
more efficient at large times. Like the solutions found in the Laplace
transform technique, the solutions for the SOV method are functions of a
summation of exponential terms, but the exponential functions are typi-
cally functions of -Cnt* where t* is the dimensionless time. The value
Cn is a function of n2 and it increases in.value as the value of the
index n increases. As dimensionless time increases, these exponential
terms decrease rapidly towards zero. The first two typical exponential

2 * 2 *

terms, e-” t and e'A' t

, of the summation found by using the SOV
method for a one dimensional slab with insulated boundary conditions,
are shown in Table 3.2 under the SOV heading. When t* > 0.500, the
2*

contribution of the first summation term, e-“ t , is less than 0.0072.

For t* greater than 0.250, five decimal place accuracy is obtained when
only the first term in the summation is retained.

Accuracy can be increased for both the LT and the SOV method at
any dimensionless time by including additional terms to their respec-
tive series. Inclusion of the additional terms not only increases the
accuracy of the solution, but unfortunately also increases the amount of
computation time necessary to arrive at the solution due to the time of
evaluation of the additional terms.

Table 3.2 shows for any dimensionless time, the second term in
the LT or SOV series is always numerically equal to or less than the
first term in the series. This means the first term will dominate the
solution if the appropriate series is chosen based on the time.
Retaining an infinite number of terms in either the SOV or LT methods
will yield an exact solution but the choice of series that results in
the most efficient method is dependent of the time.

The analysis methods used in heat transfer are typically the

iinite element (FE) method or the finite difference (FD) method and are

56

numerical in nature. These numerical analysis methods generate informa-
tion at the locations designated by the resolution of both the spatial
grid size and the time step size. A smoothing function is used to
develop information at locations that differ from the spatial grid
points or at times that differ from the time step size. Very fine
temporal and spatial resolution are needed for solution convergence at
early times for FE and FD methods; therefore, at early times these
solution techniques are not efficient in the use of computer time.

The procedure for the small time Green's functions is useful in
the continuing research to develop the unsteady surface element (USE)
method [Keltner and Beck,l981], which is a method for solving linear
transient heat conduction problems. Solutions to certain basic tran-
sient heat conduction problems, called influence or kernel functions,
are used as building blocks to solve problems of complex geometry and
problems that deal with nonlinear boundary conditions in the USE method.
The small and large time Green's functions, along with a combination of
the small and large time Green's functions, can be used to generate the
influence functions in the USE method.

Section 3.2 of this chapter gives a mathematical statement of
the general linear transient heat conduction problem and the development
of the small time Green's funct terms of a one dimensional slab.

The objective of Section 3.3 is to examine some examples of Green's
functions for semi-infinite bodies. Section 3.4 gives examples of how
the small time Green's functions are developed for finite one dimen-

sional bodies. Section 3.5 is a summary of the chapter.

 

3.2

gene
bOUI
uses
hour
the
the

the

Com

tior

57

3 . 2 lathe-atical Development of the Small Time Green's thction

The objective of this section is to develop a procedure to
generate small time Green's functions for one dimensional bodies with
boundary conditions of the zeroth through fifth kinds. The technique
uses the partial differential equation of heat diffusion, the associated
boundary conditions and Laplace transforms to generate an expression for
the Green's function that is accurate and efficient at small times. As
the dimensionless time approaches zero, the approximate solution goes to
the exact solution.

The partial differential equation for linear transient heat
conduction developed in Chapter 2 is, after the appropriate transforma-

tions,

V T + k - a at

 

in region R. (3.2.1)

The associated boundary conditions are,

31 31 _ '
ki an, + (p c b) at + hiT f1(ri,t), (3.2.2)
for'i-1,2,...,s, where the symbol 3 is the number of boundaries, and the

initial condition is,
"IQ-2.0) - F(r). (3.2.3)
The thermal conductivity, k, and the thermal diffusivity, a, are

2
Constant with position, time, and temperature, (;,t,T); V is the

Laplacian operator and n is an outward pointing normal.

be

58

For a one dimensional slab without heat generation, i.e. g(;,t)

- 0, equation (3.2.1) reduces to,

2
L} -131 0 <x <L (3.2.4)
3x a at

The one dimensional transient heat conduction equation (equation

3.2.4) is satisfied by,
2
T - [4 « c:(t:-‘r)].1/2 e-(x-x') /(4 a(t-r)) (3.2.5)

which tends to zero when r a t at all points except x', where it goes to
infinity. This solution is called the temperature due to an instan-
taneous plane source through 2! and at time 1 and of strength unity per
unit area. It is the fundamental transient solution in a planar heat-
conducting body and is actually a Green's function. Since the
describing equation and boundary and initial conditions are linear, this
solution must be included in any solution to a planar region.

A solution, T(x,t), to the transient heat conduction equation
(equation 3.2.4) is required that goes to infinity at x - x' when r a
t, but is zero for every other value of x in O < x < L when r e t and
will satisfy the boundary conditions. The method of obtaining this
solution is similar to that given in Carslaw and Jaeger[1959, pp. 359-
360].

Let the solution for the instantaneous source at r - 0 and at x'

be equal to,

2
u - [4 11' a t1'1/2 e'(""") “4 a 9. (3.2.6)

The

fer

Sim

is:

m 1:?»
x” .

The

& N ,3.»

and

Vhere

Condit

‘here D
and sin

are dEf:

59

The complete solution for the temperature, T, that satisfies the dif-

ferential equation and the boundary conditions is,
T - u + w (3.2.7)

Since T and u are described by the transient heat conduction equation, w

is also; hence w is given by the solution of a subsidiary equation,

2
2:}..iglé 0<x<L,t>0. (3.2.8)

The Laplace transformation of the subsidiary equation for w is,

- 2-
g—¥ - q w - 0 0 < x < L (3.2.9)
dx

2
for q - p/a (3.2.10)
0
and G - J e-pt w(x,t) dt (3.2.11)
t-0

Where p is the Laplace transform parameter.
A solution for w that is convenient for satisfying boundary

conditions as x -+ 0 and x _, L is of the form,
‘0 - D1 sinh(q x) + D2 cosh(q x) (3.2.12)
where D1 and D2 are constants determined from the boundary conditions

and sinh and cosh are the hyperbolic sine and cosine functions, which

are defined as,

 

6O

sinh(qL) - % [ eqL - e‘qL ] (3.2.13)
and

cosh(qL) -'% [ eqL + e.qL ]. (3.2.14)
The Laplace transform solution for T is,

i - a + e (3.2.15)

-q|x-X'|
— a 52—3-L_- + D1 sinh(q x) + D2 cosh(q x)] (3.2.16)

subject to the boundary conditions,

521 + fiT - 0 when x - 0 or x - L (3.2.17)
where £1 - 0 for the Neumann condition

- hi/k for
2
- (p c b)1q /(p c) for
2
- hi/k + (p c b)iq /(p C) for

Setting k

(second kind)

the Robin condition
(third kind)

the Carslaw condition
(fourth kind)

the Jaeger condition

(fifth kind).

1 equal to zero will result in a boundary condition of

t:he first kind (Dirichlet) and causes £1 to go to infinity.

Substituting the solution T into these boundary conditions to

find the constants D1 and D2 for all the possible cases is a tedious and

61

error-prone process when done by hand. The process is made more effi-
cient when symbolic manipulation is used. In SMP, the differentiation
operator D[$expr,{$x,$n,$pt}] forms the partial derivative of expression
$expr successively $n times with respect to coordinate $x, evaluating
the final result symbolically at the point $pt.

Differentiation of T with respect to x (for a one-dimensional
slab) is executed by using the differentiation operator D[T,{x,l,pt}]
where the location (pt) is set to zero or the thickness L. Applying the
two boundary conditions to the solution T results in two expressions and
two unknowns, D1 and 0,, as functions of sinh(qL) and cosh(qL).
Expressing the hyperbolic functions in terms of negative exponentials
and expanding the result in a series by the binomial theorem yields the
solution T as a summation of negative exponential terms.

The solution T is inverted term by term either by the Laplace
transform inversion theorem or, more simply, by a table of Laplace
transforms. Some important Laplace transforms that occur when using
this method can be found in Appendix A. Typically, a table of Laplace
transforms is all that is necessary for the inversion of T when only a
few of the terms in the series are retained.

As an example of the method, consider the X42 case, a one dimen-
sional slab with no heat generation, a nonconvective thin film boundary
condition at x - 0 and insulation at x - L. The transient heat conduc-
tion equation for this case is described in equation (3.2.4). The left
'boundary is described by equation (3.2.2) with h1 and f1(t) set to zero.
Trhe right boundary condition is described by equation (3.2.2) with h2,
(pcb)2, and f2(t) set to zero. Substituting the solution, equation

(3.2.16) for T into the left boundary condition yields,

-qx
(1 - C1qL)§q—L_ + n, - cqu 132 - 0 (3.2.18)

62

substituting the solution into the right boundary condition yields,

E-q(L-x')

D1 cosh(qL) + D2 sinh(qL) - qu - 0 (3.2.19)

Using matrix notation gives,

M H D N
11 12 1 - _1_ 1 (3.2.20)
M21 M22 D2 29L N2

where M11- 1, M12- -Cqu, M21- cosh(qL), M22- sinh(qL),
N1- -(l-C1qL)e‘qx' and N2- e-q(L-x'). The determinant of the H.matrix

is equal to,
Determinant - sinh(qL) + Cqu cosh(qL). (3.2.21)

Expanding the hyperbolic functions using equations (3.2.13) and (3.2.14)

gives,

 

Determinant - 1‘222-3—ELEEl + Cqu i-qu : §-qL), (3.2.22)
or,
Determinant - % (1 + Cqu) eqL - % (1 - cqu) e‘qL, (3.2.23)
or,

qL (1'C1qL) -2qL

Determinant - % (1+C1qL) e 1 (3.2.24)

' (1+C.qL) e

63

When 2qL is large, which is appropriate when time is small, the
second term in the brackets in equation (3.2.24) goes to zero and an

approximation to the determinant becomes,
Determinant - % (1+c,qL) eqL (3.2.25)
Cramer's Rule is used with the approximation to the determinant

(equation (3.2.25)) to solve the matrix equation (equation (3.2.20)) for

the constants D1 and D2.

_1_[‘_1fl" -q<2L+x') (”1‘11" -qx' _m—q“ -q<2L-x'>]

DI ' 2qL (1+C1qL) ° ‘ (1+c,qL) e + (1+C1qL) 3
(3.2.26)
D _ _L w) e-q(2L+X') + (_LE‘q—L) e-QX'+ 321.1; e-q(2L-x')
2 2qL (1+C‘qL) (1+C1qL) (1+C1qL)
(3.2.27)

Substituting D1 and D2 into the solution for T and expanding the
hyperbolic functions yield the Laplace transform Green's function, G,

for this case, which is,

I (1-C qL) v I
~ _ .1. -q(x-x ) ____l__ -q(X+x ) -q(2L-x-x )
Gx42 2qL [ e + (1+c,qL) e + e
(l'Cqu) - _ y _ _ I
+ m) ““2”" " ’ + e “2“ ’ J ] (31-28)

and is valid for e.2qL being small. The term qL is large when time is

small.

64

When equation (3.2.28) is expanded with partial fractions, the
coefficients of each term can be matched with a transform in a Laplace
table of transforms found in Appendix A.

The exponentials in the complete series (with no approximation
to the determinant) for boundary conditions of the first and second kind
have simple coefficients for the transforms that lead to solutions that
are valid for all times, but for large times, the solutions converge
slowly.

For the more complicated boundary condition of the third, fourth
or fifth kind, and for larger times with the first and second boundary
condition, the coefficients of the successive exponentials in the com-
plete series for the transforms become more complicated functions of
q2L2; hence only the first few terms of the series are readily used and
the solutions are valid for relatively small times. The coefficients
for the five types of boundary conditions are calculated by combining
the similar terms of the expansion. A general form is obtained below
for the small time Green's functions.

Only a few terms in the small time expressions need be generated
because the small time solution can be combined with the large time
(Fourier) solution to get a solution that is accurate and efficient for
any time. Beck and Keltner [1985] demonstrate this idea in a paper on
the time partitioning of transient heat conduction solutions. Time
partitioning of the Green's function allows the solution to transient
heat conduction problems to be more efficient since the number of terms
in the solution are tractable. The integration of a Green's function
with respect to time may have poor convergence properties when the
Green's function has not been properly partitioned.

Following the notation of Beck and Litkouhi [1985], the equation

below can be used for all cases for small times (large qL), XIJ, where I

65

represents the boundary condition on the left side (1,2,3,4, or S) and J
represents the boundary condition on the right side (1,2,3,4, or S) of

the slab ,

~ _1. .1. -qu-x'l'
GXIJ(*'x"s) a 2qL °

-q(X+x') -q(2L-x-x')

+ A(aI) e + A(bJ) e

+ B(aI’bJ) [ e'q(2L+X'X') + e'Q(2L-x+x')] ] (3.2.29)

A maximum of five terms is retained for the small time solution.
The five terms in equation (3.2.29) represents the original source and
four sources or sinks closest to the original source term as shown in

Figure 3.2. The coefficients A(-) and B(-,-) are given below.

 

 

 

_1___9__
A(c) - 2qL - qL(qL + c) (3.2.30)
_ , .1. ..l.
2qL + qL+c (3.2.31)
qL-c
- 2qL(qL + c) (3.2.32)
_1_ c + d
B(c,d) - 2qL - (qL+c)(qL+d) (3.2.33)
(qL'C) (qL-d)
' 2qL(qL+c)(qL+d) (3'2'34)
If c # d, B(c,d) can be written as,
.l. std ..l_. ..l.
B(c,d) - 2qL + c-d [ qL+c - qL+d J (3‘2'35)
If c - d, B(c,c) can be written as,
B(c,c) - -1-- -3“L—- (3.2.36)

2
2qL (qL+c)

66

+ additional reflection terms
‘(ZL-x'O-x')

_ .. {21.x .. x: ——————— x - 21.
////////////// , _ .

 

-—-x
f h

|x~x'|

_1_ T ..
//.(..\_< ,5, ///////// “ '°

+ additional reflection terms

 

 

————————_——— x-oL

+ additional reflection terms
‘(2L+x - x')

— —————— - ————— x--2L
.(ZL-n-x-o-x')

+ additional reflection terms

Figure 3.2 Reflections of Sources and Sinks in a Finite Body
and the Locations of Additional Reflection Terms.

67

The a and b values used in equation (3.3.29) are given by,

 

 

I J

a1 - a b1 - w (3.2.37)
83 - Bl b3 - 82 (3.2.39)

2 2 2 2
a, - Clq L b, - Czq L (3.2.40)

2 2 2 2
85 - B]. + Clq L b5 - 82 + Czq L (3.2.41)
where,

81 -T 132 -—k— (3.2.42)

(PCb)1 (pCb)2
c1 - (pen 0., - (ch) (3.2.43)

Equations (3.2.37) through (3.2.41) give the simplified version
of the coefficients A(-) and B(-,-) to be used in equation (3.2.29).
Table 3.3 gives a summary of the coefficients of e-q(x+x') for boundary
conditions of the first through fifth kinds and indicates the two types
of coefficients that need to be transformed.

The coefficient 5%: and the coefficient 75511—3) coupled with
the exponential term, each have a simple transform found in the table of
Laplace transforms in Appendix A. The symbol fl is a constant that
depends on the type of condition at the boundary. The Laplace Inversion
Theorem could be applied if the transform does not appear in the table

of Laplace transforms, but for small times, this is not necessary.

Table 3.4 is similar to Table 3.2 except it contains the coefficients

for e-q(2L-x-x').
Table 3.5 contains the coefficients for the two exponential
terms e-q(2L+x-x') and e-q(2L-x+x'). This table depends on the previous

two tables with the exception of nine coefficients. The nine coeffi-

cients that appear in Table 3.5 and do not appear in Tables 3.3 and 3.4

68

Table 3.3 Coefficients of e-q(x+x')

I - 0 A - 0
I - 1 A(al) ' ' 2::
I - 2 Mag) - 51..“
- - - _L _L
I 3 A(as) 2qL + QL+BI
- _ .1. - ___1___
I 4 A(a‘) 2qL qL+l/C1
5 2qL 01(51'82) QL'SI qL-82
for C1 < 2%:
where, Si - 5%— [ -1 + (1 - 43101)1/2]
1
32 " 3%: -1 - (1 - 43101)1/2}
Table 3.4 Coefficients of e-q(2L-x-x')
J - 1 A05!) _ - 5.3.1:
J .. 2 A(bz) - 531-:
- _ , .1. __1__
- - .1. _ ___1__.
J A A(b‘) 2qL qL+l/C2
- _ _ .1. __1_ 1 - 1
J 5 A(bs) 2qL + (32(83-84) [qL-Ss qL-S‘ ]
for C2 < 2%;

where, 83 - 2%- [ -1 + (1 - 482C2)1/2 ]
2

S4 ' 5% ( -1 - (1 - 4B,c,)1/2]
2

69

Table 3.5 Coefficients of e-q(2L&x-x') and e-q(2L-x+x')

 

 

J - 0 J - l J - 2 J - 3 J - 4 J - 5
I - 0 0 0 0 0 0 0
I - 1 0 -A(a,) A(a,) -A(b,) -A(b,) -A(b5)
I - 2 0 -A(a,) A(a,) A(b3) A(b,) A(bs)
I - 3 0 -A(83) A(as) B(as,b3) B(as,b‘) B(as,b5)
I - 4 0 -A(a,) A(a4) B(a,,b3) B(a,,b,) B(a,,b5)
I - 5 0 -A(85) A(as) B(as ,bs) B(as,b‘) B(as,b5)
B + B
_ .1. ..1....3 _ ___1___
where B(a,,b3) 2qL + 31 _ 32 qL + Bl qL + B, , B, i B,
23,
' 5%: ‘ 2 . B1 ' B2
(qL + Bi)
1 + 6,3,
_ _ _ ......._. , ___1.__
B(a4'b3) 2qL 1 - 0,3, [ qL + 1/0, qL + B, ] ' C131 ‘ 1
2/Ci
(qL + 1/C1)

 

B(a5,b3) - 5%: [ 1 + (31 + 52)[ 3;52 + (2-3,)(2-S,)] ]

6,8: + (B, + 3,)
(qL+S,) s,(s2 - 5,)(2 - s,)J

 

+

J 6,8: + (B, + 3,)
(qL+S,) S,(S, - 8,)(2 - 8,)

and where S, and S, are defined in Table 3.3 and C, < 1/4B,

 

 

C + C
B(a b ) - + -l--3' - --l--- c a c
" ‘ 2qL C, - C, qL + l/C, qL + l/C, ’ 1 2
2/C,
_.1.._ ’ c _ C
‘qu 1 2
(qL + 1/C, )2
B(a b )- 1 + B 1

c ,5, [31+
(qL+S) 3(82 -s )(2 - s >

 

70

Table 3.5 (cont.)

 

2
C S + B
+ 1 l 2 1
(QL+S:) 82(5‘ ‘ $2)(2 ‘ 82)

 

2
B(a b ) - —l— - (c, + c,)s, + (3‘ + B“) [ ..1._
5, 5 2qL (S,-S,)(S,-S,)(S,-S,) QL+SI
2
(C, + C,)S, + (B, + B,)
(SI'SZ)(SS'SZ)(84'S2) [

2
(C, + C,)S3 + (B, + 8,) 1
- (S,-S,)(S,-S,)(S,-Sa) [ qL+83 ]

2
(C, + c,)s, + (B, + B,) 1
- (31‘54)(S2'Se)(33‘34) [ qL+S, ]

where S3 and S, are defined in Table 3.4 and,

 

..J...]
qL+S,

 

 

B, > 0, B, > 0, c, > 0, c, > 0, c, < l/4B,, 02 < 1/432

 

 

 

2
2(CIS, + 31 ' 2 C181)
B(a5,b5) - a]; - 2 [——L—2]
($2 ' SI) (qL + SI)
2
C S
+ 1 1 [ 1 ]
(S2 ' 3,) qL+S,
2
2(6:82 + B, ‘ 2 CI82)
- 2 {—1—2]
(82 ' SI) (QL + S2)

 

2
C S
+ 1 2 [ 1 ]
($1 ' 82) qL+82
for BI - 32 , C1 - C2 , and BIC! ‘ 1/4

B(a3,b,) - B(a,,b3)

 

B, 4 B, , c, 4 c,

B(a3,bs) - B(as,b3)
B, ~ 3, , c, » c2 , s, 4 s3 , s2 4 s,

B(a‘,b5) - B(as,b‘)
B2 * B1 . C1 2 C2 . Si 2 33 . 82 T 54

71

are expanded and attached to the end of Table 3.5. Each of the nine

_L
2qL’

,, where, as before, p depends on the condition at

coefficients can be expressed as a combination of the functions

(L+fi>'or
q mL+m

the boundary. Two of these coefficients have been previously discussed
and the third coefficient can also be found in the table of Laplace
transforms in Appendix A.

Equation (3.2.29), Tables 3.3, 3.4 and 3.5, and a brief table of
Laplace transforms are all that is necessary to determine an approxima-
tion to the small time Green's function for boundary conditions of the
first through fifth kinds. Table 3.6 gives the inverse Laplace trans-
forms for the terms that are included in Table 3.3.

The four locations nearest the original source, along with the
source's location, correspond to the locations of the exponential terms
that will lead to an approximate expression of the small time Green's
function. This procedure lends to itself a simple physical interpreta-
tion. Each term in the approximate series corresponds to the solution of
a related problem for an infinite slab, see Figure 3.2, and thus the
solution for the finite region can be interpreted as the effects of
adding sources and sinks to an infinite body. Since the coefficients of
equation (3.2.29) can always be written in terms of 5%: + additional
terms, the approximate solution to problems with boundary conditions of
the first through fifth kinds will include the solution for a slab that

is insulated at both boundaries. The choice of placing a source or sink

_1.
2qL

negative. A positive value for this term gives a source at the loca-

at a particular location depends on whether the term is positive or

tion. Any additional terms associated with the location represent the

effects of a boundary condition that is not insulated. Figure 3.2 shows

72

Table 3.6 Inverse Laplace Transforms of A(-) e-qx

 

I Wafers
0 0
1 - L2 mm)
L
2 “5 EX(X.t.)
L
3 a? [ EX(X.t) - Bi ER(x,t,Bi) ]
4 - “a [ EX(x.t) - CilER(x,t,Cil) ]
L
5 - :31] EX(x.t) - ci(s,l- 5,) [ ER(x,t,S,) - ER(x,t,S,) ] 1
2
- u+*
Exam) - [a x 31'1” e “t .x* ”1‘, . t* - %
L
* 2 * *
ER(x.t.u) - e'(x ) /(“t ) rerf [ * 1 2 + u (9')”2 ]
(4:: > /

2
* ‘k
where the function rerf(z) - ez erfc(z), t - dimensionless time, and x

is the normalized coordinate.

73

the locations of the sources and sinks for a slab along with the loca-
tions of any additional terms that might occur for a boundary that is
not insulated.

The procedure may be extended to include more reflections of the
source by retaining the second term in equation (3.2.19) for the deter-
minant of matrix equation (3.2.20). A third term, C(c,d), appears in
the coefficient list of equation (3.2.29) and multiplies the reflection

- 0 _ __o
q(2L+x+x ) and e q(4L x x ).

locations e The coefficient C(c,d) is,

(qL-c)2 (qL-d) (qL-<1)2 (qL-c)
or

 

 

C(c,d) - (3.2.44)

2
2qL<qL+c> (qua) 2qL<qL+d)2<qL+c>

The third term, C(c,d), complicates the evaluation of the solution by
including additional terms and functions that are not easily trans-
formed. Recently, a paper dealing with the generalization and
application of Laplace transformation formulas in diffusion problems
[Shibata and Kugo, 1983] has eased the calculation for some of the
inverse transforms of the coefficient C(c,d), but for small dimension-

less times, the term C(c,d) is not necessary.

3.3 Green's Functions for Some Semi-infinite Cases in One Dimension

The objective of this section is to show the effects of boundary
conditions on some semi-infinite geometries. The Green's functions for
semi-infinite geometries are developed from the source solution and, for
boundary conditions of the first, second and third kinds, can be found

in Carslaw and Jaeger [1959].

74

The basic building block Green's function for infinite and semi-

infinite bodies (called the source solution by Carslaw and Jaeger [1959,

p3- 50]) is.
G I '1/2 v 2
(x,x It,r) - [a « a(t-r)] exp[-(x-x ) /(a a(t-f))]. (3.3.1)

This function represents a unit impulse occurring at time r and
at position x'. In an infinite or semi-infinite medium the impulse has
an effect on the medium.a long time after the impulse is generated.
This effect is shown in the figures below.

Figure 3.3 is a plot of the Green's function of a semi-infinite
body normalized with respect to the position of the point of interest
for boundary conditions of the zeroth, first, and second kind. The
impulse occurs at x'-O and r-O. One curve is the function l/J(4t*),
which is the leading coefficient of the Green's function for an infinite
body, where the symbol t* is the dimensionless time. The Green’s func-
tion for an infinite medium, represented by n - O, is, as expected,
about one half the value the leading coefficient curve except at small
dimensionless times because the impulse can move in two directions.
When n - 1, which means that there is a boundary condition of the first
kind occurring at the surface of a semi-infinite body, the Green's
function is zero for all times. This means the effect of an impulse at
the surface of a body with a boundary condition of the first kind at
that surface is zero.

When n - 2, which means a boundary condition of the second kind
occurs at the surface of a semi-infinite body, the curve is twice the
value of the infinite curve after a dimensionless time of five. The

boundary condition of the second kind reflects the impulse which is the

75

ON

.vcax nsooom use
.umuam .nuouON one mo mcouuuvcoo husvcaon nuu3 zoom
ouacuucu-uaom a you made manuo> soauocsh a.:oouo

m.m enamfim

 

 

wafiaucv

 

 

red

 

Ind

76

cause for the doubling. All of the curves, except when n - l, converge
very slowly as the time increases.

Figure 3.4 is a plot of the Green's function for a one dimen-
sional, semi-infinite body with a convective boundary condition (X30) at
the surface, normalized with respect to x'. The three curves show the
effect of increasing the Biot number on the Green's function. The
l/Sqrt(4 t*) curve is shown for reference. If the Biot number is very
small, the Green's function approaches the Green's function for an
insulated case (X20) as expected. As the Biot number increases, the
Green's function decreases until it becomes zero, which represents a
boundary condition of the first kind.

Figure 3.5 is a plot of the Green's function for a semi-infinite
body in one dimension with a nonconvective thin film at the surface
(X40) and normalized with respect to x'. The three curves represent the

effect of increasing the Carslaw number, C which is the thermal

1.
storage capacity of the thin film divided by the thermal storage
capacity of the solid. When the thermal capacity of the thin film
approaches zero, the Green's function approaches the Green's function
for an insulated body (X20).

Figure 3.6 are the Green's functions, normalized with respect to
the source location, x', for various positions of interest and dimen-
sionless times. Curves for an infinite medium (Figure 3.6a, X00) and a
semi-infinite medium with boundary conditions of the first (Figure 3.6b,
X10) and second kind (Figure 3.6a, X20) are shown. An important feature
of these plots is the shape of the Green's function at very small times.
Notice for t* < 0.05, the shape of the curves are independent of the

type of boundary condition that occurs at the surface when the location

of the source and the point of interest coincide.

77

 

.OHQUGEQHQQ 050 HO mflaflfi> Maowkfl> UCQ

Atom
6 u o a mo souuwnsoo muoncsom s £ua3 .
Muww«wcwm%aom a you made asnuo> coauocsh m.:oouu e n ouamwh

 

 

cH.o u .n

 

 

 

 

 

OSXb * x

0‘01 ‘OO'I ‘0t°o - ‘g

78

.auouoamuom ecu mo noSHs> nsouus> one
usax Samson ocu.mo Goduuncoo amonaaom s Jags >vom
ouHCHu:«-«aom c you mafia manuo> coauoSSh m.coouo n.n ousmum

 

 

Amiga . ..u
om m— or m o
b L L — h b b n — n n n n . — n n h P b 000
92 .. SQ
I ...o
u 84 .. .o
1 ad
26 .. ~o
I Ind
I «\Hasucv u
lall

 

 

 

 

.px

O’IXs

0.0I ‘00.! IOI.° _ to

79

 

 

.hcon ouaauwmwomw sm.m enamHm
moasouucoaun mauuo> aoauossh u.
a
new :ouuumom

.x\x

 

 

 

 

N)

80

 

.25. um“:
on”. no couuwncou Dungeon a 5:3 anon ouacfih as
you flouuamom mnoacoamcoaun nauue> :ouuocam “...—Bone n56 ounwwh
4 \fi\”

n N
_ .

 

 

 

l‘)

81

.ncax ncooom
osu mo cowufincoo huencaon a sun: anon euucah so
How coauamom mmoaaowmcoaua unuuo> coauocsm n.:uouo ow.n ousmwh

LS"

n
L

 

 

 

 

 

 

 

F)

82

Another important feature of these curves is the difference
between the curves when the time is not small and the point of interest
is at the surface (x - 0). For boundary conditions of the first kind
(X10), the Green's function is zero for all times. The Green's function
for the zeroth boundary (X00) is one half the Green's function for a
boundary condition of the second kind (X20). The result is expected and
it has been shown previously.

Figures 3.7a, 3.7b, and 3.7c show the Green's functions for a
semi-infinite body, normalized with respect to x', with a boundary
condition of the third kind on the surface (X30). The three plots
represent the Biot number increasing by factors of ten from 0.1 to 10.
When the Biot number is small, the X30 Green's function approach the
Green's function for a boundary of the second kind (X20), as expected.
When the Biot number is large, the X30 Green's function approaches the
Green's function for a boundary condition of the first kind (X10).

Similar conclusions can be observed from Figures 3.8 and 3.9
which are the Green's functions for a semi-infinite body with boundary
conditions of the fourth (X40) and fifth (X50) kind occurring at the

surface with various parameter values.

3.4 Small Time Green's Functions for Finite Bodies

The objective of this section is to use the general results of
Section 3.2 to generate approximate Green's functions that are accurate
and efficient at small times.

Two example problems will be discussed in this section. Both
problems involve a one dimensional slab with constant thermal conduc-

tivity, k, constant thermal diffusivity, a, and no heat generation. The

83

.H.o I an .ncwx nausH
one no :oHUHUCoo husncaon a and: Atom ouacum as
How Gouuunom «moasoumcoauo msmuo> :ouuo::h u.:oouo an.n ounmum

 

 

 

08x9 ¥ Ix

 

 

 

I'O - ‘a

84

.°.H I an .ncuu chunk
ecu mo couuancoo humncsom a sud: Anon ouucum as
you :oauwuom mmeafioumcoaan msmuu> coauocsm m.:oouo An.n ouswum

.x\x

._V n N F
_ _

 

 

 

 

"3

00'1 - 's

85

.o.oH I Ha .ecux guess
2.3 we acumen—woo ESE—don s 558 boon ouucum as
you couuumom mmoacoamcuaua n5uuo> couuocam rescue and spawn...—

.x\x

 

 

 

 

 

"3

00'01 - ‘a

86

.H.o I Ho .ccux :uusoh
0:» mo sowuwncoo hhsncson s suds anon ouasum as
How coauunom mmoaconSuEuo mamue> :ouuocsm m.:oeuo em.n annuam

.x\x

 

 

 

 

 

 

F)

OI'O - ‘9

87

o4 I 8 6:3. 556m
ecu mo acupuncoo huupcson s nuuz hvon uuucuh no
you coauumom mmoHcoumcoaqo nauuo> couuocsm m.:ouuo Am.n shaman

.x}.

n

 

 

 

 

 

F)

OO'I - ‘3

88

.o.o.n I H0 Jon—«x Sunfish
osu mo :oauancoo husnason u suds anon ouucum so
you :oauunom amoacoumnoaaa namuo> coauocsm m.:eouo om.n shaman

.x\x

 

 

 

 

 

 

F)

O’IXo ¥ ,x

0'01: - ‘o

89

.H.o I Hm .~.o I so .ncux smug“
ecu we sowuwvcoo husncson m Sufi: anon ouwflnh so
you couuumom mmoacoumsoaun msmuo> couuocsh micoouo sm.n ouswuh

._V m

—1—
O

 

 

 

 

N)

 

01.0 - to ong 4‘ ,X

1'0 - ‘9

9O

.c.H I an .H.o I so .ncfix suwum
0:» mo :ouuuncoo husmGSon a Sums anon uuwcwm as
you couuamom nmoaaoumcoaun annuo> couuocsm m.coouo Am.m ouswuh

.X\x

 

 

 

 

 

1')

01'0 - ‘9

0'1 - ‘s

91

 

.o.OH I ah .H.o I do .n:«x Sauna
on» we coHanCOU hucncaon a saw: anon ouafiam as .
you cowuamom mmoHGonCoawn n5nuo> :ouuocsh n.coouo om.n shaman

.xxx

 

 

 

 

oro - ‘0 “XS at .x

0'01 - ‘3

92

transient heat conduction equation that describes this case is given in
equation (3.2.4). If r i 0 in the following expressions, substitute for
the time t the value (t - r). The first example has insulated boundary
conditions on both sides (Figure 3.10a) and the second example has a
nonconvective thin film (Figure 3.10b) at x - 0 and is insulated at x -
L. The first example is called a X22 case since it has boundary condi-
tions of the second kind on either side, while the second example is
called a X42 case since it has a boundary condition of the fourth kind

on the left boundary and is insulated on the right boundary.

3.4.1 WISH).

For a slab insulated on both boundaries (X22), the coefficients

from Tables 3.3, 3.4 and 3.5 are,

-_L __l_
A[0] 2qL and B[0,0] 2qL
since a1 - 0 and b2 - 0. Adding the five terms of the Laplace transform
solution for small times for the Green's function, equation (3.2.29)

gives an approximate solution of the form,

~ _ L _1_ -qu-x'| -q(x+x') -q(2L-x-x')
6X22 a [ 2qL ( e + e + e

+ e-q(2L+x-x')+ e-q(2L-x+x') ) ] (3.4.1)

Using a table of Laplace transforms (Appendix A) or Table 3.6
for the inversion of these types of Laplace transform gives an ap-

proximation to the small time Green's function for the X22 case as,

\O
L»)

—\NX\\\\\\\\\

 

 

\\\\\\\\\\\

N

l
O
N
I
t"

Figure 3.10s Finite Body Insulated on Two Sides.

1-

3;

\j

 

\X\\\\\\\N\.i
\\\\\\\\\\\\

 

 

...—\J
-0 I
‘ x-L

Figure 3.10b Finite Body Insulated on One Side and With a Non-
convective Thin Film on the Opposite Side.

94

2 2
6x22(x.x'|t,0) ..t [a t (a t/L2)]-1/2 [e'(xL'xL) /(4 a t/L )
2 2
+ e'(Xfoi)2/(4 a t/Lz) + e'(2'xL'xi) /(4 a t/L )

2 2 2 2
+ e'(2+"L“L) /(4 " t/1‘)+ e’(2"‘L+"L) “4 a VI“ )] (3.4.2)
_ x
where xL L'
2
If the dimensionless time (a t/L ) is less than 0.1, this func-
tion is accurate to five decimal places. The Fourier or long time

Green's function to the X22 case from Beck [1986] is,

6x22(x,x' It, 0) - %[ l + 2 gle m2 x Hat/L cos(m « x/L) cos(m x x'/L) ]

(3.4.3)

A plot of the number of terms necessary for convergence to within
0.00001 versus dimensionless time (a t/Lz) is shown in Figure 3.11 for
the X22 case when the point of interest and the location of the source
term are both at the left boundary. Notice that as the dimensionless
time gets small, the number of terms in the small time solution goes
down, while the number of terms in the long time solution goes up. By
keeping the dimensionless time less than 0.025, the number of terms in
the LT method is reduced to three. When the Green's functions are
transformed into real time and used in the formalism to determine tem-
perature distributions, the integration over 1 is possible. Additional

terms leads to integration over 1 that are in closed form.

.Hococ.o u aoeuaooe .Aamxv anon
noueasmcH as you easy mmoaco«msoaaa msmuo> acuuocsm
m.:oouo one no oocomuo>soo you mamas mo nonesz HH.m ousmum

To nd We ...0 0.0

 

95

cofiuocsm mafia owes
Jlll! jlmu

K

 

 

/ ..

couuumsm mafia Hamsm «Emma

mo
scandz

o I x umououcu mo wagon mo cowusooq
o I .x oouaon mo couusood #

 

 

 

96

3.4.2 51ab_Hish_a_Qaralax_Left_fieundarx_gsnditien_end

For the X42 case, the coefficients from Tables 3.3, 3.4, and 3.5

are,
A[a ] _ .1. - --l-—-
4 2qL (qL+1/Cg)
_1_ ____l___
B[a‘,bg] ' 2qL ' (qul/Cl)

Adding the five terms to get an approximate Laplace transform solution

gives,

- _ - _ ] -q(x+x') -q(2L+x-x') -q(2L-x+x')
Gx42 Gx22 (qL+1/C1) [ ° + e + 9

(3.4.4)

where 0x22 is the Laplace transform of the Green's function of the X22
case. The additional three terms represent the reflections of the
images of the sources located at (x+x'), (2L + x - x'), and (2L -

x + x'). Using a table of Laplace transforms (Appendix A) for the
inversion of these types of Laplace functions gives an approximation to
the small time Green's function for the X42 case as,

I _ _2 I _1... o '1
GX42(x’x It,0) G L EX[x+x ,t] + LCl ER[x+x ,t,C1 ]

x22
- 2 EX[2L+x-x',t] + -l— ER[2L+x-x',t,C;1] - ; EX[2L-x+x',t]
L Lc1 L
+ -1- ER[2L-x+x',t,c;1] (3.4.5)

LCI

97
+ +
where, EX[x+,t+] - [a x c+]'1/2 ex /4t .

+ +2 _
ERIx+,t+,P] - ex P + t P erfc[ x+/./4 t+ + P J t+ ]

2
and x+- x/L, t+- a t/L .

This solution is accurate to five decimal places for dimension-
2
less time (at/L ) less than 0.1.

The Fourier solution to the X42 case from Beck [1986] is,

°° -fim 22m=/L
2 - (1/L) [t +§1[; Xm(x) Xm(x') ] ] (3.4.6)
m-l
where, Xm(x) - cos(fimx/L) - (flmCI) sin(flmx/L) (3.4.7)
No - 1 + 01 (3.4.8)
Nm - ((pmc,)2 + C1 + 1)/2 (3.4.9)
fimcot(fim) + l/C1 - 0 (3.4.10)

Equation (3.4.5) converges rapidly when (at/L2) > 0.025. Small
and large time Green's functions for the Carslaw number equal to 0.10
are shown in Table 3.7. The small time solution is composed of the
small time solution for the X22 case plus six additional terms. It is
important to notice that as the dimensionless time gets small, the
number of terms in the small time solution goes down, while the number

of terms in the long time solution goes up.

98

Table 3.7 X42 Case Cl - 0.10
Small and Large Time for Various Number of Terms

 

 

 

 

 

 

 

 

 

 

 

 

Approximate Number Number
Dimensionless time small time of Large time of
Green's terms Green's terms
function function
8.92062 1 5.64245 4
0.001 7.23578 4 6.87825 8
7.23578 8 7.16274 12
3.98942 1 4.86146 4
0.005 5.23157 4 5.22249 8
5.231524, 8 5.23150 12
2.82095 1 4.18044 4
0.010 4.27584 4 4.27565 8
4.27584 8 4.27583 12
1.26157 1 2.32323 4
0.050 2.32326 4 2.32326 8
2.32326 8 2.32326 12
0.89206 1 1.70582 4
0.100 1.70582 4 1.70582 8
1.70582 8 1.70582 12
0.39894 1 0.93720 4
0.500 0.84412 4 0.93720 8
0.91489 8 0.93720 12

 

 

99

3.5 Summary

A method for generating and organizing the small time Green's
functions is given for boundary conditions of the first through third
kind, and extended to two additional boundary conditions called the
fourth (Carslaw) and fifth (Jaeger) kinds. The small time Green's
functions developed here can be formed into fundamental building blocks,
called influence functions, since they are the response of the system to
a unit point source. The small time Green's function can be used in
conjunction with the large time Green's function to generate a Green's
function that is efficient at both small and large times. The
availability of symbolic software makes the procedure for finding the
small time Green's function both attractive and efficient.

The following equation and Table 3.8 are a summary of the small
time Green's function for boundary conditions of the zeroth through

fifth kind for dimensionless time less than 0.025.

GXIJ(x,tIx',r) - [41ra(t-1')].1/2

[_ .12 _ :+ ,12 _(21_ _ ,12

[ e 4a(t-r) + M e 4a(t-r) + N e 4a(t-r) ]

+ i [ - M 01 ER(x-x',t-r,01) - N 92 ER(2L-x-x',t-r,92)

 

 

E:
+ 1/2 { $2 ER<X+X',t'T,82) ' S1 ER(X+X'9t-fISI) }
(I'ABIC1)
E2
+ 1/2 { S‘ ER<2L'X‘X',t‘T,S‘) ' S3 ER(2L‘X'X' It'TIsa) }]I

(3.5.1)

where ER(-,-,o) and Si are defined in Table 3.6.

100

2

5 0.025

Table 3.8 Small time Green's functions for a(t-1)/L

 

 

 

 

 

 

 

2

we nu nu nu nu nu nu 1. nu nu nu nu nu 1. nu nu nu nu nu 1. nu nu nu nu nu 1. n. no nu nu nu

1

an nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu 1. 1. 1. 1. 1.

1. 1. 1. 1. 1.
2 2.2 2.2 2.2 2.2 2.2
nu nu nu nu nu a. nu nu nu no no u» nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu nu
1..l 1. 1. 1. 1.

1 111111.1.1.1.1.1.1

a. nu nu nu n. n. nu nu nu nu nu nu nu nu nu n. a. nu nu nu nu nu nu nu nu nu nu nu

N" n. nu 1. 1. 1. 1. 1. nu 1. 1. 1. 1. 1. nu 1. 1. 1. 1. 1. nu 1. 1. 1. 1. 1. nu 1. 1. 1. 1.

o .- . .- . .. . . . .

Mu nu 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
...... ...-....—-.

o

N. * m. * * * *

.mlnuuusoln34501znun mnnnunolnsm

“xxxxxxwwunxmmjmmmxxxxxxxxxnmxnx

C

 

 

* Valid for any time

CHAPTER 4

TRANSIENT ONE DIHENSIONAL CANSS PROGRAM

4.1 Introduction

A program that solves linear, transient heat diffusion problems
in the cartesian coordinate system for a variety of boundary conditions,
initial conditions, and heat generation terms in one dimension is
presented in this chapter. The program includes infinite, semi-infinite
and finite slab geometries. The boundaries of the slabs are restricted
to be parallel to the cartesian coordinate axes. The program will
calculate temperature distributions of nonlinear problems for certain
combinations of physical properties as shown in Chapter 2.

The one dimensional computer algebraic, numeric, and symbolic
solution program, (CANSS), is compiled in a computer algebraic software
system called the symbolic manipulation program (SMP). SMP [1983] is
one of a group of computer software programs designed to analyze simple
or complex mathematical problems interactively or in batch mode. A
partial list of other programs that manipulate expressions and symbols
is MACSYMA [The MATHLAB Group, 1985], MAPLE [Char, et. al., 1985],
REDUCE [Hearn, 1985] and mu-MATH [mu-MATH, 1985].

The one dimensional CANSS program is comprised of a main proce-
dure, which directs the overall analysis of the indicated problem, and a

group of specialized sub-procedures or libraries. The sub-procedures

lOl

102

contain functions for developing the small and large time Green's func-
tions and include provisions for treating the boundary and initial
conditions and the volume energy source terms. Moreover, integration
of special functions is accomplished using the sub-procedures. The
temperature distributions developed by the CANSS algorithm cover the
entire range of dimensionless time.

The CANSS program runs on a DEC MicroVax under the VMS operating
system. This device is a virtual memory machine with approximately
eight megabytes of main memory. One element of a one dimensional sym-
bolic calculation of a temperature distribution uses about 30 to 40 CPU
seconds.

Section 4.2 of this chapter explains in detail the input and
operation to the CANSS algorithm. Section 4.3 shows results of three
one dimensional problems solved using the CANSS program. Section 4.4
describes some integrals that occur during the solution of one dimen-
sional problems. Section 4.5 describes time partitioning in the one
dimensional case. Section 4.6 is a flowchart/example of the logic of

the CANSS program. Section 4.7 summarizes the chapter.

4.2 One Dimensional CANSS Program

The one dimensional CANSS program is designed to generate sym-
bolic temperature distributions that are rapidly convergent for small
and large dimensionless times in a infinite, semi-infinite, and finite
1813b. The one dimensional slab may have non-homogeneous boundary condi-
tIions of the zeroth through fourth kind occurring on the surfaces that

are functions of time but not functions of position. The slab may have

103

an initial condition that is a function of the spatial coordinate and a
volume energy heat source that is a constant in time and space.
The forcing function applied to the boundaries are of the form,

n/2
fi(t) - Ci t , (4.2.1)

where n - -1, 0, l, 2, ..., i-l is for the left hand boundary and i-2 is
for the right hand boundary. Equation (4.2.1) allows the forcing func-
tions on the boundary of a one dimensional slab to be functions of time,
constant or zero depending on the value of the constants C1 and n.

The initial condition is a function of the spatial coordinate
only and is restricted to polynomials in the CANSS program. It is noted

that the CANSS program can solve some problems that have initial condi-

tion that are transcendental and of the form,

F(x) - To sin(w x/L) or To cos(x x/L), (4.2.2)

where To is a constant temperature and the coordinate x is normalized
with respect to the length of the slab, L. In the case of a semi-
infinite or infinite slab, the coordinate x is normalized by a unit
measure of length.

The volume energy heat source is constant in both position and
time. The heat source must be constant because many of the solutions to
the integrals necessary to solve non-constant heat sources that are a
function of time or position are not available.

Tzeng and Beck [1985] describe a numbering system data base to
solutions of heat diffusion. The numbering system is used to catalog

solutions and in a numerical program to evaluate new solutions. A short

 

X

5L.

Fe.

 

104

summary of the numbering system parameters used in this thesis is ex-
cerpted and shown in Table 4.1.

The first set of parameters in the numbering system give the
boundary condition types at the surfaces of the body. The second set of
parameters gives the time variation of the boundary condition. The
initial condition is described along with additional parameters for heat
generation, fin type terms, etc. The examples discussed in the next
section show how the numbering system is applied to heat conduction
problems.

Applying this system to the CANSS program yields distinct cases
for four semi-infinite slabs (X10 ... X40) and ten finite slabs (X11 -
X44). The four distinct semi-infinite cases are allowed five types of
time variation on the boundary (BO ... B4), four types of initial condi-
tions (T0,Tl,T2,T6), and two types of volume energy heat sources
(G0,Gl). The ten distinct finite cases are allowed five types of bound-
ary condition on each side (B0 ... B4), four types of initial conditions
(T0,Tl,T2,T6), and two types of volume energy heat sources (GO,Gl). The
CANSS program will give closed form expressions for over two hundred
(204) distinct problems in heat diffusion. These cases do not include
the superposition of solutions due to the linearity of the diffusion
problem.

The fact that many types of convolution integrals are not avail-
able analytically is the limiting factor for all non-homogeneous terms
in the analytical Green's function approach for obtaining temperature
distributions. The ability of the CANSS integration algorithm of recog-
‘nizing and symbolically solving the integrals is crucial to the solution
of the problem. Many integration procedures are included in the inter-
nal SMP integrator, but due to the complex convolution nature of the

Green's function approach, the internal SMP integrator may not calculate

X0
X1

X2

X3

X4

X5

BO
Bl

B2

B3

B4

T0
T1

T2

T3

T4

T5
T6

CO
01

105

Table 4.1 Key for Heat Conduction Data Base

(Excerpted from Tzeng and Beck [1985])

Boundary Conditions (Rectangular Coordinates)

Infinite Boundary Condition
T - f(t) (Dirichlet Boundary Condition)

fig - f(t) (Neumann Boundary Condition)
k fig hT - f(t) (Robin Boundary Condition)
k ii + (pcb) gE - f(t) (Carslaw Boundary Condition)
k g: + hT + (pcb) %% - f(t) (Jaeger Boundary Condition)
Time Variation of the Boundary Condition
f(t) - 0
f(t) - To (Constant)
f(t) - Tot (Linear)
f(t) - Totn (Integer Polynomial)
f(t) - Totn/2 (Rational Polynomial)
Initial Condition
f(x) - 0
f(x) - To (Constant)
f(x) - Tox (Linear)
f(x) - Toxn (Integer Polynomial)
f(x) - Toxn/2 (Rational Polynomial)
Step Change in f(x)

f(x) - Tosin(C x) or Tocos(C x) (Transcendental)

Heat Generation Source Term

g(X.t) - 0
8(XIC) - Go

106

the indicated integrals. In such cases, the integrals can be numeri-
cally evaluated.

Special integration functions have been added to the SMP program
by means of a library of integration which resides in an external file
called gngdimiint and is called by the main CANSS program to assist the
internal SMP integrator. The internal SMP integrator is invoked when an
integral is not evaluated by the external CANSS integrator. The in-
tegral is returned to the user if it is not evaluated by either the
CANSS external integrator or the SMP internal integrator. The procedure
is halted until a solution to the integral is found. If the user can
evaluate the unknown integral, the solution is placed into the external
CANSS integration procedure and the problem is restarted.

The one dimensional CANSS program begins with the loading of
external procedures that query the user for information regarding the
right and left boundary condition. These procedures set up the geometry
for the description and calculation of the Green's function. Also, the
procedures that query the user for the initial condition and the volume
energy heat source are loaded.

Two assumptions are made at this point of the algorithm. The
partition time, t1, which is defined as the dimensionless time when the
convergence characteristics of the small time Green's function require
more than three reflection terms, is set to 0.025. The value of the
partition time may be increased up to 0.10 if boundary conditions of the
third and fourth kinds are not applied to the boundaries. The second
assumption is that the forcing functions for all the boundary, initial,
and volume energy term begin at time 1 - 0.

The first step in generating a temperature distribution for a
slab is to determine the left and right boundary conditions thus

describing the correct Green's function. The one dimensional CANSS

107

algorithm allows the user to choose from five type of boundary condi-
tions; boundary conditions of the zeroth through the fourth kind. The
boundary condition of the fifth kind is not included in the one dimen-
sional CANSS algorithm because specific numerical conditions must be met
and the solution will not be symbolic.

The next step in the algorithm is the input of the initial
condition. The initial condition condition in the CANSS algorithm, by
definition, is only a function of position. For some problems, the
initial condition may be composed of functions of sines and cosines
where the argument is a function of position. The heat generation term
is input next and must be a constant.

The program has the information necessary to develop a tempera-
ture distribution for the body. The conditions of the problem are
displayed to the user and a library of small and large time Green's
function is loaded. The small and large time Green's function for the
stated problem are completely described by the boundary conditions and
displayed for the user. The Green's functions are infinite summations,
but, the summation sign is not available for display in the CANSS en-
vironment at this time.

The summation for the small time Green's function goes from
minus infinity to plus infinity. Only the most dominant terms of the
series (n - 0,:1) are retained because of the way the time is parti-
tioned. The summation of the large time Green's function begins at one
and goes to plus infinity. Again, only the most dominant terms of the
summation series (m - l,2,3,...,N) will be retained due to the way in
which the time is partitioned. The number N is always less than ten
with proper time partitioning.

An additional equation called an eigencondition must be solved

for each term in the large time temperature distribution to obtain the

108

eigenvalue. This means an additional function must be solved for each
term in the large time solution. Mikhailov, et. al., [1983] considers
the safe and efficient calculation of the eigenconditions for the eigen-
values to be tricky, and in many cases, avoided and ignored. A Newton-
Raphson iterative procedure may be used to estimate the eigenvalues from
the eigencondition but, for boundary conditions of the third, fourth,
and fifth kinds, the iteration process adds additional calculations and
decreases the efficiency of obtaining the temperature distributions. An
additional problem with the Newton-Raphson iteration process is the
choice of the initial guess. The eigenconditions are functions of
transcendental functions and the appropriate initial guess is critical
for obtaining the eigenvalue. Beck [1986] lists some procedures for
estimating or approximating a good guess for the eigenvalues in heat
transfer problems.

The number of terms in the small time Green's function is kept
small because the integration of certain convolution integrals in time
is difficult or not possible for some terms in the series. The large
time Green's function allows integration on convolution integrals for
all terms in the series but will typically need an additional calcula-
tion of the eigencondition for the eigenvalues.

Once the Green's functions for small and large time have been
determined, the CANSS routine begins to direct the integration of the
four terms in equation (2.3.14) that will lead to the temperature dis-
tribution. The integrands of the integrals are checked against the
external CANSS integrator and a temperature distribution is generated.
If the external CANSS integrator can not recognize or match the in-
tegrand, the internal SMP integrator is tried. If both integrators
fail, the procedure halts because the integration cannot be performed in

the CANSS environment and the integral is returned to the user for

109

evaluation. A list of special integrals included in the CANSS program,
but not readily available, is provided in Appendix B.

The user must decide at this point if the integration is pos-
sible, and if so, extend the CANSS external integration library to
include the new integral. The procedure is then restarted from the
beginning. When a solution is presented, it is a function of the posi-
tion and time and has convergent expressions for small and large times.
The temperature distribution is not exact when only a few of the terms
in the appropriate series are retained but it is accurate to within four
decimal places of the exact distribution. This should be sufficient for
most applications.

The one dimensional CANSS program treats the infinite and semi-
infinite geometries as special cases. The infinite and semi—infinite
problems have only one form of the solution and are not time parti-
tioned. There are no convergence problems and the solutions to the
infinite and semi-infinite geometry problems are exact. The small time
solutions for the finite cases reduce to similar expressions as for the

semi-infinite body for small times.

4.3 One Dimensional CANSS Examples

Three example problems solved by the CANSS program are presented
in this section. The first example is a semi-infinite slab with a
constant heat flux, qo, at the surface x-O. There is no internal heat
generation in the body and the initial temperature condition is zero.
The Green's function for this problem is assigned the number X20 because
of the boundary condition of the second kind on the left hand side of_

the body and a natural or zeroth boundary condition on the right hand

110

side; see Figure 4.1a. Beck and Litkouhi [1985] and Tzeng and Beck
[1985] define a numbering system that assigns the number X20B1T0 to the
solution of the example.

The second example problem is a finite slab with a constant heat
flux on the left surface, qo, and a constant temperature, To, on the
right boundary; see Figure 4.1b. The initial condition is a transcen-
dental function of the space variable x and no heat is generated in the
slab. The Green's function is assigned the number X21 and the solution
is assigned the number X21B10T6.

The third example problem is a finite slab with heat flux a

c 1/2

function of the square root of time on the left surface, qo (E ) , and
o

a temperature condition equal to zero on the right hand side; see Figure
4.1c. The initial temperature is zero and there is no volume energy
heat source in the slab. The Green's function for this case is X21.

The solution is assigned the number X21B30T0.

4.3.1 Semi-infinite Example Problem (XZOBLTO)

The mathematical description of the semi-infinite problem is,

2
fl—% - 1 Q1 , for x > 0 and t > 0, (4.3.1)
0 at
6x
where,
- 3% - qo , for t > 0, (4.3.2)
x-O

lll

 

Figure 4.1a

Semi-infinite Body
With Constant Heat Flux
on the Surface (XZOBITO).

 

 

Initial Condition

1

cos(n x/L)

T
o

 

 

W...

Figure 4.1b

Finite Body With Constant
Heat Flux on the Surface
and a Transcendental Initial

Temperature (X21B10T6)

v

q ' q

 

: ...................
uueeuuu
x/L

 

 

 

SHH

/\.

——4II-

_—

P

q ' flak/co)15 *1 Figure 4.1c
Finite Body With a
- I B t Flux Condition as a

ea

L_.-J

Function of Time (X22330T0).

 

\\\'\\ \\V

 

 

t"

112

T(x,0) - 0 , for x > 0. (4.3.3)

Carslaw and Jaeger [1959, pg. 75, #6] show the exact solution to

be,
1 2

2 go (at) /

v(x,t) - ____k___— ierfc[ 2—7%;E) ]. (4.3.4)
The CANSS program returns the temperature distribution,
T(x,t) - Lfl—L (19")1/2 Gamma[l] IErfc[ Z—XL; , 1 ]. (4.3.5)
Jt
*
In the CANSS program, t is the dimensionless time, -%, and the symbol
L

xL is the dimensionless position, i.

The symbolic solution from the CANSS routine matches the solu-
tion found in Carslaw and Jaeger exactly since Gamma[l] - F(l) - l and
IErfc[y,1] - ilerfc(y) - ierfc(y). A plot of the dimensionless tempera-
ture distribution in the semi-infinite slab is shown in Figure 4.2.

This figure shows the temperature decreases at a specific time when the
point of interest is moved into the body. The effect of the heat flux

is felt instantaneously at all points in the semi-infinite slab, though
when the point of interest is far from the source at the boundary, the

effect is insignificant.

The repeated integrals error function are important in heat
diffusion problems because they appear often. The integral of the error

On 0
function, 1 erfc, is

inerfc(y) - I in'lerfc(£) d5, (4.3-6)
Y

113

.xSHm use: uceumcoo new: Assn ouucumsu-uaom a you
came muoacoamcoaqa manuo> eunusuodaoa «meadowmaoaun ~.¢ oHSMuh

 

A
Adwqu I «u eaqa emoacounsoauo
0. — 0.0 0.0 *0 «.0

D D D '— D D D - D D D - D D D - D
1
.. o; ... x

0n.0 u x
r. m~.0 u x
To
0.0 u x

r

D L D n D D D h D D D - D D D - D D

 

rm.0 A<Aevvu

 

 

114

for n - 1,2,3..., and ioerfc(y) - erfc(y), where erfc(o) is the com~
plementary error function. Integrating equation (4.3.6) by parts and

letting n-l yields,

2

Y
ilerfc(y) - ierfc(y) - £7;— - y erfc(y). (4.3.7)

A general recurrence relationship for the integral error func-

tions [Carslaw and Jaeger, 1959] is,
.n .n-2 n-1
2 n 1 erfc(y) - 1 erfc(y) - 2 y i erfc(y). (4.3.8)

A plot of the integral error function, inerfc(y), as a function of n and
y is shown in Figure 4.3. The integral error functions approach zero
quickly as the argument y increases and the order of the index n in-
creases. Solutions to heat transfer problems which contain integral

error functions tend to converge quickly.

4.3.2 Finite Slab with arTranscendental Initial Condition (X21B10T6)

 

The second example problem is a finite slab of length, L,
without heat generation. The slab has a constant heat flux, qo, at the
boundary x - 0 and a constant temperature of zero at x - L. The initial
condition of the slab is a periodic function of space. The mathematical

description of this case is,

o;
r-i

, for 0 < x < L and t > 0, (4.3.9)

QIH
°’|
r1

6x

115

 

.: mo couuossm s we :ouuucsh uouum Heywoucm

n.q ousmHh

 

 

 

 

116

where,
- 31 - qo , for c > 0, (4.3.10)
x
x-0
T(L,t) - 0 , for t > 0, (4.3.11)
and,
T(x,0) - Tocos(« x/L) , for 0 < x < L. (4.3.12)

The temperature distribution from the CANSS routine for this
case is composed of two parts. The first portion, denoted X21B10T0, is
the effect of the non-homogeneous boundary condition on the temperature
distribution. It is composed of a small time function, which is ac-
curate when the dimensionless time is less than the partition time, t1,
and a large time function, which is accurate when the dimensionless time
of interest is greater than the partition time. The small time tempera-

ture distribution is,

qoL * 1 2 ” |o.5 (2 n + xL)|
T(x,t) - 2 [ -k— ] (t ) / F(l) E IErfc[ (t*)1/2 , l ].

 

n--w

(4.3.13)

while the large time temperature distribution is,

 

q L m [0.5 (2 n + x )]
T(x,t) - 2 [ -%f-] (651/2 r(1) E IErfc[ L , 1 ]

(cf 1/2
90
+ 2 [ _k_ ]

n--m

l"

E Cos[(m - 0.5)n xL]
m-

l

117

2 2 2 2
g-(m - 0.5) x cf _ e-(m - 0.5) x c*

 

2 2 , (4.3.14)
x (m - 0.5)

where xL is the position of interest normalized with respect to the
length of the slab. The index n in equations (4.3.13) and (4.3.14) need
not extend from minus infinity to plus infinity. Only three terms are
necessary (n - -l, 0, 1) for the distribution to converge when the
partition time is chosen correctly. The index m in equation (4.3.14)
may be kept small (m s 10) for efficient calculation of the temperature
distribution at large times.

The second term in the solution to the temperature distribution
(X21T6) is caused by the non-zero initial condition. Time partitioning
is not always needed for the initial condition. The second term, gener-

ated by the initial condition, is,

co 2 2*
2 To E e-(m ' 0'5) " t Cos[(m - 0.5)« x] .
m-l

(4.3.15)

grub—17.7.1271]

The temperature distribution of the slab is determined by adding
either equation (4.3.13) or (4.3.14) to (4.3.15) depending on whether
the dimensionless time is greater than or less than the partition time.
The temperature distribution of the second example problem has not been
previously determined. The temperature distributions of the slab for
equations (4.3.13) and (4.3.14) are plotted in Figure 4.4a for the
appropriate dimensionless time. The partition time is set to 0.10 and
splits the region into two areas for problems with boundary conditions

of the first and second kind.

118

.couu :o
humocsom can you 038 oamsmxm cu zoom oUHCwm “chow

oafih mmoH:0umcoEuo msmuo> ouzuouoaaoh amoucoamCoaaa ac.c ousmwh
A

AMMMHN I a» mafia mmoHao«mCoaua

04 md od 30 Nd 0.0

— \h In — o
p n . .‘b b n n L u 00

 

... r to

”.2
m.nwnu
Inc
.. . rmd
.. Hod
.. ....E
.. . Ind

:owuavaoo bump—30m I a. O

 

 

 

 

 

L b n — n n u n p n u — n p n b n n - OoF

Ax\qouva

119

Equation (4.3.15) is plotted in Figure 4.4b for the constant To

QOL
- _k- so that both curves can be plotted on the same scale. The tem-

perature distributions for the boundary source (Figure 4.4a) look
similar to the distributions for the semi-infinite case in Figure 4.2
but the temperature distributions in Figure 4.4a converge to a steady
state temperature due to the imposed boundary conditions.

Figure 4.4c shows the temperature distribution in the slab when

the boundary condition solution is added to the initial condition solu-

QOL
tion and To - -k—' The nature of the initial condition solution's

effect is quickly damped as the time increases.

The temperature distribution in Figure 4.4c when the location of
the point of interest is at fo 1 shows a small error. The exact solu-
tion at the point xL- 1 should be zero but the graph shows a value of
0.0036. The error occurs when the small time function is required to
calculate the approximate solution without regarding the effect of the
boundary on the right hand side. Since the point of interest is the
boundary on the right hand side the small time function picks up a small
amount of error because additional terms in the small time series become
significant. The approximation at this boundary can be readily improved
by reducing the value for the partition time or adding more terms to the

small time series.

4.3.3 Finite Slab With a goundary Condition a Function of the
S u r Root 0 Time X B30 0

.The third example, denoted X21830T0, is a finite slab of length,
L, with no heat generation, zero as the initial condition, heat flux a

function of the square root of time on the left boundary and temperature

120

.cossuccoo

HmwuwcH osu you 038 ofiaswxm cw anon ouficwm m you

oEHH mmoacofimcuswa msmuo> ousumuomsoh mmoHCOMMCoEMQ nq.q ouswam

 

 

m.o I x

o.« I x

cesusncou Hmsuscs

 

 

Ax\s°ccw
.Am4Mqu

121

.umzuowOF copc< COwUMpcoo Aumcczom cam
HmauscH was you ozs usaemxm as soon mousse m too

oEfiH mmoacomeoEwn msmuo> ousumquEoe mmoaconcoEwa

m.o I x

 

A
N I U

Auquw *

0.0 To Nd

o.H I x

acaumoCoo humpcaom + acuuwpcoo HouuwcH

. L . . . _ p . . _

oq.¢ ouawwm

0.0

10.0

 

 

 

Ax\qoav~

122

set to zero on the right hand boundary; see Figure 4.1c. The mathemati-

cal description of this case is,

2
fi-§ - 1 i1 , for o < x < L and c > 0, (4.3.16)
6x a at
where,
-k g1 - qo(t/to)l/2, for t > 0, (4.3.17)
x
x-O
T - o , for t > 0, (4.3.18)
and,
T(x,0) - 0 , for 0 < x < L. (4.3.19)

The constant to is an arbitrary unit of time.
The partitioned temperature distribution cannot be calculated by
the CANSS program in closed form. The integral expression for the

temperature at small dimensionless times that remains is,

(I)

qo L * I0.S (2 n + xL)I
T(x,t) - 4 [ * ] F(l.5) t E IErfc [ * 1/2 , 2 ]
k Jto (t )

 

 

n--m

(4.3.20)
and for large dimensionless times is,

2
-(2n+xL)

no t .
‘1 L l 1 4 A
T(x’t) ' [ o * ] { Jr E I (J ' ”U2 6 /A ‘u

 

123

2 2 *
-(« (m-.5) c,)

Q cos(«(m-.S) ) e
+ [ 2 E XL

3 3
m-l t (m-.5)

[ («(m-.5)) (c* - cf)1/2 - Daw[«(m-.5)(t*- cf)”2 ] ] },

(4.3.21)

where Daw(-) is the Dawson integral and is discussed in the next sec-
tion, t: is the dimensionless partition time, and t* is the
dimensionless time of interest.

The integral in the large time expression for the temperature
distribution has not been evaluated in closed form. Many approaches may
be used to evaluate this integral approximately. The most significant
approach for evaluating the integral in a symbolic context is to define
the integral as a new function. This new function would be evaluated
for all values of the parameters and would become a well known function.
The exponential integral function, En(z), is an example of defining an
integral that cannot be expressed in closed form and is considered a
well known function.

Another approach for the evaluation of the temperature distribu-
tion at small times is to approximate the unknown integral using
numerical integration. The solution will be dependent on the specific
parameters associated with the problem. Since the solution can be
partitioned in time, the integral in the large time expression can be
evaluated using a small number of terms in the infinite series.

A third approach for the solution of the integral in the large
time expression is to approximate the forcing function over the range of
integration for small times. Figure 4.5 plots the forcing function f(r)

' (T/to)1/2

from the third example problem for a time of interest

124

mamuwousn oEHH HHmEm now mCOuumafixouaa< me awake 038 m.¢ ouswfim

 

 

=o_uas_xotgq<
ucmpmcoo

 

 

:owmom ms_h __cEm IIIIy

 

IIIIIIIInu|.:o_mmm me_b magma

N\3ios\.c . 1.3t

:o_ume_xocaa<
Lamcwm

 

     
  

 

 

125

greater than the partition time. The small time portion of the forcing
function could be approximated by a constant value of the function
averaged over the range of integration or more accurately by a linear
function. A closed form for the solution of the small time integral can
be found for a constant or linear forcing function.

The small time portion of the solution may be approximated using
only the first three terms (n - -1, 0, l) of the series when the dimen-
sionless time is small. The n - 0 term will dominate the solution at
small times. The index m, for large times, will also remain small
because the first few terms of the large time series will dominate the
solution.

The large time expression for the temperature distribution
represents a new solution that has not appeared in the literature pre-
viously. The second term in the large time solution for temperature
distribution will represent the total solution when the partition time
is set to zero. This solution is exact but needs many terms in the

series to converge.

4.4 Sale Integrals Used in.0ne Dimensional Problems

The objective of this section is to discuss two types of in-

tegral functions that occur during the calculation of the temperature

distribution of a one dimensional slab. These integrals are the Dawson

integral and an integral which occurs calculating large time solutions.

4.4.1 The Dawson Integral

126

The Dawson integral, F(x), discussed and analysed by Dawson
[1898], was first analysed in heat diffusion problems by Gordon and

Miller [1931] and is of the form,
x

F(x) - e"x I eu du. (4.4.1)
u-

The Dawson integral can be approximated for all ranges of its
argument as a function of the confluent hypergeometric function.
Lebedev [1965, pg. 20] shows the Dawson integral satisfies the linear

differential equation,
F'(x) + 2 x F(x) - 1, (4.4.2)

with the initial condition, F(O) - 0. When a series expansion of F(x)
is substituted in equation (4.4.2) and the coefficient of similar powers

of the argument x are collected, the expansion yields,

F(x) - (-l)m 2m x2m+l
1-3-0-(2 m + l)‘
m-0

 

-m < x < co (4.4.3)

Equation (4.4.3) can be expressed as,

m 2m

C-llm 2 x
“’0 ' x E 1-3---(2 m + 1)' (“'4'“)

 

m-O
Which is a form of the confluent hypergeometric function, ¢(.,-,-),

3

F(x) - x ¢(l,2,-x2). (4.4.5)

127

The properties of the confluent hypergeometric function can be found in
Abramowitz and Stegun [1964, Chapter 13].

The Dawson integral provides the solution to one dimensional
diffusion equations for finite bodies when the non-homogeneous forcing
function on a boundary is a function of the square root of time. A plot
of the Dawson integral versus its argument is shown in Figure 4.6. A
disadvantage of this integral is its slow convergence as the argument
gets large but since a difference in time is needed, the difference
converges quickly. The temperature distributions obtained for finite
bodies with boundary conditions a function of time to the one half power
or, more generally, to the 3 power, where n is an odd index, have not
been previously discovered and represent a new type of solution. The
argument to the Dawson integral is a function of the eigenvalue of the
diffusion problem that increases in value as more terms are added to the
series. This helps the solution convergence.

When the forcing function is tl/Z, t3/2,..., t<2n-1)/2, the

Dawson integral is involved and it is divided by the eigenvalue to an

odd power other than one. The general equation to be solved is,

t-t1
In(t,t1) - I 1(2 n - 1)/2 exp( - flza(t-r)/L2 ) d7 .
r-O
for n - 0, l, 2, .... (4.4.6)

When n - 0, the solution of the integral is,

_fi. .13. '1/2 -328 * *1/2
a 5m 0 flm

(4.4.7)

r-S'I'

.Hmumoucn Cowman cé 95m“..—
x

.4. m 2 o
b s , p _ 0.0

 

128

 

Tfio

Tad

rnd com

 

 

129

at
* a E 1
where t - 2 and t1 - “—3 .

L
The general solution to equation (4.4.6) when n - l, 2, 3,

is,
2 2 .iZn;ll 2 * 2 * * £22311
In(t,t1) _ [ _L_; ] [ _L_; ] 2 e-flmtl [ flm (t - t1)]
a fim a fim
- (Zn-l) In_1(t,t1) (4.4.8)
4.4.2 A n a te 0n D me oble s

An exponential integral that occurs in solving the one dimen-
sional heat diffusion problem using a large time Green's function
approach when the forcing function is a function of time raised to an
integer power is presented in this section. The general form of this

integral is,

t'tl
2 2
Zn(t,t1) - I 7“ exp( - fima(t-r)/L ) dr , (4.4.9)
r-O
where n - O, l, 2, ..., flm is the eigenvalue associated with the

eigencondition, and t1 is the boundary partition time between small and
large times.

As an example, consider the case of n - l,

t-t1
Zl(t,t1) - I 1 exp( - flza(t-r)/L2 ) dr , (4.4.10)
r-O

130

may be written as,

2 * t‘tl 52 21
-6 t m 2
21(t,t1) - e m r e L dr , (4.4.11)
r-O
*
where t - 2—% . Integrating by parts, letting,
L
1928-;
u - r and dv - e m L dr
2 2 21
1 5m 2
du - dr and v - [ 2 ] e L d1
a flm

yields,

l"

2a=_r.
t-tl — I eflIn 2 dr ].
r-O r-O

 

(4.4.12)

Performing the indicated operations gives,

 

(4.4.13)

* * O 0
where t and t1 have been described in the prev1ous section. The number
of terms necessary for this expression to converge to six decimal place

accuracy is three.

131

A general expression for the solution of integrals of the type

in equation (4.4.9) is the recursion relation,

2 2 n -fi2t* 2 * * n
Zn(t,t1) -[-L—2] [L2] e m1 [5m (t-t1)]
a 5m a 6m
- n zn_1(t,c,) ] (4.4.14)
where,
2 -fiztf -fl t*
20(c,c1) - [-—L—2 [e m -e m ] (4.4.15)
a flm

While the general expression for the integral in equation
(4.4.9) does not represent a new type of solution, it does allow the

expression to be presented in a compact form and is easy to program.

4.5 Time Region Partitioning for One Dimensional Problems

Solutions to linear, transient heat conduction problems in the
cartesian coordinate system must be split into solutions that are con-
vergent for the dimensionless times that are specified. The lack of a
convergent solution causes an unnecessary use of computation resources.
A "tuned" solution, one which is optimized for speedy convergence
characteristics, is much preferred over a solution that converges
slowly, i.e., required hundreds or even thousands of terms.

The solutions to heat diffusion problems for temperature dis-
tributions in cartesian coordinate systems can be tuned readily using
the Green's function approach. The solutions using the Green's func-

tion approach will be tuned for small and large times because the

132

Green's function for heat diffusion problems have already been placed in
the convergent form for small and large times.

The Green's function solutions are generated for one dimensional
cartesian coordinate systems, but it has been shown in Chapter 2 that
the one dimensional Green's functions may be multiplied together to
obtain multi-dimensional Green's functions for certain boundary condi-
tions. This creates some minor adjustments in the procedure to
partition the dimensionless time for the spatially multi-dimensional
solutions which wiil be discussed in the next chapter.

In a one dimensional problem of heat diffusion, the dimension-
less time where both the small and large time Green's function yield
acceptable convergence rates will be called the partition time, tf,
where the star (*) denotes dimensionless time. There exists a Green's
function solution that is tuned for convergence for each time parti-
tioned region. In the one dimensional case, the two time regions are
separated by the partitioning time, cf, such that in the small time

region,

*
t -t1_<.r 5:, (4.5.1)

*
O s r < t - t1, (4.5.2)

*
where t is the dimensionless time of interest, see Figure 4.7. The
integration of the Green's function for times that fall in the region of

small time is,

.oauh mo
:ofiuoczm m can :owuocsh m.:oouo ozu mo cauusao>coo 5.: Guzman

 

133

$334.50

4 I CC

Asuu..x.xvc

A

 

seamed oaub omumd "1 Asvm

 

.IIIIII.:o«mom mafia aaoam

 

 

 

 

 

134

*
t
* * * *
Is - I GFS dr for t - t1 5 r S t , (4.5.3)
* * *
T -t -t1

and the integration of the Green's function for times that fall in the

large time region is,

* * 'k
*
I2 - I GFS dr + I GFL dr for 0 S r < t - t1, (4.5.4)
* * * *
f -t ’t1 7 -0

where CPS and GFL are the small and large time Green's functions respec-
tively that have been integrated over the surfaces or volume of the body
and include forcing functions that are functions of time. Of particular
interest is that the small time Green's function is quickly convergent '
for dimensionless times at or close to zero but slowly convergent for
large times.

Due to the choice of ti, the integral Is is quickly convergent.
The large time Green's function converges very slowly when the time, t*,
approaches zero, but since the time integration of the second term on
the right hand side of I! does not involve zero or a time close to zero
by the definition of CT, the integral converges quickly. It is noted
that for constant forcing functions, and if the time of interest is

*
greater than the partition time, t1, the first integral on the right

'hand side of 12 becomes a constant and needs to be calculated only once.

14.6 One Dimensional CANSS Flowchart/Example

135

This section contains the flowchart/example for the one dimen-
sional CANSS program. The example problem is a finite body, insulated
on the right hand side, zero as the initial temperature, and zero as the
volume energy heat source. The left hand boundary has a linear heat
flux of the form q - qot. The symbol <CR> means carriage return.

Section 1 of Appendix C shows the input and output from running

this example.

0 Di ion CANSS Flowcha t

Flowchart Enter Explanation
Start <"canss.prg"<CR> Begin program while in SMP
Program

The subroutines exp.int,

Load grab.int and the external
subroutines CANSS integration routines are
loaded.
Display <CR> A banner is displayed that
Environment defines the CANSS environment

Load the left and right

Load <CR> boundary input conditions. The
B.C.,I.C.,& initial condition and heat
HG terms generation input routines are

loaded and note special values.

Input 2<CR> Input the left boundary
LB Qo<CR> condition(0-4), constants, and
Condition l<CR> the index on the time variable
Input 2<CR> Input the right boundary
RB 0<CR> condition(0-4), constants, and
Condition 0<CR> the index on the time variable
Input 0<CR> Input the initial condition.
Initial 0<CR> Type constants and polynomial

Condition power of coordinate x.

Input
Source

Display
Status

Load
GF
Library

Generate
GF

Calculate
small

boundary

solution

Calculate
large

boundary

solution

Calculate
initial
condition

solution

Calculate
volume
source

solution

STOP

4.7 Summary

0<CR>

0<CR>

<CR>

<CR>

<CR>

<CR>

<CR>

136

Input the volume energy source
term as a constant.

Display the description of the
problem.

Load the Green's function
library.

Calculate the Green's function
based on boundary conditions
for small and large times.

Calculate and display the small

time boundary solution.

Calculate and display the large
time boundary solution.

Calculate and display the
initial condition solution.

Calculate and display the
volume energy source solution.

End calculation and display
calculation time in CPU sec.

A one dimensional program for the symbolic solution of transient

heat diffusion problems is presented in this chapter. Three example

problems have been examined and checked against known solutions when

possible.

Some integrals that occur during the calculation of the

temperature distribution have been examined. The importance of time

137

partitioning for one dimensional problems is discussed. A flowchart of

the one dimensional CANSS algorithm is presented.

CHAPTER 5

TRANSIENT TWO DIMENSIONAL CANSSZD PROGRAHL

5.1 Introduction

The computer algebraic, numeric, and symbolic system called
CANSSZD calculates temperature distributions for transient, two dimen—
sional heat diffusion problems using a Green's function approach. The
program CANSSZD is written in the language of SMP [1983] and generates a
symbolic solution for the temperature distribution in a homogeneous
plate having boundary conditions of the zeroth, first or second kind on
any surface. A distinctive feature of the two dimensional CANSSZD
program is the ability of the program to allow nonzero, but constant,
boundary conditions to cover only part of a surface. The remaining
portion of the boundary condition at the surface is set to zero.

Special forms of nonlinear problems can be addressed for certain com-
binations of the physical parameters as was shown in Chapter 4.

The Green's function approach to the solution of the temperature
distribution for linear, transient heat diffusion problems in two dimen-
sions leads to integrals that have not previously been discovered or
integrals that are not well known. These integrals are discussed later
in this chapter or appear in Appendix B. Again, for certain combina-
tions of physical parameters, nonlinear problems may be solved for

temperature distribution.

 

 

139

Many of the integrals to be evaluated for two dimensional
problems are found in Appendix B of this thesis. Typically, the in-
tegrations of a one dimensional system provide integrals that are not
well known but have been studied extensively. The two dimensional
systems generate integrals that are not well known and have not been
studied as extensively.

Most temperature distributions for two dimensional heat diffu-
sion problems that appear in references and textbooks are left as a
function of an integral on time. Ozisik [1980], Carslaw and Jaeger
[1959] and other references generally avoid generating explicit, closed
form two dimensional temperature distributions and leave the distribu-
tions in terms of integrals. This is due to the complexity of the
evaluation of the integrals in closed form. The two dimensional CANSSZD
program returns functions that are recognizable and can be evaluated.
The CANSSZD two dimensional program solves the integrals and presents
the temperature distribution in terms of partitioned dimensionless time
so that the temperature distribution for each time region converge
quickly.

The Laplace transform technique, developed in Section 3.2,
produces Green's functions that are accurate and efficient at small
dimensionless times. Green's functions expressed in terms of Fourier
expansion and developed using the separation of variables technique for
finite bodies (see Churchill and Brown [1978]) can involve infinite
series that converge slowly at small times. In many cases, the tempera-
ture distributions obtained using the Green's functions developed by the
‘Laplace transform technique involve integrals that are unfamiliar, have
‘not been evaluated, or not been tabulated.

The Laplace transform approach is used in this thesis to deter-

Inine unknown integrals caused by the convolution of time. Arpachi

"II

 

140

[1966] describes the theory of convolution integrals in his text of
conduction heat transfer. Doetch [1961] relates the convolution in-
tegral to an effect (forcing function) multiplied by a weighting
function (Green's function). The individual integrands in the integral
are transformed into Laplace transform space, multiplied together, then
inverse Laplace transformed. Typically, all that is necessary for the
inverse Laplace transform is a short table in Appendix A.

A description of the CANSSZD algorithm for obtaining temperature
distribution in a plate and the flowchart for the CANSSZD algorithm are
presented in Section 5.2 of this chapter. Two example problems solved
using the CANSSZD program are discussed in Section 5.3. Some integrals
that arise from the calculation of temperature distribution in two
dimensional problems are reviewed in Section 5.4. Two and three dimen-
sional time partitioning is discussed in Section 5.5. Section 5.6

summarizes this chapter.

5.2 Transient Two Dimensional CANSSZD Program

The objective of this section is to present a program that will

calculate the temperature distribution in a two dimensional body. A

flowchart/example of the CANSSZD program appears at the end of this

section.

5.2.1 CANS Pro am

CANSSZD is a computer program designed to generate two dimen-

sional symbolic temperature distributions in a homogeneous plate using

 

141

symbolic manipulation. The temperature distributions symbolically
calculated by the algorithm are partitioned with respect to time. The
numerical evaluation of the temperature distribution expressions are
quickly convergent in each time region. No heat generation is allowed
in the plate and the initial temperature of the plate is constant at
zero. Only constant boundary conditions of the zeroth, first, or second
kind may occur at the surfaces of the plate. The describing partial
differential equation is given in equation (3.2.4) and the boundary
conditions are given in equations (2.2.11) and (2.2.13).

The non-homogeneous portion of the boundary condition may extend
to any percentage of the surface length but the boundary cannot have a
mixed condition. This means, for example, a boundary condition of the
first kind cannot coexist with a boundary condition of the second kind
on the same surface. The products of Green's functions that include
boundary conditions of the fourth and fifth kinds are not allowed in the
Green's function method for two dimensional bodies as was shown in
Chapter 2.

The temperature distributions are generated by the Green's
function approach described in Chapter 2, equation (2.3.15), where both
the initial condition and volume energy heat source are zero. The

equation for the temperature distribution is,

('1‘

S

 

fi(xi.yi.r)
T(x,y,t) - a I I k1 G(x,y,t|x',y',r)x,_x dSi dr.
T-O S1 i-l y'-yi
(for boundary conditions of the second kind)
t s'
- a I I E f.(x!,yi,r) fig dS. dr, (5.2.1)
. J J J anj x'-x
r-O Sj J-l y,_yg

142

(for boundary conditions of the first kind)

where s is the number of Robin conditions and s' is the number of
Dirichlet conditions on the boundaries of the plate. The boundary
condition of the zeroth kind is automatically included in equation
(5.2.1).

The CANSSZD program for obtaining temperature distributions of
spatially two dimension plate problems begins by loading five external
files. The names of procedures that accompany the SMP program begin
with capital letters and are underlined. The names of procedures with a
lower case letters and underlined were written by this author.

The first two external files loaded help simplify expressions so
that the integration routines do not waste time expanding or searching
through expressions. The Repgxplg external file of SMP performs re-
placements on terms that involve exponentials and logarithms. The
grab.int file removes constants from expressions so that the integration
routines do not need to parse through full expressions.

The next external file loaded by the CANSSZD routine is called
ex 'nt. It contains the procedures used to integrate the forcing
functions over space. The fourth file loaded in the CANSSZD algorithm
is the temporal integration procedure called z2d.int. This file con-
tains procedures to calculate the integration with respect to the time
variable. When an integral is not recognised by the external integra-
tion routines over space or time, the internal SMP integrator is
invoked.

Some constants are loaded along with some mathematical rules
that define properties of the transcendental functions. A banner is

displayed defining the environment for which the program is run.

 

mflrt“?“" ‘1

143

The first step in calculating the temperature distribution in
the plate is to describe the boundary conditions. The CANSSZD algorithm
can use boundary conditions of the zeroth, first and second kinds that
have a forcing function which is either zero or a constant. The input

data routine also asks for the length of the plate in the x and y direc-

L
tion and calculates the aspect ratio (r+ - L!) for the plate. The
x

aspect ratio is important in describing the partition times for each of
the three time regions and will be discussed in a following section.
The input routine also queries the user for the forcing function placed
on each boundary based on the type of boundary condition that occurs on
the surface. The CANSSZD routine has all the information necessary to
describe the temperature distribution for the stated problem.

The CANSS2D program begins the calculation of temperature dis-
tribution by deciding if the two dimensional problem is actually a
special case of a one dimensional problem. This may occur when a bound-
ary condition is semi-infinite. A special routine is loaded to handle
problems with boundary conditions of the zeroth kinds. This problem may
not be strictly two dimensional depending on the boundary conditions on
the other surfaces.

If the problem is not a special case, the Green's functions data
base is loaded and the appropriate Green's function for the three time
regions are calculated and displayed. The next step is to integrate
each of the Green's function with respect to the spatial variables and
display the results for the three time regions. The pattern that is
used in the integration is 1) check the integral against the external
CANSSZD integration library and if a match is not found, 2) let the SMP
internal integrator operate on the integral and if a solution is not

obtained, 3) return the integral to the user for further evaluation.

144

The CANSSZD program will generate closed form solutions for cases as-
sociated with the limitations placed on the boundary and initial
conditions, and the volume energy heat source.

The last calculation step is to generate the temperature dis-
tribution for the appropriate time region by integrating with respect to
time. The temperature distribution is displayed and the program stops
and displays the calculation time in CPU seconds. A flowchart of the
CANSSZD is shown in the next section and includes the input to the
example problem of a partially heated plate examined in the next sec-

tion.

5.2.2 Tw ens o a S owcha t Exam 1e

This section contains a flowchart/example of a two dimensional

plate partially heated on one side. The input and output from the

CANSSZD program is shown in Section 2 of Appendix C.

F wc art Enter Explanation
Start <"cansde.prg"<CR>
Program
Load The subroutines for the Green's

Subroutines function, grab.int, exp.int and the
external integration routines are
loaded.

Display <CR> A banner is displayed that defines the

Banner environment for 2-D problems.

Display <CR> Give the default values and define

Defaults some constants.

Enter
Boundary
Condition

Numbers

Display
GF Number

Input
Lengths

Input
Bottom
Forcing
Function

Input
Left
Forcing
Function

Input
Top
Forcing
Function

Input

Right
Forcing
Function

Load
2-D
Routine

Display
Small-Small
GF

Display
Small-Large
GF

Display
Large-Large
GF

Calculate

2<CR>
2<CR>
l<CR>
l<CR>

l<CR>
2<CR>

y<CR>

n<CR>
Qo<CR>
0<CR>
.5<CR>

y<CR>

y<CR>

<CR>

145

Enter the types of boundary condition
for the bottom, left, top, and right
sides.

Display the Green's function number
for this problem.

Input the length of the plate in the
x and y direction.

Input the forcing function on the
bottom of the plate. For zero, type y

Input the forcing function on the
left side of the plate. For zero,

type y.

Input the forcing function on the
top of the plate. For zero, type y.

Input the forcing function on the
right side of the plate. For zero,

type y.

Check for BC of zeroth kind. Load
appropriate subroutine.

Display the two dimensional Green's
function for the small-small time
region.

Display the two dimensional Green's
function for the small-large time
region.

Display the two dimensional Green's
function for the large-large time
region.

Calculate the spatial integration for

146

Spatial <CR> the three time regions.
Integration <CR>
Calculate <CR> Calculate the temporal integration for
Temporal <CR> the three time regions.
Integration <CR>
STOP Stop program and display the

calculation time in CPU seconds.

 

 

5.3 Two Dimensional CARSSZD Examples

Two example problems of the two dimensional CANSSZD algorithm
are presented in this section. The first example problem is a plate
with an aspect ratio of 1/2. A constant heat flux and temperature occur
on two of the four boundaries (X21B10Y21B01T0); see Figure 5.1a. The
second example problem is a a plate insulated on the bottom and with
zero temperature on the top and right hand side. The left boundary has
a constant heat flux over half the boundary and insulated otherwise

(X22B(l,0)0Y22T0), see Figure 5.lb.

5.3.1 A Two Dimensional Plate with Heating

The first example problem solved using the CANSS2D program is a
plate with a constant heat flux, qo, on the left surface, a constant
temperature, To, on the upper surface, insulated on the bottom and
having the temperature on the right hand side zero; see Figure 5.1a.

The aspect ratio of the plate is 1/2 or the length of the plate in the y
direction is L and the length of the plate in the x direction is 2L.

1}“; terms x and y, which represent distances in the x and y direction,

147

 

 

 

 

///////////////

Figure 5.1a Thin Plate With Heat Flux and
Temperature Conditions

[

y = 2L

 

\\\

l l:\\
T.

YEL—

 

 

 

 

/////

Figure 5.lb Partially Heated Thin Plate

 

 

148

are made dimensionless in the CANSSZD algorithm by dividing by the
respective lengths of the plate in the x and y direction. The Green's
function for this problem (X21Y21) is split into three convergent ex-
pressions for three regions of dimensionless time.

The third expression of the temperature distribution, the part
applicable when the dimensionless times for the x and y direction are
considered large, is examined first because an analytical solution is
available. Beck [1984a, pg. 1242] has obtained an analytical solution

for the temperature distribution to be,

 

 

co an 0 fl x fl y
T(XIYIt) - a E E [ 1 ' e-at[ ] ] (308(fi ) COS(—II':‘-) o
m—l n-l
n 5 2 q L
$31) 2 {4 To -° (~1)m - B 0k }, (5.3.1)
(3m + 4 fin) m n

where the term flm is an eigenvalue and m is the summation index for the
Green's function in the x direction, the term fin is the eigenvalue and n
is the summation index for the Green's function in the y direction. The
term a is the thermal diffusivity and k is the thermal conductivity of
the plate. The term [o], when the length of the plate in the x direc-

tion is 2L and the length in the y direction is L, is given by,

[.1 -[§%]2+[%] (5.3.2)

*
The CANSSZD algorithm splits the dimensionless time t into
*
three regions. In the third region, t2 is defined as the dimensionless
*
partition time between the second and third time region and t is the

dimensionless time of interest. A description of the partition times

 

1m‘r.‘ - _‘1 ‘fl

 

149

for two and three dimensional heat diffusion problems is given later in
this chapter.

The thermal conductivity and the thermal diffusivity default to
one in the CANSSZD program and give the change in temperature distribu-

tion in the third time region as,

” ° fl qL fix fly
AT(x,y,t) - 8 E E { 2 To En (-l)In - Eg— } cos(Emi) cos(‘E‘) .
m n

m-l n-l

 

' e

3492 + 4 22)
n -
('1) [ e 2 m “ (5.3.3)

* 2 2
m<g+459
2 2
mm+ag>
*

If the dimensionless partition time t2 is set to zero, the
temperature distribution from Beck and from the CANSSZD program match
exactly.

In the second time region, the change in temperature distribu-
tion that is efficient for times greater than the dimensionless

* *
partition time t1 and less than the dimensionless partition time t2 is

given as,

m+n+l B y

qoL a m (-D cost?)
AT(X.y.t)- T E E 2 A

n-l m--o fin

 

I3 B '4fifl
{ [ erfc(23n(tf)l/2 -(;$93/2> - erfc(Zfln(t*)l/2 -(;$?I/2)] e “ m°
1

 

 

fl 3 453
+ [ erfc(23n(cf)1/2 + m° ) - erfc(2fin(t*)1/2 +( m: )] e “ m° }
C

(cf)1/2 * 1/2

 

150

 

m m B Y
To (_1)m+n+1 cos(_1£_) *
+ .5 E E B [ coshe(t ’fln’flm+)
n-l m--w n
* * *
+ coshe(t ’fin’fim-) - coshe(t1,fln,flm-) - coshe(t1,fln,flm+) ].

(5.3.4)

where the term fin - x (n - 0.5) is the eigenvalue in the y direction.
The terms fimo - 0.5 (2 m + xL), fim+ - 0.5 (2 m + xL + 1), flm- - 0.5 (2 m

+X

L - l), and,

41919
coshe(t,fln.fimi) _ l; [ e 29 mi erfc[ 3 (t )1/2 2.1/2]
fin (t)
25 5 +' * fl 2 fl
- g r13. m L] , fit [L]
2 erf°[ fin<t ) +(t*>1/2 e n erfc (6‘)”2 1
(5.3.5)

The term coshe(-,-,-) will be examined more closely in the next section.

The term qo is the constant heat flux, L is the length of the
plate in the y direction and k is the thermal diffusivity. The most
dominant terms in the summations are m - l and n - 0. The term fin is
the eigenvalue for the large time Green's function in the y direction.
The term flmi can be thought of as similar to an eigenvalue for the small
time Green's function in the x direction. The term 3m: is not an eigen-
Value.

When the dimensionless time in both the x and y direction is

small the temperature in the first time region is,

q L co co
T(x.y.t) - 2—ok— E E (-1)m+“-

m--m n--w

151

 

 

 

 

 

 

2 ~1-
* 1/2 -fi /t 5
LE-L-- e mo [ erf( ‘2: ) ' erf( .2; ) ]
{ J; <c*)1/2 (6‘)“2
a 22+4f 52+4;
+ -E2 [ expi( l _mg__;___ni ) - expi( 1, mo * n- ) ]
t t
fimo fimo fimo fimo
* fimo [ ”((t*)1/2' 2 p ) “((t* 1/2' 2 B + ) ]

To °° °° fl 219 19 2fi
+_ §§(1)m+n[a(_29_’_119) H(_mo’ ‘10)],
a m—-¢ n--m (t*)1/2 Bm+ (t*)1/2 5m‘

(5.3.6)

where the function H(-,-) is discussed by Litkouhi [1982] and in
Appendix B. The symbols fino - 0.5 (2 n - yL + l), fin+ - 0.5 (2 n + yL +
1), fln- - 0.5 (2 n + yL - l), and the term expi(-,-) is the exponential
integral. The exponential integral is defined in Abramowitz and Stegun

[1964] as,

-Xt

 

expi[n,x] - En(x) - I e
1

n dt. (5.3.7)

(1

Other symbols in the temperature distribution equation, when the
dimensionless times for the x and y direction are considered small, have
been defined previously. The convergent equations for the change in
temperature distribution in the first and second time region have not

been previously investigated.

5.3.2 A Partially Heated Plate

152

The second example of the CANSSZD program, see Figure 5.1b, is
for a plate with zero temperature on the top and right hand side, insu-
lated on the bottom, a heat flux, qo, applied to the lower half of the
left hand side and insulated on the top half of the left hand side. The
aspect ratio for this case is two.

The temperature distribution in the plate for the first time
region generated by the CANSSZD program (when the dimensionless time in

the x and y directions is small) is,

qL m w
T(x,y,t) - 53; E ) <-1>”+“

m--w nF-m

 

2
* 1/2 -a /c 23 26 _
{ iE—1-—- e m° [ erf( -;9f72) - erf( 'IEI72)]
JR (t) (t)
fl 2+4fi2 fl2+4fi2
+ _E_ [ expi( 1, mo * n+ ) - expi( 1’_E2__;___B; ) ]
t t
l9 B l9 l9
+ 5 [ H( -mQ——-, -—mQ- - H( —99——-, m° ) ] . . (5.3.8)
m° (t*)1/2 2 fin- (331/2 2 fin+

The terms 5n+ - 0.5 (2 n + yL + %) and 3n_ - 0.5 (2 n + yL - %) for this
problem and flmo has been defined in the previous example.

The change in temperature distribution for the second time
region generated by the CANSS2D program, where one dimensionless time is
large and one is small, is,

5 x

.J1_)
L [ coshe(t*,flm,2 fin-)

 

AT(x,y,t) -

NIH

qoL m a (-l)n cos(
7))
m-l n--m fl

2
m

153

+ coshe(t*,flm,2 fin+) - coshe(tt,fim,2 fln+) - coshe(tf,fim,fln_) ],

(5.3.9)

where the term fin is the eigenvalue for the x direction, the terms fin:
were described in the previous time region and coshe(t*,fim,flni) will be
described in the next section.

The change in temperature distribution for the third time region

generated by the CANSSZD program, when both dimensionless times are

 

large, is,
... ... fix fly fl
qoL cos(—%—) cos(ifli) sin(zn)
AT(x,y,t) - 16 -k_ E 2 2 2 .
mpl npl fin (4 fim + fin)

- e

* 2 2 * 2 2
[ -t2(fim + fin/h) -t (am + fin/a) ]
e .

(5.3.10)

The change in temperature distribution when the time in the x
and y direction are considered large is shown in Table 5.1 for the

position x - % and a partition time t2 - 0.8.

5.4 Two Dimensional Integrals

Three integrals that often appear in calculating the temperature
distribution in two dimensional diffusion problems are presented in this
section. The first integral examined occurs in the calculation of
temperature distribution in the first time region. The next integral
occurs in the calculation of the temperature distribution in second time

region. The two integrals are difficult to solve and evaluate because

154

Table 5.1 Dimensionless Change in Temperature From Equation
(5.3.10) When the Partition Time Is Set to 0.1 and
the x Coordinate Is Set to L/2.

Position time - 1.0 time - 1.5 time - 2.0
y
0.0 3.482508-4 6.68336E-4 7.36753E-4
0.1 3.47147E-4 6.66239E-4 7.34445E-4
0.2 3.438463-4 6.59963E-4 7.27537E-4
0.3 3.38373E-4 6.49553E-4 7.16078E-4
0.4 3.30769E-4 6.35082E-4 7.001488-4
0.5 3.21092E-4 6.16653E-4 6.79859E-4
0.6 3.09412E-4 5.94393E-4 6.55349E-4
0.7 2.95815E-4 5.68458E-4 6.26788E-4
0.8 2.803993-4 5.39023E-4 5.94368E-4
0.9 2.632733-4 5.06288E-4 5.58306E-4
1.0 2.445548-4 4.70470E-4 5.18841E-4
1.1 2.24369E-4 4.31803E-4 4.76228E-4
1.2 2.02851E-4 3.90537E-4 4.30743E-4
1.3 1.80138E-4 3.46933E-4 3.826728-4
1.4 1.56373E-4 3.012623-4 3.32315E-4
1.5 1.31703E-4 2.53807E-4 2.79982E-4
1.6 1.06275E-h 2.048558-4 2.25991E-4
1.7 8.02387E-5 1.54700E-4 1.70667E-4
1.8 5.37464E-5 1.03638E-4 1.14338E-4
1.9 2.694948-5 5.19707E-5 5.73369E-5
2.0 0.00000 0.00000 0.00000

155

of the types of functions involved and the appearance of the convolution
of time. Other special integrals which may occur in two dimensional

problems are given in Appendix B.
5.4.1 A e ra me e o ne
An integral that often occurs when the dimensionless time in the

x and y direction is small and the coordinates xL and yL are normalized

with respect to the lengths Lx and Ly is,

 

erfc( ) dr (5.4.1)
r-O J (t") Jha(t-r)/Ly2

t 2 40 L 2
1 'J g- L/< (c-r)/ x > Y1.
1
Dimensionless groups used to eliminate the convolution on time

are defined to be,

+
x x y r
w _ L X _ L ’ Y _ L

J4a(t-r)/Lx2 JAac/Lx2 Jaat/Lx2

+ . . .
where (r ) is the aspect ratio. The integral I1 becomes,

 

w 2
“ e.w Y
I1 - 2/t X I 2 erfc( _i_ w ) dw , (5.4.2)
w
w-X

where t is the time of interest. Litkouhi [1982, pg. 123] has shown

that this integral can be written as,

 

 

156

_ { _ _x2 Y 2 2
I1 - 2 J: Jr ierfc(X) - e erf(Y) - _ E,(x + Y )

 

+ J; x H( Y,-§— }, (5.4.3)

where H(-,-) is an integral discussed by Litkouhi and obtained from a

text by Rosser [1948]. When X - 0, integral I1 becomes,

1, - 2 J? { erfc(Y) - -€E—- E,( Y2 ) }. (5.4.4)

1?

A plot of the integral when X - 0 is shown in Figure 5.2. This figure
shows the immediate effect on the temperature distribution of the in-
tegral at the surface Y - 0 as time increases. As the point of interest
is moved deeper in the body (i.e. Y - 0.25), the effect on the distribu-
tion is smaller and rises slower.

When Y - 0, integral 11 becomes,
1, - 2 j? { J; ierfc(X) }. (5.4.5)

A plot of this function is shown in Figure 5.3. This figure shows the
effect of the integral on the temperature distribution as the time
increases. Due to the nature of the ierfc function, a small increase in

the argument to the function leads to a large decrease in the value of

the function.

5.4.2 Two Integrals in Time Region Two

157

 

 

.0 I x co:3 oEwH
mmo~coMmcoEH0 msmuo> A¢.¢.nv COwumzvm cw douwoucH ~.n obsmwm

00.0 00.0 *00 No.0 00.0
h u h n b p _ b p u — p n - - ...-0.0

mN.0 n >

 

 

 

 

.‘

 

 

158

“I, LASHEI; . I

.0 n ... cos: 2:2.
mmoficofimcosfio mzmuo> 3.050 :owumscm cw Hmumouzn afim Housman

 

 

 

 

(a‘x'xflx

 

 

159

Two integrals that often occur in two dimensional problems when
the dimensionless time in one direction is small and in the other direc-

tion large are,

 

t*
2 c
- 2
I2 - I e 0,0 erfc( _ ] d0 , (5.4.6)
o-cf J9
and,
t*
2 2
I, - I 9‘1/2 e'(C1 9 + Cz/o) d0, (5.4.7)
o-cf

where 0 - Q—i%;11 is the convoluted time variable and C1 and C2 are
L

constants.

The solution to the integral 12 is found in Cho [1971] to be,

2 0 J9 J0
1
-C2 C2 t*
- 2 e 1 erfc( -— ) ] * . (5.4.8)

By defining a new function coshe(t,C1,Cz) to be,

- - C2
coshe(t,C1,C2) - 1 2 [ e 2 Cl C” erfC(C1Jt - 7; )
2 c1

 

02

C2 erfc(CIJE + 7:

2 c
- e2 C1 ) - 2 e'C1 t erfc( 7% ) ], (5.4.9)

160

the integral I2 may be written in compact form as,
* *
12 - coshe(t ,C1,C2) - coshe(t1,C1,C2). (5.4.10)

Two plots of the coshe function versus dimensionless time for
various values of the constants C1 and 02 are shown in Figures 5.4 and
5.5. Figure 5.4 plots the coshe function when the value of the
parameter C1 - g. This value corresponds to the first eigenvalue of a
Green's function for a X21 case. The difference between any two values
of dimensionless time is positive. The difference between two values of
the coshe function decreases as the value of the parameter, C2, in-
creases.

Figure 5.5 plots the coshe function when the value of the first
parameter C1 - x. This value corresponds to the first eigenvalue of a
Green's function for a X22 case. This figure shows the coshe function
decreases rapidly as the parameter C1 is increased. Additional eigen-

values are not necessary for computations involving the coshe function.

A summary of the coshe function is:

 

l) t + 0 coshe(0,C1,C2) 4 0
2) t 4 w coshe(m,C1,C2) + 0
3) C1 # 0
4) C1 + m coshe(t,w,C2) * 0
2
e-Clt
5) C2 4 0 coshe(t,C1,0) * - 2
C1

6) C2 4 m coshe(t,Cl,w) 4 0.

161

 

 

 

 

 

 

 

0.00..s,-..,-..,- ,.3
-0.05d
-O.10d
on
m2. -O.15-
.” .
u
75 -O.20~
.c:
m
8
-o.25-4
-o.:sod
-o.35....-.,...,. ,
0.0 0.2 0.4 0.6 0.8
* gfg-r)
t - 2
L

 

1.0

Figure 5. 4 The coshe Function Versus Dimensionless Time When

C1- 2 and Various Values of 6,.

 

 

 

 

 

 

0.00 . r
J 0.5 1
‘0.01" -'
IN 0.25
0? 4 d
k
* .
3 -O.02'-l q
2
8 4 .
U
-0003“ q
-o.04 - . ,
0.0 0.8 1.0

Figure 5. 5 The coshe Function Versus Dimensionless Time When

6,- x and Various Values of C2.

 

162

The solution of the next integral, 13. which appears regularly
in the second time region can be found in Abramowitz and Stegun [1964,

pg. 304] to be,

c c t

J« 2C1C2 .3. '2C1C2 —3 1
13 - 2 C; e erfc(C1/0 + Jo) + e erfC(2C1/9 + J9) 6-t*
(5.4.11)

The solutions to both integrals, I, and 13. include a term that

may cause numerical instability of the total solution. The term,

C2
ezclc2 erfc(C1/9 + 7;), (5.4.12)

may not be easy to evaluate due to the positive argument in the exponen-
tial function. The complementary error function can be approximated for

large values of the argument ( > 1.7 ) and is,

02 2 2
erfc(Cljo + J6) fr e e

 

 

 

1 l 3 ]
- + . ... ’
[ (C1/9 + C2/j0) 2(01/9 + C2//0)3 4(C1/9 + C2/J9)5

(5.4.13)

No approximation is necessary when the argument to the com-

plimentary error function is small ( < 1.7 ).

163

5.5 Time Partitioning in.Two and.Three Dimensions

The time partitioning scheme for two or three dimensional heat
diffusion problems is more complicated than for a one dimensional heat
diffusion problem because the dimensionless time variables for the
additional coordinates are not equivalent unless the lengths of the
sides of the plate are equal.

The CANSSZD program for a two dimensional plate bases the dimen-
sionless time on the x-coordinate length of the plate. This choice is
arbitrary and need not be followed in general. The generalized dimen-

sionless time for the plate is,

c - 2 . (5.5.1)

In the y-direction, the dimensionless time based on the length in the y

direction is,

t - 2 . (5.5.2)

Defining an aspect ratio (r+) as,
(r+) - L /L (5 5 3)
y X, . .
the generalized dimensionless time in the y-direction can be written as,

*
t -t (r) . (55.4)

164

The time region is split into three components. When the aspect
ratio (r+) is greater than one, see Figure 5.6a, the first time region
is defined as the XSYs region because the Green's function for this
area, in both the x and y direction, are for small time. The dimension-
less time goes from zero to tlx' The second time region, XIYS, is when
the Green's function in the x-direction is for large times and the
Green's function in the y-direction is for small times. The time region
for region 2 goes from tlx to t2x' The third time region is valid for
times greater than tgx and is designated XlYg. Both the Green's func-
tions, in the x and y direction, for the third region are based on the
large time.

When the aspect ratio is less than one, see Figure 5.6b, a
change in the Green's function occurs in the second time region. The
Green's function in the x-direction is for small times and in the y-
direction is for large times, or XSYI.

* *
Once a partition time is chosen, the values of t

1x and t2x may
be calculated. These values depend on the aspect ratio and the boundary
condition applied to the plate. When the aspect ratio is greater than

one, see Figure 5.6a,

* *

tlx - c, (5.5.5)
and,

* * +2

t2x - t1 0 (r ) , (5.5.6)

*
where t1 is the dimensionless partition time. The partition time is
defined as the time when the small time Green's function can be ex-

pressed with less than four reflection terms.

165

l

 

 

 

 

I
l I
Temperature '
Small x Large x ' Large x
Small y | Small y | Large y
l l
I
‘ l
l
l I
l l
o l l -
....
O tlx-t1 tZX-t1 r
Figure 5.6a Time Partitioning for an Aspect Ratio
Greater Than One for a Two Dimensional Body
l , ,
Temperature I |
l
Small x I Small x | Large x
Small y I Large y Large y
l
I I
l I
l
I
l
I I
o _I + I -__
0 tlx-t1 r t2x=t1

Figure 5.6b Time Partitioning for an Aspect Ratio
Less Than One for a Two Dimensional Body

166

When the aspect ratio is less than one, see Figure 5.6b,

-k 9: 2

tlx - t, - (r+) , (5.5.7)
and,

* *

t2x - t1. (5.5.8)

When the aspect ratio is equal to one, i.e. the plate is square,

* * *
tlx - t2x - t1, (5.5.9)

and only two time regions appear.

The time regions for the three dimensional case follow in the
same manner as the two dimensional case. For three dimensional cases,
it is convenient to make the assumptions that the direction of the
shortest length of the body is in the x direction, the next shortest
length is in the y direction, and the longest length is in the z direc-
tion. Assigning these coordinate directions leads to the following

three relationships for the aspect ratios,

£1 +

L - (r )yx > 1. (5.5.10)
x

Lz +

L— - (r )zx > 1, and (5.5.11)
x

L

.2._ +

L (r )zy > 1. (5.5.12)

167

If, for example, the partition time in the x direction is
defined to be 0.1, and the lengths of the body in the x, y and z direc-
tion respectively are Lx - 1, L.y - 2, and L2 - 4, the aspect ratios for

+ + +
the previous equations are (r )yx - 2, (r )zx - 4, and (r )z - 2.

Y
Choosing the x direction as the characteristic dimension yields,

2 - 0.1 , (5.5.13)

as the definition of the partition time. When the time in the charac-

teristic direction is 0.1, the associated time in the y and z direction

2 2
L L
- - - 1
QL%_L1 _ el%_zl _g _ a 2 T -§ - 0.1 (Z) - 0.025, (5.5.14)
L L L L L
y y x x y
and,
L2 L2
21%;11 _ 5-1 _§ - 21%:Ll -32 - 0,1 (‘%g) - 0.00625 (5.5.15)
L L L L L
2 Z X X Z

This means that when the time based on the characteristic direction x

forces the program to switch the x direction Green's function from small

to large time, the local time in the other two directions are still
small and the small time Green's functions are appropriate. A descrip-
tion of this example is shown in Figure 5.7. Notice that no more than
four time regions may appear for the three dimensional problems. The
number of time regions is reduced by one if the length of two of the
sides of a three dimensional body are equal. If the three dimensional

body is square, only two time regions appear.

 

“...-any-

 

168

x

o N o XN 0
~ A. +L use H A. +L a A. L

x»
+

mgmg: atom _a:o_m:me_a «mesh a not m:o_mom :o_uwugaa mew» ~.m mtam_m

 

 

 

 

 

 

".0 . mNo.0 mN00.0
.III- A 1/\7/ w “
«.0 “.0 m~0.0
0.~ v.0 ~.o
. _ _ _
_ _ .
. _ .
.
_ _
_ _ _
. _
_
_ _
_ _ _
N omens ~ __osm _ ~ __asm _ N __oEm
» omens . a mason _ » __osm . » __aEm
x mate; _ x mates _ .x muses _ x __oEm ocauotanmh
_ _ .

 

169

5.6 Summary

A program that successfully generates symbolic temperature
distributions for plate geometries for three types of boundary condi-
tions is presented. Two example problems were chosen and calculated by
the CANSSZD algorithm. Some parts of the two example problems appear in
the literature and the CANSSZD program matches the analytical solutions
exactly. Some integrals that appear in the calculation of the tempera-
ture distribution of a two dimensional plate are discussed extensively.
The concept of time partitioning in two or three dimensions is dis-

cussed.

CHAPTER 6

SUMMARX’AND CONCLUSIONS

6.1 Sun-ary

A computer algebra program called SMP is applied to transient heat
diffusion problems using a Green's function approach to obtain analyti-
cal temperature distributions in one and two dimensions. The method
restricts the describing partial differential equation, initial condi-
tion, boundary conditions, and volume energy heat source to be linear.
The symbolic algorithm is written in the SMP programming language for
the cartesian coordinate system.

The Green's function approach is an ideal candidate for computer
algebra programs. The formalism for the approach, developed in Chapter
2, relies on the complex mathematical concepts of differentiation and
integration. The CANSS programs applies the symbolic integration and
differentiation procedures that are internal to SMP and that have been
developed for this thesis to the Green's function formalism and calcu-
lates symbolic temperatures. The computer aids the calculation by
repeated application of the Green's function and the integration
knowledge base to the formalism described in Chapter 2.

The Green's function approach operates on transient, linear, multi-
dimensional heat diffusion equations and is developed in Chapter 2. The

boundary conditions are linear and include the traditional boundary

170

171

conditions of the first, second and third kinds. Three additional

boundary condition types are discussed in this chapter. The three
addflflonal types of boundary conditions are a boundary condition of the
zeroth'kind, or a "natural" condition, a boundary condition of the
fourth kind, a non-convective thin film condition, and a boundary condi-
tion of the fifth kind, a convective thin film.

A major difficulty of the Green's function approach is obtaining
Green's functions for small dimensionless times that converge quickly.
A new analytical technique for generating Green's functions efficient
for small times is given in Chapter 3. An equation and a group of
tables are presented in Chapter 3 to generate the Laplace transformed
small time Green's functions. A small table of Laplace transforms
(Appendix A) is used to transform these functions back into real time.

The Green's functions that are efficient for large times are
developed by the use of the separation of variables technique and have
been investigated in many texts and papers for boundary conditions of
the first, second and third kinds.

The effect of different types of boundary conditions on the surface
of a semi-infinite body in one dimension is described in Chapter 3.
Boundary conditions of the zeroth, first and second kinds are examined
and boundary conditions of the third, fourth, and fifth kinds with
«different values of their associated properties are examined. The

figures show that when the location of the source and point of interest

< the Green's

coixuzide and.the dimensionless time is small (t*_ 0.05),
function is the same for all boundary conditions.

.A one dimensional program called CANSS which uses the techniques of
the Green's function method to generate symbolic temperature distribu-
tions for a variety of boundary conditions, initial conditions and

'volinme energy heat sources is presented in Chapter Four. Three one

172

dimensional example problems are examined. The first example problem is
well documented and the CANSS distribution matches the known solution
exactly. The second and third example problems have not been examined
previously and represent new solutions. Two types of integrals that
appear during the calculation of the one dimensional temperature dis-
tribution, called the Dawson and convoluted exponential integral are
examined. The efficient evaluation of these integrals are described and
shown in the figures.

Sixteen distinct cases of the geometry can be accessed by the CANSS
program. Boundary conditions that are functions of polynomials in time
(i.e. l, t, t2...), initial conditions that are functions of the spatial
coordinates (i.e. l, x, x2...), and volume energy heat sources that are
polynomials in both time and space can be split up by superposition,
calculated independently, and added to obtain temperature distributions
because of the linearity of the Green's function approach.

Time partitioning of the one dimensional problems is necessary for
the efficient evaluation of the temperature distribution. The dis-
cussion and figures in Chapter Four describe the method to be used for
partitioning the time variable and suggests the point of time for the
partition. A flowchart for the one dimensional CANSS routine is
presented.

Symbolic temperature distributions for problems in two dimensions
using the Green's function approach are presented in Chapter Five. A
symbolic algorithm called CANSSZD calculates the temperature distribu-
tion boundary conditions of the zeroth, first, and second kind. The
solutions to many two dimensional problem in the literature and texts
are often left in terms of the initial functions, but the CANSSZD
program calculates the symbolic solutions in terms of the coordinates

and an appropriate range of time. There are a minimum of ninety three

173

distinct cases for the CANSSZD program. This number is lower than
actually possible because the forcing function on the boundaries may
extend over any portion of the surface which would lead to an infinite
number of cases handled.

Two example problems are presented in Chapter Five. The symbolic
temperature distribution for a portion of the first example can be found
in the literature. The temperature distribution output from the CANSSZD
program matches the known distribution exactly. New expressions for the

efficient calculation of temperature distribution are calculated for the
small and medium time ranges. These solutions have never appeared in
the literature.

The second example problem in Chapter 5 shows the strength of the

CANSSZD program. Exact temperature distributions are generated for a
plate with only partial heating occurring at one surface. None of the
solutions for this case have appreaed in the literature.

Two integrals that often appear in calculating the solutions to two
dimensional problems are presented. Each integral is described and
examined to obtain solutions that are efficient.

Time partitioning in two and three dimensions is explained and the
technique used in the CANSSZD program is shown in some figures. Some

insight to the extension of the time partitioning method is given for

three dimensional problems.

6 .2 Conclusions

The use of a symbolic manipulation program to generate a knowledge

based expert system is a new development for a computer algebra system.

174

The SMP language was used to develop a knowledge base of Green's func-
tions and integrals to be used in conjunction with two programs called
CANSS and CANSSZD. A numbering system associated with the Green's
functions [Beck and Litkouhi, 1983] has been used extensively by the
CANSS program based on the Green's function formalism. The numbering
system eliminates the need in the CANSS program for calculation of the
appropriate Green's function based on the problem geometry and boundary
conditions The Green's function data base contains all the possible
combinations of one dimensional cartesian Green's functions. The
strength of the CANNS program is not in the controlling program, but in
the knowledge bases over which it has control.

A symbolic manipulation program operates well in an problem en-
vironment which includes a rigid structure. The mathematical formalism
of the Green's function approach is particularly appropriate for use in
a computer algebra system. The symbolic manipulation program takes
advantage of the structure by eliminating the need for calculating the
appropriate Green's function from the homogeneous Green's function
equation. Instead, a data base of Green's functions is used to deter-
mine the Green's function. Also, a small data base of integrals is
included that keeps the program from searching through extensive exter-
nal libraries, wasting memory space and computer time.

It is important to note here that the smaller the data base, the
faster the calculation of the temperature distribution. In most computer
algebra systems, any defined function, such as a Green's function, is
loaded into the memory and retained until the system is halted or the
definition is removed. This extra baggage has a detrimental effect on
the speed of the calculations because larger pieces of memory will be
swapped to the storage device. The number of page faults (the number of

times the memory is filled and must be dumped to a storage device) can

175

become enormous if a check is not kept on the data base size. For

example, in the CANSS environment, the Green's function data base is

split into three libraries. Each library can be loaded into the CANSS

program but since all the definitions are not necessary, a large amount

of dead storage is not carried with the program and swapped in and out

of memory core .

Computer algebra programs can also be useful for calculations that

do not include a broad class of problems. SMP was use to calculated and

collect coefficients from the partial fraction expansion. These cal-

culations are straightforward, but are very repetitious and tedious when

performed by pencil-and-paper. The probability of human error in these

computations is high when done by hand but low when delegated to a

computer algebra program.
The advantage of using SMP over other computer algebra software is

the small amount of memory initially allocated loading the program. For

small and simple problems, many users may interact with SMP simul-

taneously. The kernel of functions in SMP is small when compared with

other computer algebra systems. External files, which contain special

functions and procedures written in the SMP language, are only loaded
when necessary. This causes the SMP basic kernel of functions to grow

by adding to the memory and storage, thus slowing down the calculation

speed.

The program SMP is written in the C language which makes it port-

able to many hardware systems. Other computer algebra systems are

written in the computer language LISP and are machine and hardware

dependent.

SMP has a very poor user interface. The natural form of equations

is two dimensional, but SMP can input only strings of text. This weak-

ness in the user interface is common to all computer algebra systems and

176

nmch work is necessary to strengthen the area. Some recommendations
regarding the computer algebra user interface are given later in this
chapter.

Most computer algebra software developers claim a simple yet effec-
tive language for their system. For many simple problems, the
interactive nature of the software is direct and efficient. When it is
necessary to program complicated problems in mathematical physics, the
interactive language may retard the progress of the solution to the
problem. This is due to the lack of a debug function that could catch
major flaws in the input of the procedures and functions.

Good program documentation is essential for developing functions
and procedures in new programming languages. Documentation of the
capabilities of computer algebra systems and, most importantly, limita-
tions of the program are poorly described. The documentation available
for computer algebra systems typically consists of a library of func-
tions available to the user. A few vague examples of simple functions
are examined in the reference guide and for simple or small problems,
this type of documentation is sufficient.

Examples of well documented procedures written in the SMP language
cannot be found in the SMP reference guide or primer. The SMP reference
and primer direct the user to poorly documented external files in the
SMP library for instruction. This may be compared to learning French
armed with only a dictionary and a simple set of phrases. The structure
and intent of a command may be different from what is expected even
though the words may be correct.

The use of computer algebra software for the analysis of problems
in heat transfer is appealing because of the ability of the software to

obtain symbolic solutions. Numerical solutions to heat transfer

177

pmoblems are acceptable as long as the solution at a particular coor-
dinate and time is desired. Even so, numerical instabilities of the
method may cause the numerical solutions to deviate or, in a worse case,
to oscillate from the correct solution.

The computer algebra program CANSS can be useful identifing areas
in conduction heat transfer problems that cannot be solved in closed
form. Identification of these areas can lead to detailed examination of

the types of integrals which must be evaluated to obtain symbolic solu-
tions.

The new concept of time partitioning is used to optimize the sym-
bolic temperature distributions calculated by the two CANSS programs.
These symbolic solutions are quickly convergent for the time regions of
interest. This is the first application of the time partitioning scheme
introduced by Keltner and Beck [1985] for symbolic temperature distribu-
tions.

Expert systems using computer algebra systems are best suited for
fields that exhibit narrow, specialized domains. Almost every technical
area can benefit from the use of computer algebra based on this descrip-
tion. For example, the field of acoustics is a narrow, specialized
domain of a broad class of wave mechanics. Sommerfield [1949, Section
27, Appendix II] shows the Green's function approach can be applied to
problems in acoustics and results in a formalism that is similar to the
formalism for heat transfer described in this thesis. The Green's func-
tions for the wave equations can be cataloged and stored in a knowledge
base. Special integrals that occur during the calculation of the solu-
tions may also be stored in the knowledge base. The formalism and the

knowledge base can lead to an expert system for wave equation problems

in acoustics.

178

A floppy disk containing the CANSS and CANDDZD programs can be
obtained from the author. All libraries, procedures and sub-procedures

can be loaded to a VAX/VHS microcomputer and run with SMP.

6.3 Recs-lendations

l. The small time Green's function method should be applied to
other orthogonal coordinate systems such as the radial, spherical, and
elliptical systems. The approximations in the radial coordinate system
are similar to the approximations made in the cartesian coordinate
system except the functions that are approximated are different. Care
must be exercised in the application of the Green's function approach to
other coordinate systems. In the radial coordinate system for example,
only a limited number of one, two, and three dimensional cases
can use the Green's function formalism because of the differences in the
partial differential equations between the rectangular and radial coor-
dinate systems. A data base of small and large time Green's functions
need to be calculated for the radial coordinate system. Beck [1986] is
continuing to catalog small and large time Green's functions for rectan-

gular, radial, and spherical geometries for various boundary conditions.

2. The kernel or influence functions developed in the small
time Green's function approach should be applied in unsteady surface
element [Keltner and Beck, 1981] and finite element methods. In the
unsteady surface method, small time Green's functions are particularly
important because they allow calculations at small and large times are
made speedily. The small time influence functions developed can be used

as accurate trial functions in finite element methods which will improve

179

the convergence of the solution. Green's function based influence
functions can be used for trial functions for non-linear problems. The
efficiency of both numerical procedures will be enhanced through the use

of functions that are defined for specific ranges of time.

3. The integration procedures used in SMP should be extended to
recognize and include more functions. The work of Cherry [1985] and
Knowles [1986] in the area of symbolic integrators should be examined
closely and incorporated into the internal SMP integration routines.
Their work extends the algorithm of Risch [1969] to include specified

logarithmic and exponential functions.

4. The two dimensional CANSS program should be extended to
include additional types of boundary conditions, initial conditions, and
volume energy heat sources. An additional boundary condition that can
be readily applied to the CANSSZD program is the boundary condition of
the third kind. Initial conditions that are polynomial and transcenden-
tal function can be investigated. Volume energy heat sources that are

functions of the spatial coordinates and time can be examined.

5. The one dimensional CANSS program should be extended to
include boundary conditions of the fifth kinds, convective thin films.
Due to the non-symbolic nature of the Green's functions associated with
boundary conditions of the fifth kinds, more variables will be gener-
ated, which will cause a shortage of storage space, and greatly slow

down the calculation of the temperature distribution.

6. The Green's function approach to the solution of heat diffu-

sion problems is not serial but parallel. The formalism and linearity

 

 

 

 

 

 

180

of the method allows each piece of the solution to be calculated inde-
pendently. Each processor in the parallel computer could be assigned a
portion of the problem and independently and simultaneously execute the
operations necessary for the solution. The SMP program has anticipated
the move into parallel processing by including functions that takes

advantage of the parallelism.

7. A major flaw in computer algebra systems is the inability of
the program to interface with the user. Equations are two dimensional
objects with structure but computer algebra system treat them as strings
of text. A user interface in needed for both the input and output that
can make use of bit-mapped video screens, menus, and a mouse. Bit-
mapped displays allow the user to point and pick objects such as
equations and portions of equations. Also, bit-mapped screens can draw
the special symbols that appear in mathematical equations. Draw down
menus can speed the input of complicated expression by drawing mathe-
matical symbols, along with the mathematical operation, to the input

line.

8. Output expression display should make use of windows or
horizontal scrolling. The expressions generated by a typical computer
algebra problem are, in many cases, longer than the 80 characters avail-
able on a normal video screen. Windows would allow the user to split
the input or output expressions into small portions. Horizontal scroll-
ing would be an alternative method of displaying large output

expressions.

9. An on-line status area is recommended for computer algebra

systems. The on-line status area can inform the user of the complexity

181

of a calculation by indicating the amount of time in the calculation,
the amount of storage space and memory used by the calculation, and an

indication of the progress of the calculation.

10. When an integral in the CANSS environment is not evaluated,
it may represent a new type of function. Numerical integration could
completely describe the unknown integral. The new expression for the
unknown could be entered into the integration procedures of the CANNS

environment and the solution could be determined.

11. A data base of convolution type integrals associated with
diffusion problems is recommended. The Risch-base integration proce-
dures cannot evaluated these integral at the present time. Computer
algebra systems can use their expertise in handling Laplace transforms

and inverse Laplace transforms to begin this data base.

12. The output of symbols and equations from the computer al-
gebra systems should be in a form that is natural to the investigator.
The computer algebra output should be able to generate symbols such as
integral signs, summation signs, and partial derivitive symbols to name
a few. More work using the graphical capabilities of bit-mapped screens

is necessary for the natural output of symbolic expressions.

APPENDICES

APPENDIX.A

LAPLACE TRANSFORHS FOR SHALL‘TIHES

A short table of Laplace transforms is given below. See Abranowitz

and Steguns [1959, pp. 1021-1026] for a comprehensive table of Laplace

 

 

 

 

 

 

transforms.
f(s) F(t)
l
s 1 (A.1)
1: t (A.2)
s
—(—)-P V (v > 0) tV-l (A.3)
sl/
1 -at
s + a e (A.4)
l at
s _ a e (A.5)
ls l 2 a Daw
- (a ft) (A 6)
S + a? r(1/2) t1/2 r(1/2)
l 2 a
Daw(a ft) (A.7)
/s (s + a2) F(l/Z)

182

-a/s

-a/s

e

Js

 

-a/s

 

e
./s (b + J5)

e-a‘lS3
s (b + J8)

 

 

183

i-az/(a c)

g
2 2

r(1/2) t3/

E-a2/<4 t)

r(1/2) c1/2

2
-b t + a b a
e erfc( 7(Z_E) + b/t ]

 

l a
erfc[ J(4 t) ]

2
-b t + a b
e

OWH‘

erfc( 7(%_E) + b/t ]

(A.8)

(A.9)

(A.10)

(A.ll)

APPENDIX.B

SOME USEFUL INTEGRALS

B.1 Introduction

The purpose of this Appendix is to collect and present special
integrals important in diffusion problems. Many of the integrals found
in this appendix are used in the CANSS programs. A definition that

appears in following sections is,

-(x-x')2
K(x-xv,t-r) - K(-x+x',t-r) _ (4 I a (t T))-1/2 e4a(t - T)

(3.1)

184

185

B.2 Integration Over Space with Respect to x' for Small Time Functions

B.2.

Int E uat'on

b

I K(x-x',t-r) F(x') dx' - I K(-x+x',t-r) F(x') dx'

b
a
F x'
6(xo-x')
l
a;
L

a

Integral

K(x-xo,t-r), a < xo < b, otherwise zero

(B.2)

 

 

NIH

erfc[ x - b ] - erfc[ x - Q45 ] }
{ (40 (ti-THU2 (4a (t-r))1/2

(B.3)

 

 

_x_ erfc( x ' b ] - erfc[ x ' é ] }
2 L (4a (c-r)>1/2 (40 <t-r>)1/2

+ 2.94% { K(X’a,t'7) - K(X'b,t‘7) }
L

(B.4)

 

erfc[ x - b ] - erfc[
{ (4a (t-r))1/2

 

l}
(4a (t-r))1/2

+ 23—L%;11 { (x + a)K(x-a,t-T) - (X + b)K(X'b't’7) }

L

(B.5)

186

B.2.2 Integral Eguation

 

 

 

b b
I K(x+x',t-r) F(x') dx' - I K(-x-x',t-r) F(x') dx'
a a
F x' Integral
6(xo-x') K(x+xo,t-r), a < xo < b, otherwise zero
(B.6)
Ll {-J—L ] [ l}
1 erfc - erfc
2 (4a (t-mL/2 (40 mm”2
(B.7)
x' x x + b x + a
-— ‘-— erfc[ ] - erfc( ' ] }
L L L (4a (t-r))L/2 (4a (c-¢))L/2

+ 2L.(_12=;T_l{1((x+a,t:-r) - K(x+b,t-T) }
L

(B.8)

 

 

{ erfc[ x +él 2] - erfc( X + b 1 2] }
(4a (t-r)) / (4a (t-r)) /

+ Lg-L%;Ll { (x - b)K(x+b,t-r) - (X ' a)K(X+a't'T) }
L

(B.9)

187

B.3 Exponential and Error Integrals

B.3.l Inte x n

 

n Ingggzél
-3 (4 «)1/2 [ 7i— ierfc(C Ja) - 7%— ierfc(C Jb) ]
(3.10)
-l (7r/C2)1/2 [ erfc(C Ja) - erfc(C lb) ]
(8.11)
2 2
o -L; [ e‘c a - e'C b ] (3.12)
C
2 a
l [ -1% e'C 0 + -1% erfc(C J0) ]
C 2C 9=b
(3.13)
3/2 3/2 2 a
3 [ [ ;_27_ + i‘; ] e.C 0 + 21% erfc(C J0) ]
2 C C 4C 9-b

(B.l4)

B.3.

 

2

188

Inte r res on
b

2
I on/2 e-c /9 do

a

Integral

(1r/C2)1/2 [ erfc(C/Jb) - erfc(C/Ja) ]
(B.15)

(4 «)1/2 [ (b)L/2 ierfc(C jb)
- (0)1/2 ierfc(C Ja) ] (3.16)

2 b
[ 03/2 e’c /9 + 2 02 J(«o) ierfc(C/JO) ]

umv

0-a

(B.l7)

2
1%. [ 03/2 (3o - 4(32)e'C /9 + 4 C4 /(«0) ierfc(C//9> ]

(B.18)

b

0=a

 

189

B.3.3 Inte ress’on
b 2 2
I on/2 e-c1 o -c,/o do
a
n Integral
-3 11‘ e201C2 erfc(C J0 + 23)
202 1 J0
C2 b
- e-2c,c, erfc(C1/0 - 73) ] (3.19)
0-a
C
A. 2C1C2 _2
-1 2C1 e erfc(C1/9 + J0
C2 b
+ e-2C1C2 erfc(C1/0 - 73) ] (B.20)
0-a
0 13‘ ezclczerfc(C 0 + SE)
4C1 1 0
C2 b
+ eO2C1C2 erfc(Clfl - 5‘) ] (B.21)
6-a
11_ 1 2c 0 C2
1 2 "' - C2 e L 2 erfc(C1/0 + -)
2C 261 Je
1
C2
+ [ 5%: - c2] 8-20102 erfc(C1/6 - 73)
2 2 a
+ 2 ff 6'01” ‘ 02/0 ] (3.22)

6-b

190

30 C
2 2 CC 2
IE; [ [ C2 + -;3 - EE- ] e2 1 2 erfc(C1/0 + 7;)
2C1 4C1 2
BC C
2 3 2 -2C C 2
+ [ C2 + “—3 - ——— ] e 1 2 erfc(Cljfi - “)
401 2C1 J6

a

2 2
.12[_:L+Hcl]e-clo 'C2/9]
0-b
(8.23)

191
8.3.4 Inte ra Ex ression

b

n
I 0 erfc(C 9) d0
a

n I e a
b b 2 2
-1 [ £n(6) erfc(co) ] + 29 I £n(0) e'C 9 do
O-a "
a
(3.24)
o % [ ierfc(C a) - ierfc(C b) ] (3.25)
1 b
1 ——3 [ erf(C 0) - 2 0 C ierfc(C 0) ]
AC 0-a

(8.26)

192

8.3.5 Int a e ion

b
I a“ erfc(C/J0) do
a

.___n___ te a1

b
[ erfc(C/la) - 2 c/o ierfc(C/JG) ]

0-a

(8.27)
2 2
[ 2; erfc(C/J0) + 29311 ierfc(C/Jfi)

c 3 1/2 -cz/0 b
- 3 (6 /n) e ]

0-a

(8.28)

193

B 3.6 .13223131_£zn;2§§123

b
I an” erfc(Cljo + 02/10) do
a

 

___£L__. .13£2312L_
0 [ ‘1; [ erfc(C1J9 + C2/Jfl) + e-4C1C2
401
ii . a
erfc(C1/0 - Cz/Jfl) ] + C 1erfc(C1/0 - Cz/Jfl) ]
1 6-b
(3.29)
2 2 1 -(c2 - c2)o
(c1 - c3) 2 2 [ 2 e 1 3 erfc(Cljo + cz/jo) do
2(c1 - c3)
0
- (1 + 51) e‘2(c1 ' C3)C? erfc(C3/0 + cz/Jo) d9
2
c b

+ (1 - 5%) e’2(C1 + C3)C2 erfc(C3/6 - cz/Ja) do ]
6-a

(8.30)

194

8.4 Integration Over Time for Some Error Function Integrals

An integral which occurs often in two dimensional problems is,

I - 51; I w'3 erfc[x w] erfc[y w] dw (3.31)
w-1/(4ac)1/2

Integrating this equation by parts, letting,

v - erfc[x w] erfc[y w]

du - w-3dw
yields,

3.[ -2 J”
I - w erfc[x w] erfc[y w]

a “ 1/(4at)1/2

z -2 2
_ _ w exp[-(xw) ] erfc[y w] dw
2“/“w-1/(aac)1/2
_ ——1:_ v-2 exp[-(yw)2] erfc[x w] dw (8.32)
Za/w w-l/(Aat)1/2

Note that the first and second integrals in the above equation has the
same form but different parameters.
Substuting u - x w in the first integral of equation (8.32),

the integral can be written as,

KN

m
2
II - —:5 x2 I u‘2 e'u erfc[u/p] du (3.33)
" x

195

where X2 - x2/(4at), Y2 - y2/(4at), and p - x/y or X/Y.

Litkouhi [1982] shows that this integral can be expressed as,

- 2
II - 2:E x2 [ 1% ierfc[X] - e’x erf[Y]/X
In
_1_ 2 2 .. -<p w)2
- _ 81(X + Y ) + 2p I e erf[w] dw ] (8.34)
p/« x

The second integral in equation (8.32) can also be placed in the form of

equation (8.34) as,

2

III - 2:3 Y2 [ 1% ierfc[Y] - e'Y erf[X]/Y
/«
2 2 2 2 m -(w/p)2
- _ 81(X + Y ) + E I e erf[w] dw } (8.35)
1l'
Y

Litkouhi replaces the integrals in equations (8.34) and (8.35)
with a function called H(-,-) which is examined in a text by Rosser
[1948]. The solution to the following four integrals below use a

similar technique and functions.

196

t
[erf( —3—] erf[ —'Y-—] dr -
-0 J4 a (t-r) J4 a (t-r)

1.

2 l egficfixz erfc(!) erfc(X) ggfich)
t - - +
2 2 2 2
[ X ierfc(X) + Y ierfc(Y) ]

- [x zerfc(X,Y) + Y zerfc(Y,X) ]//Z (3.36)

t
I erfC[ :] erfc[ 121—'4] dr -
-0 J4 a (t-r) J4 a (t-r)

1'

2 t erfc(X) erfc(Y)
2

- [ X zerfc(X,Y) + Y zerfc(Y,X) ]/J; (8.37)

 

Q
2
where X - _§_ , Y'- ':E:- , zerfc(X,Y) - X I u'1 e.u erfc(xig)du
J4at J4at X
— -x2 Y 2 2 -
or zerfc(X,Y) - Jn ierfc(X) - e erf(Y) - : EI(X + Y ) + Jw X H(Y,X/Y)
«

and where H[-,-] is the integral discussed in Litkouhi [1982] and Rosser

[1943].

197

E
Jr erf[ t] erf[ _—1—__] dr -
-o 14 a (t-r) 14 a (t-r)

T

2 4 4
2 t a: - 8 (i erfc(X) + i erfc(Y) ) + ——L)———(—lerf° X “erfc Y

2 2
(1 +§ ) zerfc(X,Y) - (1 +1; ) zerfc(Y,X)

=III><
§llh<

J
2 2 2 2 2 2
+272 [(11 +Y)3,(x +Y)-e'(X+Y)]

2 2

+ -—l- [ X e.x erf(Y) + Y e.Y erf(X) ]

(8.38)

 

2 e fc er c Y
2 t 4
2 2
- é (1+%)zerfc(X,Y) - E (1+§)zerfc(Y,X)
Jr Jr
2 2
2 2 2 2 _
+ ¥[(X+Y)E,(X+Y)-e(x+y)]
2 2
+ 1_ [ X e.X erf(Y) + Y e.Y erf(X) ]
6 Jr

(8.39)

where X, Y, zerfc(X,Y), and H(X,Y/X) have been defined on the previous

page.

198

8.5 Some Integrals for Small Tile Green's Functions

C
I r“ K(x,t-r) dr - 3%; r( g + 1 ) (4 c)(“+1)/2 in+lerfc[ --i¥i—-— ]

 

 

 

 

r O (4at)1/2
n - 1-,o,1,2,... (3.40)
t n 3:2
I r“ 315311 K(x,t-r) dr - 2—32 r( g + 1 ) c 2 .
L L
r-O
1/2 | I | I
{ [ 22% ] in+3erfc[ x 1 2 J + in+2erfc[ x 1 2 ] }
x (4at) / (4at) /
n--2,-l,0,l,2,... (3.41)
t $2
I r“ erfc[ lxl 1/2 ] - r( 3 + 1 ) (4c) 2 in+2erfc[ -—-i§if7§ ]
r-O (4a(t-r)) (4at)
n--l,O,l,2,... (3.42)
t
1/2
I K(w,t-r) dr - [ :t ] ierfc[ lwl 1/2 ]
r-t-At (40 At)
(3.43)
t t 1/2 W2
2 I a(t-r) K(w,t-r) dr - I [ 2L£:Ll ] e4a(t-r) d7

r-t-At 13C-At

199

 

 

 

 

 

 

3 1/2 w
-—w—[erfc[—v Hm—H [1- ...],
6 a (40 At)1/2 « w w
(8.44)
t
I erfc( w ] dr - 4 At 12erfc[ jg ]
r_t_At (4a(t-r))1/2 (4a At) 1/2
(8.45)
t
I (t-r) erfc[ Y ] dr -
r-t-At (4a(t"))1/2
2 2 w 4 w
4(At) { i erfc[ ] - 4 i erfc[ ] (8.46)
l/2 (4a At)1/2

(4a At)

2
t W

I [ a w 1/2 ] e4a(t-r) dr - erfc[ -——J3d-——']

3
T_t_At [4n (a(t-r)) 1 (4a At)1/2

 

 

(8.47)

AI2IIDIX C
CAISS AID CAIBSZD PIDGIAI llalPLlS
Welcome to Paul Zang's version of SMP.
SMP 1.5.0
28-FEB~1987 13:54:49.05

qullz: <‘canss.prg'

The CANSS procedure takes about 30 seconds to load!!!!
Please wait a moment ......

The [zang.smp.c2d]grab.int function is loaded.
The [zang smp.c2d]exp.int library is loaded.
The CANSS utilities are loaded.

ccccccccc AA N N ssssss ssssss
c AA rm N s s
c A A u u N s s
c AAAAA u n 1: ss 53
c A A u an s s by
ccccccccc A A u N 3555555 sssssss p. H. Zang

This procedure will calculate the temperature distribution in a one
dimensional slab that has boundary conditions of the first through
fourth kind at any surface. Constant heat generation may occur
and the initial temperature of the slab is a polynomial function of x.

Press (Enter) or <Return> key to continue:

The CANSS left boundary input routines are loaded.
The CANSS right boundary input routines are loaded.
The CANSS IC & HGT and Status are loaded.

The CANSS function to generate CF‘s is loaded.

NOTE
Small dimensionless time is defined as being < 0 025
Assume the forcing functions on the surfaces begin at tau - O

200

201

Press <Enter> or <8eturn> key to continue:

Left Boundary Condition

What type of boundary condition do you have:

- Semi-infinite Condition
Temperature Condition

- Heat Flux Condition

- Convective
Non-convective Thin Film

user—O
I

>
I

Enter : 2

Form: dT/dx - constant*time‘(n/2) n - -l,0,l.2,...
Enter the constant value.Qo

Enter the value for n. n can equal -l.0.l,2.... 1

Right Boundary Condition

What type of boundary condition do you have:

- Semi-infinite Condition
- Temperature Condition
Heat Flux Condition

- Convective

- Non-convective Thin Film

bearer-IO
0

Enter : 2

Form: dT/dx - constant*time‘(n/2) n - -l.0.l.2....
Enter the constant value 0

Enter the value for n.O

Initial Condition

202

The initial condition term is constant with respect to time.

Form :: F(x) - Constant * x“$n , where Sn - positive integer.
Enter the constant value for the initial condition.0

Enter the value for $n.0

Enter the constant value for the heat generation.0

Status of the CANSS routine at this time.
You are dealing with the X( 2 2 ) case.
The initial condition is equal to 0
The heat generation is constant and equal to 0

Loading the [zang.smp.gflibs]eigen.con library.
Hang on.... This takes about 10 seconds!!
Press <Enter> or <Return> key to continue:
Loading the [zang smp.gflibs]gfld.lib library.
Hang on.... This takes about 30 seconds!!
The [zang.smp.gflibs]gfld.lib library is loaded!!
The small time Green's function is ::
2
«0.25 (2nn + x - xp) -O.25 (2nn + x + xp)
0.5(Exp[ --------------------- ] + Expl --------------------- ])
theta theta

theta Lx Pi °
The large time Green's function is ::

l + 2Cos[mm x Pi] Cos[mm xp Pi] Exp[- mm theta Pi ]

oooooooooooooooooooooooooooooooooooooooooooooooooooo

Press <Enter> or <Return> key to continue:

Integral e 3

The integral of 0 is not in the library.
we will try the internal integrator.

The small time boundary solution is 2:

2 0.5(2nl + x)
4tl Lx Qo Gamma[1.5] IErfC[ ------------ .2]
0.5

203

0.5
alpha k
NOTEzz n1 or n2 are summation indexes that may go from
minus infinity to plus infinity.

Press <Enter> or <Return> key to continue:

0.5
The integral of tau is not in the library.
We will try the internal integrator.
The integral of 0 is not in the library.

we will try the internal integrator.
The large time boundary solution is ::

0.5
-DawsonsInt[ml Pi (-tl + t2) 1

2 2
* Exp[- ml t1 Pi ]
2Cos[ml x Pi] ( --------------------------------
ml Pi

2 2 0.5
2 + Epr- ml tl Pi ] (-tl + t2) )
Lx Qo ( ------------------------------------------------------
2 2
ml Pi

1.5
. 0.666667 t1 + 0.666667 t2 )

0.5
alpha k
NOTEz: ml or m2 are summation indexes that may go from 1 to infinity.

Press <Enter> or <Return> key to continue:

The integral of 0 is not in the library.
We will try the internal integrator.
The initial condition solution is 2: 0

NOTE:: ml or m2 are summation indexes that may go from 1 to infinity.

Press <Enter> or <Return> key to continue:

The integral of 0 is not in the library.
we will try the internal integrator.
The small time solution for the heat generation term is: 0

NOTE:: nl or n2 are summation indexes that may go from
minus infinity to plus infinity.

Press <Enter> or <Return> key to continue:
Err[6l.35.0]

31(2]:: <end>

2(14

welcome to Paul Zang's version of SHP.

SMP 1.5.0
28-FEB-l987 13:59:01.05

eI[l]:: <'canss2d.prg'

Loading the XRepexplg external file...

Loading the [zang.smp.c2d]grab.int external file...
The {zang.smp.c2d]grab.int function is loaded.
Loading the [zang.smp.c2d]exp.int external file...
The [zang.smp.c2d]exp.int library is loaded.
Loading the [zang.smp.c2d122d.int external file...
Loading the [zang.smp.z]de.int library

The [zang.smp.c2d122d.int library is loaded.

ccccccccc AA N N ssssss 555555
c A A NN N s s
c A A N N N s s
c AAAAA N N N ss 55
c A A N NN s s by
ccccccccc A A N N 5555555 sssssss 2. N. Zang

This procedure will calculate the temperature distribution in a two
dimensional plate that has boundary conditions of the first or
second kind at any surface. No heat generation will occur in the plate
and the initial temperature of the plate is zero. You may either use a
constant forcing function at the surfaces or zero at the surfaces. The
forcing functions may extend to any percentage of surface length.

Press <Enter> or (Return) key to continue:

Some values must be given before we can continue. The default values
of the constants are:

thermal conductivity (k) - 1

thermal diffusivity (alpha) - 1

NOTE
Small dimensionless time is defined as being < 0.025
Assume the forcing functions on the surfaces begin at tau - 0

Do you wish to use the default values? (Y or N)
Enter : y

 

 

205

Now we determine the boundary conditions. The options are ::
- Semi-infinite boundary
1 - Temperature Condition (Dirichlet)
2 - Heat Plux Condition (Neumann)

"What type of condition occurs on the bottom? (O.l.2)'2
"What type of condition occurs on the left side? (0.1.2)'2
"What type of condition occurs on the top? (0,1,2)'l

”What type of condition occurs on the right side? (0.1.2)'l

 

This is the X( 2 1 )Y( 2 l ) case

This is the input data routine. You will be asked
to enter the length of the plate in the x and y
direction. Then you will be asked to input the
forcing function that occurs at the boundaries.
Please continue .....

 

What is the length of the plate in the xodirection?l
What is the length of the plate in the y-direction?2

Is the forcing function on the bottom zero?
Enter : y

Is the forcing function on the left side zero?

Enter : n
What is the value of the heat flux forcing function
on the left side of the plate?

Forcing function on this surface - Qo
Where does it start? Dimensionlessly (0 - l) - 0
Where does it end? Dimensionlessly (0 - l) - .5

Is the forcing function on the top zero?
Enter : y

Is the forcing function on the right side zero?
Enter : y

Loading the [zsng.smp.gflibs]gf.lib library.
The [zang.smp.gflibs]gf.lib library is loaded!!
4*e******eeeee*************************************e***********e**

mink **
** You have chosen a two-dimensional problem **
** **

w*****e**teete*****************************************e**********
We begin by setting up the 2-D Green's function for time region 1.

2
nl + n2 -0.25 (2n1 + x)
0-5 ('1) Q0 Exp[ ---------------- ]
theta

206

2 2
- (2n2 + y - yp) - (2n2 + y + yp)
* (Expl ----------------- l + Exvl ----------------- 1)
theta theta
(0. -------------------------------------------------------------- ,0.0)
theta Pi

The GP for time region 1 is completed

Begin setup of GP for time region 2.
n2 2 2
(~l) Qo Cos[x Pi (-0.5 + ml)] Expl-theta Pi (-0.5 + ml) ]

2 2
- (2n2 + y - yp) - (2n2 + y + yp)
* (EXPI ----------------- l + Earpl ----------------- ])
theta theta
(0, ------------------------------------------------------------ ,0.0)
0.5 0.5

theta Pi

The GP for time region 2 is completed

Begin setup of GP for time region 3.
{0,200 Cos[x Pi (-0.5 + ml)] Cos[0.5y P1 (-0.5 + m2)] Cos[0.5yp P1 (-0.5 +
m2)]

2 2 2 2
* Exp[-theta Pi {-0.5 + ml) ] Epr-O.25theta Pi (-O.S + m2) 1, 0.0)
The GP for time region 3 is completed

Now we have finished entering the data. We move to integration
over space...
Now I am calculating the spatial integration for region 1....

2
n1 + n2 ~0.25 (2nl + x)
0.25 (-l) 00 Exp[ ---------------- ]
theta
-0.5(l - 4n2 - 2y) 0.5(1 + 4n2 + 2y)
* (Erfc[ ------------------ ] - Erfc[ ----------------- l)
0.5 O 5
theta theta
0 5 0.5
theta Pi

Finished with region 1!!!

Now I am calculating the spatial integration for region 2....
n2 2 2
0.5 (-1) Q0 Cos[x Pi (-0.5 + ml)! Exp[-theta Pi {-0.5 + ml) 1

-0.5(1 - 4n2 — 2y) 0.5(1 t 4n2 r 7y)
* (Erfcl ------------------ ] - Erfc[ ----------------- })
0.5 0.5

 

2(17

theta theta
Finished with region 2!!!

Now I am calculating the spatial integration for region 3....

2 2
4Qo Cos[x 21 (-0.5 + .1)] Cos[0.5y 31 (-0.5 + -2)] Exp[-theta Pi (-0.5 + ml)
1

2 2
* Exp[-0.25theta Pi (-0.5 + m2) ] Sin[0.25Pi (-0.5 + m2)]

Pi (-0.5 + m2)
Finished with region 3!!!

Here is the temperature for time region 1 ...

n1 + n2 0.5 0.5
0.25 (-l) Q0 (2 cl Pi
2
-0.5(l - 4n2 - 2y) -0.25 (2nl + x)
oErfI ------------------ ] Exp[ ---------------- 1
0.5 tl
ti
e ( ..............................................
0.5
Pi
2
0.25 (2nl + x)
0.SExpi[l, ---------------
t1
(1 4n2
2
0.25 - 2y)
+ -------------- ]
tl

0.5
tl Pi
0.5(2nl + x) -(2nl + x)
0.5HH[ ------------ , ------------ 1
0.5 l - 4n2 - 2y
:1

* (2nl + x)
4. ................................

0.5

 

2CH3

0.5(2nl + x)
+ lErfcl ------------ l)
0.5
tl
0.5 0.5
- 2 tl Pi
2
0.5(1 + 4n2 + 2y) -0.25 (2nl + x)
-Erf[ ----------------- ] Exp[ ---------------- ]
0.5 (:1
tl
* ( .............................................
0.5
Pi
2
0.25 (2nl + x)
0.5Expi[l. ---------------
t1
(1 + 4n2
2
0 25 + 2y)
+ -------------- ]
tl
* (2nl + x)
0.5
ti Pi
0.5(2nl + x) 2n! + x
0.5HH{ ------------ . ------------ 1
0.5 l + 4n2 + 2y
tl
* (2nl + x)
+ ................................
0.5
tl
0.5(2nl + x)
+ IErfc[ ------------ )))
0.5
cl
0.5
Pi
More is the temperature for time region 2
n2
0.5 (-1) Q0 Cos[x Pi {-0.5 + m1)!
-0.5(l - 4n2 - 2y) 2 2
-(-Erfc[ ------------------ ] Expl-tl Pi (-0.5 1 ml) g

0.5
:1

.0.5(1 - 4n2 - 2y) 0.5

 

209

- 0.5Erfc[ ------------------ + :1 91 (.o.s + n1)]
0.5
:1

* Exp[-P1 (-0.5 + II) (1 - hn2 - 2y)]

0.5(1 - hnz - 2y) 0.5
+ 0.58rfc( ----------------- + :1 P1 (-0.5 + n1))
0.5
:1
* Exp[P1 (~O.5 + m1) (1 - én2 . 2y)])
* ( ...............................................................
2 2
P1 (-0.5 + .1)
-0.5(1 - 6n2 . 2y) 2 2
-Erfc[ ------------------ I Expl-tZ P1 (-O.5 + m1) ]
0.5
:2
~0.S(1 - an2 - 2y)
- 0.5Erfc[ ------------------
0.5
t2

0.5
+ :2 Pl (-0.5 + m1)]

* Exp[-P£ (-0.5 + m1)
* (1 - anZ - 2y)]

0.5(1 - an2 - 2y)
+ 0.5Erfc[ .................
0.5
t2

0.5
+ t2 P1 (-0.5 + m1)]

* Exp(P1 (-0.5 + m1) (1 . 6n2 - 2y)]

Pi (-0.5 9 m1)

0.5(1 + AnZ + 2y) 2 2
-Erfc[ ----------------- ] Epr-tl Pi (-O.S + m1) ]
0.5
t1

-O.S(1 + An2 + 2y)
+ 0.5Erfc( ------------------
0.5
:1

0.5
+ t1 Pi (-O.S + m1)]

* Exp[-Pl (-O.5 + m1)

211)

* (1 + anZ + 2y)]

0.5(1 + hn2 + 2y)
- 0.58rfcl -----------------
0.5
:1

0.5
+ :1 Pi (-0.5 + n1)]

* ExplPi (-0.5 + .1) (1 + hnZ + 2y)]

P1 (-0.5 + n1)

0.5(1 + hnZ + 2y) 2 2
~Erfc[ ----------------- ] Expl-tZ P1 (-0.5 + m1) ]
0.5
:2

-0.S(1 + an2 + 2y)
+ 0.SEr£c[ ------------------
0.5
:2

0.5
+ t2 Pi (.0.5 + m1)]

* Exp(-P1 (-0.5 + n1)
* (1 + hn2 + 2y)]

0.5(1 + AnZ + 2y)
- 0.SErfc[ -----------------
0.5
:2

0.5
+ t2 P1 (-0.5 + m1)]

* ExplPl {-0.5 + n1) (1 + 6n2 + 2y)]

P1 (-O.S + II)
Here is the temperature for time region 3 ...
aQo Cos[x P1 (-0.5 + n1)] Cos[0.5y P1 (-0.5 + n2)]

2 2 2 2
* (-Exp{t2 (- P1 (-0 s + m1) - 0.25 P1 (-0.5 + m2) )1

2 2 2 2
+ Exp[t3 (- 91 (-0.5 + n1) - 0.25 P1 (-0.5 + m2) )1)

* Sin[0.25P1 (00.5 + 32)]

P1 (-O.5 + n2) (- P1 (-0.5 + .1) - 0.25 P1 (-0.5 + m2) )
DDD....DDDDD ...... DDDDDDets’ Ill folks'!!!!!l

LIST 0? REFERENCES

LIST OF REFERENCES

Abramowitz, M. and Stegun, I. A., editors, 1964, Handbook of
Mathematical Eungtiong with Formulas, Graphs and Mathematical
Iables, National Bureau of Standards, Applied Mathematics Series
55, U.S. Printing Office, Washington, D.C.

Aizen, A. M., Redchits, I. S., and Fedotkin, I. M., 1974, "On Improving
the Convergence of Series Used In Solving the Heat Conduction
Equation," Journal of Engineering Physics, Vol. 26, pp. 453-458.

Arpachi, V. 8., 1966, Conduction Heat Iransfer, Addison-Wesley, Reading,
Mass.

Beck, J. V., 1986, ME817 Class Notes, Mechanical Engineering Department,
Michigan State University.

Beck, J. V. and Litkouhi, B., 1985, "Heat Conduction Numbering System
for Basic Geometries," ASME Winter Annual Meeting, November, Paper
No. 84-WA/HT-63.

Beck, J. V. and Keltner, N. R., 1985, "Green's Function Partitioning
Procedure Applied to Foil Heat Flux Gages," ASME 85-HT-56,
Presented to the National Heat Transfer Conference, Denver, CO,
August 4-7.

Beck, J. V., 1984a, "Green's Function Solutions for Transient Heat
Conduction Problems," Int. J. of Heat and Mass Transfer, V0127,
pp. 1235-1244.

Beck, J. V., 1984b, "Green's Function and Numbering System for
Transient Heat Conduction Problems," AIAA Journal, Vol 24, No. 2,
pp. 327-333.

Butkovskiy, A. C., 1982, Green's Eunction and Transfer Handbook, Wiley,
New York

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

211

212

Char, B. W., Geddes, K. 0., Gonnet, G. H., Watt, S.M., 1985, Mapie
User'g Guidg, WATCOM Publications Ltd., Waterloo, Ontario,
Canada.

Char, B. W., Geddes, K. O., Gonnet, G. H., Marshman, B. J., Ponzo, P.
J., 1986, "Computer Algebra in the Undergraduate Mathematics
Classroom,” Proceedings of the 1986 Symposium on Symbolic and
Algebraic Computation, July 21-23, University of Waterloo, Ontario,
Canada.

_Cherry, G. W., 1986, "Integration in Finite Terms With Special
Functions: The Logarithmic Integral," SIAM J. Comput., Vol 15, No.
1, pp. 1-21.

Chester, M., 1963, "Second Sound in Solids," Physical Review,Vol. 131,
No. 5, pp. 2013-2015.

Cho, 8. W., 1971, "A Short Table of Integrals Involving the Error
Function," unpublished, Department of Mechanics, Korean Military
Acadamy, Seoul, South Korea.

Churchhill, R. V., and Brown, J. W., 1978, Fourier Series apd Bopndary
Value Erpbiems, Third Edition, McGraw-Hill, N.Y.

Cole, K. D., 1986, "Conjugated Heat Transfer With the Unsteady Surface
Element Method," Ph.D. Thesis, Mechanical Engineering Department,
Michigan State University.

h 2
Dawson, H. C., 1898, "Numerical Value of I ex dx ," Proceedings of the
0

London Mathematical Society, Vol. 29, pg. 519.

Doetsch, G., 1961, Guide to the Application of Laplace Transforms, D.
Van Nostrand Company, London, 255 pp.

Greenberg, M. D., 1971, Application of Greep's Functions in Science and
Engineering, Prentice-Hall, New Jersey, 132 pgs.

Haji-Shiekh, A., and Lakshminarayanan, R., 1986, "Integral Solution of
Diffusion Equation with Boundary Conditions of Second and Third
Kinds," ASME Winter Annual Meeting, Paper No. 86-WA/HT-83, Anaheim,
California, December 7-12.

213

Hassanein, A.M. and Kulcinski, G. L., 1984, "Simulation of Rapid Heating
in Fusion Reactor First Walls Using the Green's Function Approach,"
Journal of Heat Transfer, Vol. 106, pp. 486-490.

Hayes, J. E., and Michie, M., 1984, intelligent Systems; rhg
Unprggeggptgd Qppprrppiry, J. Wiley and Sons, N.Y., N.Y.

Hayes-Roth, F., Waterman, D. A., Lenat, D. B., (eds.), 1983, Buiiding
firpgrr Systems, Addison-Wesley, Reading, Mass.

Hearn, A. C. (ed.), 1985, REDUCE User's Manuai, Version 3.2, The Rand
Corporation, Santa Monica, California.

Keltner, N. R., and Beck, J. V., 1981, "Unsteady Surface Element
Method," Journal of Heat Transfer, Vol. 103, pp. 759-764.

Knowles, F., H., 1986, "Integration of Liouvillian Functions with
Special Functions," Proc. 1986 Symposium on Symbolic and Algebraic
Computation, University of Waterloo, Ontario, Canada, July 21-23,
pp. 179-184.

Lebedev, N., N., 1965, S ec a u t'on and e'r A lications,
Prentice-Hall, Englewood Cliffs, New Jersey, 308 pp.

Litkouhi, B., 1982, "Surface Element Method in Transient Heat Conduction
Problems," Ph.D. Thesis, Mechanical Engineering Department,
Michigan State University.

The MATHLAB Group, 1983, MACSYMA Reference Manual, Version Ten,
Laboratory for Computer Science, Massachusetts Institute of
Technology, Cambridge, Massachusetts.

Maxwell, J. C., 1867, "On the Dynamic Theory of Gases," Phil. Trans.
Roy. Soc., Vol 157, pp. 49-88.

Mikhailov, M.D., and Ozisik, M.N., 1983, "Diffusion in Composite Layers
With Automatic Solution of the Eigenvalue Problem," Int. J. Heat
Mass Transfer, Vol 26, No. 8, pp. 1131-1141.

Miller, W. L. and Gordon, A. R., 1931, "Numerical Evaluation of

Infinite Series," J. Physical Chemistry, Vol. 35, No. 2, pp. 2785-
2884.

Moenck, R., 1977, "On Computing Closed Forms for Summations,"
Proceedings of the 1977 MACSYMA Users Conference, NASA CP-2012,
NASA Scientific and Technical Office, Washington, DC.

214

Morse, P. M. and Feshbach, H., 1953, Mgrnong pf Ibeoreticai Physics,
McGraw Hill, New York.

mu-Math, 1983, EHLEQSD Rgfprenng fining, Soft-Warehouse, Hawaii.

Ng, E. W., 1977, "Observations on Approximate Integrations,"
Proceedings of the 1977 MACSYMA Users Conference, NASA CP-2012,
NASA Scientific and Technical Office, Washington, DC.

Ozisik, M. N., 1980, Heat Conduption, Wiley, New York.

Ozisik, M. N., 1984, "Propagation and Reflection of Thermal Waves in a
Finite Medium," Int. J. Heat Mass Transfer, Vol 27, pp. 1845-1853.

Rand, R., 1984, Compnrgr Alggnrn in Applng Marhenntics, Pitman
Publishing, Marshfield, Mass.

Risch, R., 1969, "The Problem of Integration in Finite Terms," Trans.
Amer. Math Soc., Vol 139, pp. 167-189.

Ritt. J. F., 1948. W. Columbia University

Press

Roach, P. J. and Steinberg, S. , 1984, "Symbolic Manipulation and
Computational Fluid Dynamics," AIAA Journal, Vol. 22, pp. 1857-
1865.

Rosser, J.B., 1948, " eo a cat on o

O‘——;N
(D
Q.
X
n)
:3
Q.

2 2 2

Z Y
-p y 'X n
e dy e dx , Mapleton House, Brooklyn, N.Y.
O 0

Shibata, T. and Kugo, M., 1983, "Generalization and Application of
Laplace Transformation Formulas for Diffusion," Int. J. Heat Mass
Transfer, V01. 26, No. 7, pp. 1017-1027.

SMP, 1983, SMP Reference Guide, Inference Corp., California.

Sommerfeld, A., 1949, Pnrrini Differential Equations in Physics,
Academic Press, New York, 335 pgs. '

215

Stakgold, I., 1979, Qrggn's Eunntipng nnd Boundary Value Problems,
Wiley, New York.

Tzeng, L, and Beck, J. V., 1985, "Data Base for Solutions in Transient
Heat Conduction," M.S. Project, Michigan State University,
Department of Mechanical Engineering, MSU-ENGR-85-018.

Vernotte, P., 1958, "Les Paradoxes de la Theorie Continue de 1'equation
de la Chalear," Comptes Rendus, Vol. 246, pp. 3154-3155.

Vick, B. and Ozisik, M. N., 1984, "Growth and Decay of a Thermal Pulse
Predicted by the Hyperbolic Heat Conduction Equation," Journal of
Heat Transfer, Vol 105, pp. 902-907.

Walters, A. C., 1949, "Solution of the Transient Diffusion Equation by
Means of Green's Functions", Proc. Camb. Phil. Soc., Vol. 45, pp.
69-80.