n u..-

u.

...n..

 

 

, .THEscs

kg
‘3
0

()5

lllllllillllfllllihlWWill”|H|H|HH|1||HINHIIHN

3 1293 02074 799

This is to certify that the

thesis entitled

A PRIORI TIME STEP ESTIMATE FOR
SPHERICAL AND RADIAL FIELD PROBLEMS

presented by

Munawar Hussain Chaudry

has been accepted towards fulfillment
of the requirements for

Master's Jame in Mechanicai Eng

 

 

MW

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

 

LIBRARY
Michigan State
University

 

 

 

PLACE IN RETURN Box to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11m Mass-p.14

A PRIORI TIME STEP ESTIMATE FOR SPHERICAL AND RADIAL FIELD
PROBLEMS

By

Munawar Hussain Chaudry

A THESIS

Submitted to
Michigan State University
In partial fulfillment of the rquirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

1998

ABSTRACT
A PRIORI TIME STEP ESTIMATE FOR SPHERICAL AND RADIAL FIELD
PROBLEMS
By
Munawar Hussain Chaudry
The objective of this study was to develop an a priori time step estimate for three
single step methods used to solve the system of ordinary differential equations associated
with radial and spherical field problems.

The hypothesis was that the a priori time step estimate has the general form

(AUX, = C for the unconditionally stable methods and (A011,... = C for

conditionally stable methods, where C is a constant to be determined by numerical
experimentation, and l] and km” are the lowest and highest eigenvalues in the system of
ordinary differential equations.

Numerical solutions of step change problems were used to determine the
coefficient C for each solution procedure in time and each type of physical problem. The

final a priori time step equations developed in this study were

Central Difference Method: (ADM = 0.050 for N 2 11
Backward Difference Method: (ADM = 0.025 for N 2 11
Forward Difference Method: (Ammax = 1 for N 2 11

where N is number of nodes in space. Each equation can be used for both the spherical
and the radial problem.
The time step equations were validated by using different problems involving a

different set of material properties and boundary conditions.

ACKNOWLEDGEMENT

The author wishes to express his deepest gratitude to the continued
encouragement, guidance and support of Dr. Larry J. Segerlind. His trust as mentor is
greatly appreciated.

Appreciation is expressed to Dr. John B. Gerrish, Dr. Craig W. Somerton and
Dr. John J. McGrath for serving on the committee.

A great appreciation to all my friends at Michigan State University, Mechanical
as well as Agricultural Engineering Department.

Special thanks goes to my family in Pakistan, specifically my parents, my in-laws
and my children for their special support and encouragement.

Deepest gratitude is expressed to my wife, Nusrat for her love, her patience, her

support and her efforts in looking after our affairs. She has made my life meaningful.

iii

TABLE OF CONTENTS

LIST OF FIGURES ................................................................................ v
CHAPTER 1

INTRODUCTION ................................................................................. 1
CHAPTER 2

REVIEW OF THE LITRATURE ................................................................ 5
CHAPTER 3

OBJECTIVES ..................................................................................... 19
CHAPTER 4

THEORETICAL CONSIDERATIONS ....................................................... 21
CHAPTER 5

METHODOLOGY ............................................................................... 3 1
CHAPTER 6

RESULTS: UNIFORM GRID .................................................................. 36
CHAPTER 7

RESULTS: EQUAL VOLUME GRID ....................................................... 59
CAHPTER 8

EVALUATION OF THE TIME STEP ESTIMATES ...................................... 78
CHAPTER 9

DISCUSSION AND CONCLUSION ......................................................... 84

LIST OF FIGURES

FIGURE 5.1

SAMPLING POINTS IN SPACE ............................................................... 34
FIGURE 6.1

SPHERICAL UNIFORM GRID 6 NODE CD ................................................ 38
FIGURE 6.2

SPHERICAL UNIFORM GRID 11 NODE CD ................................................ 40
FIGURE 6.3

SPHERICAL UNIFORM GRID 21 NODE CD ................................................ 41
FIGURE 6.4

SPHERICAL UNIFORM GRID 6 NODE BD ................................................. 42
FIGURE 6.5

SPHERICAL UNIFORM GRID 11 NODE BD ................................................ 43
FIGURE 6.6

SPHERICAL UNIFORM GRID 21 NODE BD ................................................ 44
FIGURE 6.7

SPHERICAL UNIFORM GRID 6 NODE FD .................................................. 46
FIGURE 6.8

SPHERICAL UNIFORM GRID 11 NODE FD ................................................. 47
FIGURE 6.9

SPHERICAL UNIFORM GRID 11 NODE FD ................................................. 48
FIGURE 6.10

RADIAL UNIFORM GRID 6 NODE CD ......................................................... 50
FIGURE 6.11

RADIAL UNIFORM GRID llAND 21 NODE CD ........................................... 51
FIGURE 6.12

RADIAL UNIFORM GRID 6 NODE BD ........................................................ 52
FIGURE 6.13

RADIAL UNIFORM GRID 11 AND 21 NODE BD .......................................... 53

FIGURE 6.14

RADIAL UNIFORM GRID 6 NODE FD ..................................................... 55
FIGURE 6.15

RADIAL UNIFORM GRID 11 AND 21 NODE FD ....................................... 56
FIGURE 7.1

SPHERICAL EQUAL VOLUME GRID 6 NODE CD ................................... 60
FIGURE 7.2

SPHERICAL EQUAL VOLUME GRID 11 AND 21 NODE CD ....................... 61
FIGURE 7.3

SPHERICAL EQUAL VOLUME GRID 6 NODE BD ................................... 62
FIGURE 7.4

SPHERICAL EQUAL VOLUME GRID 11 AND 21 NODE BD ...................... 63
FIGURE 7.5

SPHERICAL EQUAL VOLUME GRID 6 NODE FD ................................... 65
FIGURE 7.6

SPHERICAL EQUAL VOLUME GRID 11 AND 21 N ODE FD ...................... 66
FIGURE 7.7

RADIAL EQUAL VOLUME GRID 6 NODE CD ........................................ 67
FIGURE 7.8

RADIAL EQUAL VOLUME GRID 11 AND 21 NODE CD ........................... 68
FIGURE 7.9

RADIAL EQUAL VOLUME GRID 6 NODE BD ........................................ 69
FIGURE 7.10

RADIAL EQUAL VOLUME GRID 11 AND 21 NODE BD ........................... 70
FIGURE 7.11

RADIAL EQUAL VOLUME GRID 6 NODE FD ....................................... 72

FIGURE 7.12

RADIAL EQUAL VOLUME GRID 11 AND 21 NODE FD ........................... 73
FIGURE 7.13

SPHERICAL PROBLEM MAXIMUM EIGENVALUES ............................... 74

vi

FIGURE 7.14

RADIAL PROBLEM MAXIMUM EIGENVALUES ................................... 75
FIGURE 7.15

SPHERICAL PROBLEM MINIMUM EIGENVALUES ............................... 76
FIGURE 7.16

RADIAL PROBLEM MINIMUM EIGENVALUES ..................................... 77
FIGURE 8.1

RADIAL VERIFICATION AVERAGE ERROR CD ................................... 80
FIGURE 8.2

RADIAL VERIFICATION AVERAGE ERROR BD ................................... 81
FIGURE 8.3

RADIAL VERIFICATION MAXIMUM ERROR CD .................................. 80
FIGURE 8.4

RADIAL VERIFICATION MAXIMUM ERROR BD .................................. 81

vii

CHAPTER ONE

INTRODUCTION
No other field of mathematics has shown a recent increase in importance to the
engineers comparable to that of numerical methods, nor has any other field developed as
rapidly. The main reason for this evolution is the developments in digital computers.
Indeed, each new generation of computers invites new tasks in numerical analysis; in this
connection even a small improvement in the algorithm may have great impact on time,
storage demand, accuracy and stability. This opens up a wide area of research with a
view toward improving accuracy of the software/techniques used for numerical solutions.
Mathematical modeling of physical problems is an important tool in engineering
analysis because it provides the opportunity to study a problem and obtain an
approximate solution without going into expensive and/or time consuming physical and
manufacturing processes. Most of the time-dependent problems in engineering and other
branches of science are modeled in the form of Partial Differential Equations (PDEs).
One group of these equations is referred to as Parabolic or Difiusion Equations, which

have the general form

8U
at

where c is the capacitance coefficient, k is the conductivity/stiffness coefficient and U is

 

= kV . (V U) (1.1)

the unknown variable, that is, temperature, moisture contents, pressure head, and so on.
Equation (1.1) applies to transient heat conduction in solids, gas diffusion/drying
of granular materials, flow of fluids, and transport of solutes in a porous media. Many

engineering and mathematics books deal with the derivation of PDEs and their solution

of the above problems. Powers (1987), Ozisik (1980), Patankar (1980), and
Churchill (1987) are a few examples.

The analytical solution of a partial differential equation is very difficult to obtain
for complicated field problems. Partial differential equations are often converted into a
system of ODEs by applying numerical procedures like the finite element method (FEM)
or the finite difference method (FDM). This conversion of time and space-dependent
partial differential equations (PDEs) into a time-dependent system of ordinary differential
equations (ODEs) has been discussed in many books dealing with numerical solution of
PDEs; some of them are Segerlind (1984), Smith (1985), and Narasimhan (1978). A

system of ODEs has the general form

[C]{U'}+[K]{U}—{F}={0} (1.1)
where [C] is the capacitance matrix coming from the transient term in the PDEs, [K] is
the stiffness matrix coming from the second partial derivative with respect to space and
{F} is the forcing function. Since the forcing function, {F}, in the partial differential
equations is often zero, {F} is zero until the boundary conditions are incorporated.

Finite element or finite difference methods are used to solve (1.2) in the time
domain. The FEM shows clear advantages over the FDM in the space domain in solving
(1.1). This advantage, however, does not extend to the time domain, Segerlind (1984).
There are numerous finite difference schemes available in the literature for solving (1.2)
in the time domain. Different schemes require a different criterion to ensure numerical
stability and to minimize oscillations. Various authors have discussed the solution
procedures in detail but have always based the size of the time step on their art and

experience. The authors seldom discuss the entity of time steps with respect to accuracy.

There is no clearly defined technique, available to select the time step needed to reach an
accurate solution, particularly in two and three-dimensional problems.

Mohtar (1994) pioneered the development of empirical equations that can be used
to estimate the time step required to solve (1.2) accurately when using one of Euler’s
forward difference method, the central difference method or the backward difference
method. Mohtar deveIOped equations to compute optimal time steps using the lowest
eigenvalue of the system of ODEs as the basic parameter. He compared numerical results
with analytical solutions to establish the empirical equations. Time step prediction
equations were developed for one-dimensional problems and two-dimensional problems
where the grid consisted of square elements. Each prediction equation given by Mohtar

had the general form hAt = CNb where his the lowest eigenvalue for the system of
ODEs, N is number of nodes in the region, C and b are empirically determined
coefficients and At is the time step.

Tan (1995) extended the work done by Mohtar into radial coordinates, as a first
step towards solving axisymetric problems. Tan also developed empirical time step
estimate with a form similar to that of Mohtar (1994). Tan pioneered the technique of
using a numerical solution with a highly refined grid in space and very small time steps to
generate reference values. Tan used the central difference method to generate reference
values because it is second order accurate, Gear (1971). Since most complex field
problems either do not have an analytical solution or the analytical solution is also based
on a series solution with truncated terms, use of a numerical technique to generate a set of
reference values seems appropriate. Tan’s approach to generating the reference set

simplifies the research procedure.

The general objective of this study is to extend of the work done by Mohtar

(1994) to spherical shapes and to re-look at the radial field problems studied by Tan

(1995) in the light of recent redefinition of numerical methodology.

The study of transient heat transfer in spherical as well as radial coordinates is

applicable to numerous engineering problems including:

1.

Development of instant heat and its study in the gun barrels, rocket tubes, and
missile launchers, during and after fire. It can also assist us in determining the
optimum rate of fire for a weapon.

The study of a gun shell, movement of projectile in the air and its terminal
ballistics can be facilitated.

Heat dissipation study in piston, cylinder, crankshaft and other components of
automotive engines exposed to combustion or frictional heat.

Cooling or heating of a large number of natural products and the cooling of
processed products in food containers.

Grain drying is governed by a diffusion equation. Accurate numerical schemes are
critical in the optimal design of grain dryers. Parameters such as the time needed
to dry the grain and the rate of drying are critical in determining the dryer

specification. The same can be said about the drying of other organic products.

The results of this study should make it easier to perform a numerical study of these

subject areas.

CHAPTER TWO

REVIEW OF LITRATURE

The solution of the time dependent field problems using the finite element method
was discussed only briefly in early finite element books. Heubner (1975) discusses the
derivation of the capacitance matrix [C], for transient heat transfer but never discusses the
solution of the resulting system of ODEs. Zienkiewicz (1971) and Segerlind (1976)
discussed the numerical solution of the system of ODEs but did not discuss any of the
problems that can arise during the solution process and they did not compare the different
types of elements.

Recent books cover the time dependent problem in more detail but may mislead
an inexperienced analyst. Allaire (1985), concentrates most of his discussion on Euler’s
single step explicit method with one-dimensional problems. This method is known to be
unstable and is not the most accurate of the single step methods. Allaire does not discuss
any solution in two or three dimensions and makes no comparison between linear and
quadratic elements in the one-dimensional case.

Segerlind (1984) discusses some practical aspects of the numerical methods
related to oscillations and physical realities. He warns the reader to avoid using the
quadratic elements because of physical reality problems but does not describe the exact
significance of the errors.

Until recently, most application oriented books in heat transfer and ground water
flow focused their discussion on numerical solutions using finite difference methods and

did not go into detail on the finite element method. Jaluria and Torrence (1986) discuss

the three node triangular element for solving heat transfer problems but do not make any
comparisons with a two-dimensional finite difference solution. These authors do not
discuss any of the other types of two—dimensional elements. Their discussion could lead
one into thinking that the three node triangular element is the most appropriate for a
numerical scheme. Segerlind (1984) indicated that the four node quadrilateral element is
superior to the three-node triangle. The presentation in J aluria and Torrence (1986) can
be contrasted with Patankar (1980) who limits the discussion of the finite element method
to two pages and does not give any equations for the method. Patankar recommends the
use of the backward difference scheme in time due to its "friendliness" for all values of
time and grid size and advocates a control volume approach for the space dimensions.
Patankar (1991) presented a heat transfer “computer program” called CONDUCT. This
program uses the backward difference scheme to solve a heat transfer problem in time,
but Patankar never discusses selecting a time step when solving transient problem.

Shih (1984) has a chapter on accuracy and error bounds. Most of his discussion
analyzes the error bounds for different orders of the finite element method. Shih does not
discuss any estimate for At when solving time dependent problems. He does, however,
have a chapter on the comparison between finite difference and finite element methods.
He covers smoothness of the basic function, numerical instabilities, higher order accurate
discretization schemes and the incorporation of mixed boundary conditions. Shih does
not give any numerical results and concludes with the statement, "Much work remains in
comparing these two powerful methods in a rigorous and conclusive manner".

Shih (1984), Jaluria and Torrence (1986) and Segerlind (1984) avoid explicit

numerical evaluation of the time step. They discuss stability and numerical oscillation

problems but none of the authors gives a procedure for estimating the time step as it
relates to the accuracy of the computation. The typical scenario is to present a numerical
solution procedure and compare it with an analytical solution of the PDE using one or
more time step values. The authors, however, never said how they determined what
numerical value of the time step to use. Dhat and Touzot (1984) comment that the time
step value that eliminates stability and numerical oscillations may not produce accurate
calculations. They also do not give any suggestions on how to select the time step value.
Gear (1971) and Stoer and Bulirsch (1980) discuss the mathematical approaches
to determine a time step value. They define an error as being the difference between two
solutions with time steps of At and At/2 and use this error to determine an appropriate
step size. This approach however, does not give much information on how to select a
starting value for At and requires two or more solutions before a time step is defined.
Myers (1977) discusses the critical time step, applicable to two—dimensional heat
conduction transient problems. His discussion, however, centers on estimating the
maximum eigenvalue for use in the Euler stability criterion or the Crank-Nickolson

oscillation criterion of ammo S 2 where 7cm, is the maximum eigenvalue. Myers does
not discuss the determination of At as it relates to the accuracy of the integration.
Another approach for selecting At is to limit the maximum change in any nodal

value to a certain percent of its previous value. This approach is used in some commercial
finite element software when solving nonlinear problems. This method suffers from the
need to repeat the calculations if the time step is too large and also does not give any

information on how to select a starting value for At.

Reddy (1984) has a section on time dependent problem, that is consistent with
much of the mathematical literature. Reddy describes the stability in terms of the roots of
the characteristic equations and the eigenvalues of the global system. Roots of that
equation should be bounded by one to avoid numerical oscillations. Reddy gives a time

step estimate for structural dynamics problems. He states that At = Tmin/rt, where Tmin is

the smallest period of natural vibration associated with the approximate problem gives an
accurate solution. According to Reddy, another estimate to At can be obtained from the
condition that the smallest eigenvalue of the characteristic equation be less than one.

Smith (1985) discusses the explicit Euler’s method for solving the non-
dimensional form of (1.1). Smith rearranged the difference equation and defined a term
r=8tl(5x)2. During the discussion of stability, Smith stated that the explicit method is
stable for r with values less than 0.5. The implicit Crank-Nickolson has the advantage of
being stable for all values of r. Smith recommends r = l for an accurate solution for the
Crank-Nickolson method. Smith also discussed convergence and stability for some time
stepping schemes and gave a time step expression that satisfies both criteria. No criterion
for selecting a time step based on accuracy was given. The term r defined by Smith does
not include material properties since the thermal diffusivity coefficient, cp/k, was defined
as one.

Allaire (1985) called Smith’s r term the Courant Number. Allaire’s variable
included the material properties. In addition to illustrating stable and non-stable schemes,
Allaire defined an oscillatory stable scheme as having spatial oscillation that eventually
dies out with the solution converging to the correct steady state values. Allaire showed

the following criteria to be true for the single step methods:

0 < r S 0.25 No oscillation
0.25 < r $0.5 Oscillatory and stable

0.5 < r Unstable (Euler’s method only)
The solutions given by Allaire have no indication of instability for values of r < 0.5.
Allaire discussed a “weighted explicit-implicit scheme”. His scheme reduces to the

explicit method and has stability criterion of r < 0. 5 when his parameter 0 equals zero.

Allaire showed that the Crank-Nickolson method and the fully implicit methods are
accurate for values of r up to 1.335.

Jaluria and Torrance (1986) defined Allaire's (1985) Courant number as the
Froude number, F0. These authors suggested using values of F0 less than 0.5 for the
implicit method although lower values gave better accuracy. They never give any
example of what the lower values should be.

Wood and Lewis (1975) studied seven different finite difference time marching
schemes. They comparedimethods based on an accuracy criterion. The authors related
accuracy to oscillations and stability. They determined the critical non-oscillatory time
step for the Crank-Nickolson (C-N) based on the maximum eigenvalue. They showed
numerically that when increasing the time step beyond a critical time step, oscillations
occurred. Wood and Lewis observed inaccurate values in backward difference scheme for
some time step values. They did not state that accuracy is a separate consideration in the
numerical solution of parabolic equations that needed to be addressed and adjusted
accordingly.

Wood (1990) gives an extensive list of time stepping schemes. His list included

most of the known schemes and some new ones. He studied stability, consistency, and

oscillations where the term "time step" was mentioned at several places. For many of
these schemes, numerical results were tabulated using various time steps and the
corresponding error was presented. The author showed that these methods were
consistent with the analytical solution. Wood also refers to the use of time step
adjustment where the size of the time step changes after every set of calculations but

never give any formula for determining a time step value.

Ortega (1990) defined and discussed three types of errors that are all associated
with the time step. The discretization (global) error, convergence error, and rounding
error. He did not indicate how to define the numerical value for the time step that will

minimize these errors.

Rushton and Tomlinson (1971) used the alternating direction approach as a
numerical scheme. They studied stability and found that for different boundary conditions
the Courant number, C, that generates accurate time steps changes. For a sudden change
of pressure head on the boundary, C should be less than 1.0. For a draw down at a well, C
should be less than 0.05. For a sudden change in discharge at a well, C should be less
than 0.5. The authors suggest that a trial and error procedure is still required for selecting

the optimal At value.

Henrici (1977) had an extensive discussion about the error propagation for the
difference methods in solving the PDE. His theoretical treatment did not include

discussion of the time step size needed for accurate results.

10

Williams (1980) and Fried (1979) both studied the numerical solutions of PDE
and used the time step criteria that satisfied stability requirements. Williams used a term
equivalent to the Courant number and stated that it should be less than 0.5. Fried used the
stability criteria (Amman...)

Haghighi and Segerlind (1988) solved the coupled heat and mass transfer
equations using the finite element method. They used Maadooliat’s (1983) non-oscillation
criteria as well as the physical reality conditions that Segerlind (1984) discussed in his
book.

Nripendra and Kunze (1991) presented a finite element solution for temperature
distribution in storage bin. They used the Crank-Nickolson scheme for the time domain.
They presented comparisons between numerical and exact solutions. There was no
mention of time step in their paper.

Irudayaraj (1991) and Irudayaraj et. a1. (1990) applied the finite element method
to the solution of a coupled heat and mass transfer problem. Both papers used the stability
criteria for selecting the time step. There was no check whether this time step ensured
accurate results. The author’s calculated results did not agree with experimental data in
the literature. The same stability criterion was followed by Liu et. a1. (1984). These
authors used a modified Runga—Kutta method to solve the parabolic system. Their work
did not discuss solution accuracy.

Peraire et. a1. (1988) studied the finite element solution of fluid flow. They used

the Courant stability criteria of (AtSK*hJu+c), where c is local speed of sound, he is the

average element length, u is the velocity, and K is a constant.

11

Alagusundaram, et a1. (1991) applied the finite element method to model the
diffusion of carbon dioxide in grain bins. Their calculated results did not compare well
with the measured values. They listed several reasons for this discrepancy. They did not
mention how they determined the time step. They did not state what time step value they
used and did not state whether the size of the time step might be one reason for the
inaccuracy of their calculation.

Cleland and Earle (1984) studied the freezing time of food material using six
finite difference methods. They ensured accuracy by reducing the time step until the
numerical results converged to a consistent value. They encountered a stability problem
and a physical reality violation that they called "jumping" and said it was related to the
latent heat. Although there is evidence of accuracy in their solutions, there is no
evaluation of a time step expression that could be translated to other problems.

Abdalla and Singh (1985) simulated the thawing of food using the finite element
method. They presented comparisons between analytical and predicted values but they
did not state what time step value they used.

Segerlind and Scott (1988) were among the first to deal with the time step
estimates from the accuracy perspective. They presented a time step estimate for one and
two-dimensional problems that produced accurate results. They did not give any
derivation for their estimate and stated that much of it was based on their experience.
They did not show any evidence that their time step estimates really work. However, they
have stated an important observation that a time step based on the oscillation criterion is
conservative. The time step could exceed this criteria by a factor of two before

oscillations were observed.

12

Ne-Zheng Sun (1989) studied numerical solutions for the coupled ground water
flow and advection-dispersion equation. He applied a variation of the linear finite
element method and compared his solution with analytical ones. No indication was given
as to what time step was used in his analysis.

Scientists reporting new time stepping schemes seem to discuss stability and
oscillations only. Yu and Heinrich (1987), Sega] and Praagman (1986), Fong and Mulkey
(1990), Riga] (1990). Schreyer (1981) used the stability time step requirement (At < C
h2/2), where C is the thermal capacitance when performing a numerical solution for the
heat conduction equation. None of these authors discussed the accuracy of the solutions.

Shu-Tung Chu and Hustrulid (1968), and De Baerdemaeker, et al. (1977) did not
define a time step estimate when they discussed the numerical solution of the diffusion
equation. Scott (1987) uses the following arbitrary accuracy criteria At=(time to steady
state)/ 100. In other words Scott assumes that running the problem for 100 time steps
should be sufficient to ensure an accurate solution. Although this estimate might be a
good starting point for some problems, no justification for its use was given.

Maadoliat (1983) studied stability and physical reality oscillations of the finite
element numerical solution. He concluded with a set of conditions that must be satisfied
in order to avoid both numerical problems and recommends a time step estimate
accordingly. He did not consider the accuracy criteria.

Mohtar (1994) was among the first researchers to define the time step value based
on an experimental accuracy criterion. He investigated the one-dimensional problem and
two-dimensional problems consisting of square elements. The general procedure

developed by Mohtar was to:

1. Define a measure of the error,

2. Convert the PDE to a system of ODE using the finite element method in space.

3. Solve a problem using several different values of the time step and several
subdivisions of the problem in space,

4. Plot the error value against the number of nodes and select the time step value, At,
that produced a specified error,

5. Empirically fit an equation to the time step data using lowest eigenvalue as the
basic parameter, and

6. Checked the equations by solving a different set of problems.
In one-dimensional problems, Mohtar (1994), defined the accuracy ratio as

= 2;. 2:.INODE. -APDE..I
23:1 ELI/AODE.) —APDEU|

 

e (2.1)

where NODE is the numerical solution for the system of ODEs, APDE is the analytical
solution for the PDEs, AODE is the analytical solution for the system of ODEs, and n and
m are the number of sampling points in the space and time domain respectively. The
dynamic time step equation developed by Mohtar for the three single step methods
applied to one-dimensional problems are:

Forward difference in time

 

—1.6
At = 0.27 N (2.2)
Central difference in time

—l.18
At=1.13N (2.3)

 

Backward difference in time

-3.91
At = 30.6 N

 

(2.4)

In each equation, At is the time step value, N is the number of nodes in space and
M is the lowest eigenvalue for the system. The time step estimates were validated using

four different problems with analytical solutions: A sine wave variation and a linear
variation in the initial conditions with boundary temperatures known and two problems
with uniform initial conditions and derivative boundary conditions. The problems were
solved using fractions or multiples of the calculated time steps. Time step values of one-

half, two and three-times At were used along with the error ratio

.”' NODE,
e = :1 (2.5)
ZHAPDE,

 

The accuracy ratio for At/2 and At were equivalent. The results for multiples of two and
greater were less accurate than the results for At.

In two-dimensional problems, Mohtar (1994), defined the accuracy ratio as

.. 2m1 lNODEU - APDE,

-- ._ NODE.

e: j-I 1—1 I] (26)
mn

 

 

where NODE is numerical solution for the system of ODEs, APDE is the analytical
solution for PDE, m is the number of sampling points in the space domain, and n is the
number of sampling points in the time domain. The sampling points in the time domain
were at 9.5, 19, 28.6, 38.1, 47.6, 57.1, 66.7, 76.2, 85.7 and 95.2 percent of the time to

steady state defined by tss=(4/lowest eigenvalue).

15

The accuracy ratio, (2.1), used with one-dimensional problems became too
difficult to evaluate for larger two-dimensional problems. Mohtar (1994) restricted the
two-dimensional grid to square elements to allow a comparison of the finite element and
finite difference formulations in space. Using the error estimate (2.5) and a five-percent
error in the calculated values when compared to the analytical solution of the PDE,
Mohtar developed the empirical time step estimates for two-dimensional square grid
given below. Equation (2.7) through (2.9) are for the finite difference formulation in
space while next three are for the finite element formulation in space.

Finite difference method in space

Forward difference in time:

 

 

 

-1.01
At = 1.19 N (2.7)
Central difference in time:
—0.55
At = 1.6 N (2.8)
211
Backward difference in time:
-0.1
At = 0.05—N— (2.9)
11
Finite element method in space
Forward difference in time:
——l.04
At = 1.8N (2.10)
A

16

Central difference in time:

 

—0.55
At=1.6N (2.11)
Backward difference in time:
—0.1
At = ODS—NA.— (2.12)

The above equations are valid when the number of nodes used with the finite
difference formulation in space is equal to or greater than nine. The equations for the
finite element method in space are valid when the number of nodes is equal to or greater
than twenty-five.

Tan (1995) carried forward the work of Mohtar (1994) and developed empirical
equations for calculating the time step required to numerically solve the system of ODEs
related to time dependent radial field problem. The specific objectives of his study were
to develop an empirical time step estimate for the three single step integration methods
that satisfy an accuracy criterion and validate the time step equations by solving different
set of problems.

Tan (1995) used the central difference solution scheme with a very small time
step size and a highly refined grid in space to generate a set of reference values instead of
using an analytical solution. The analytical solution of radial problems involves
numerical evaluation of the series that define Bessel functions. Tan thought it was more
appropriate to simply use a numerical solution to generate the reference values.

Tan’s study used the linear radial element and the three single step integration

methods in time. He used FEM in space and the lumped formulation for the capacitance

17

matrix. The time step equations were presented using the same format as used by Mohtar

(1994). The equations are as follows:

Forward difference method in time:

-2.13
At: 3.91 N

 

Central difference method in time:

N—l.81
A

 

At=.12

Backward difference method in time:

—1.66
At = 5.28 N

 

18

(2.13)

(2.14)

(2.15)

I

CHAPTER THREE
OBJECTIVES

After studying a large amount of literature on the solution of time dependent field
problems, the need for a priori time step estimate seems obvious. An a priori time step
estimate would eliminate the present trial-and-error procedure. A very good time step
estimate would allow the user to generate an accurate solution while satisfying numerical
stability and oscillations criteria. The specific objective of this study was to develop an a
priori time step estimate for three single-step methods used to solve the system of
ordinary differential equations associated with the radial and spherical field problems.
This study uses the finite element method in the space domain and lumped formulation in
time.

The general hypothesis was that the a priori time step estimate has the general
form
(ADA = C
for the unconditionally stable methods and
(ADAM = C
for conditionally stable methods. These equations presented are similar to those
developed by Mohtar (1994) and Tan (1995). The number of nodes in the space
dimension, N, has been deleted as a basic parameter. It is well known that increasing the
number of nodes in space increases the solution accuracy. It is also well known but often
forgotten that the physical parameters that occur in PDEs are generally accurate to two
significant digits and occasionally to three. (of. tables in Perry et a1. 1984). A large

number of nodes in the space dimension significantly increases the solution time and

19

generates accuracy beyond that justified by the number of significant digits in the basic

parameters.
Some of the limits on this study included:

1. Radial and spherical field problems were investigated.

2. Finite element solution was adopted in space domain by using linear one-
dimensional elements.

3. In the time domain single-step solution procedure was adopted. The numerical
schemes used were: The forward difference, central difference and backward
difference methods.

4. The capacitance matrix was formulated by using lumped (diagonal) formulation.

20

CHAPTER FOUR
THEORETICAL CONSIDERATIONS

Galerkin’s finite element formulation was used to obtain the element matrices for
both the space and time domains. The finite element technique is preferred over the finite
difference method because both of the resulting matrices, [C] and [K], are positive
definite, symmetric, and their eigenvalues are real and positive. The stiffness matrix [K]
is singular before boundary conditions are imposed. The finite difference method in space
produces an unsymmetrical stiffness matrix. The global matrices, [C] and [K], are built
from element contributions using the direct stiffness method, Segerlind (1984). The
coefficients in the element matrices depend on the type of interpolation function used to
solve the problem. There are two types of formulation for the capacitance matrix [C]:
consistent and lumped. There are some disadvantages associated with the consistent
formulation (Visser, 1965, Wilson and Nickel], 1966, Brocci, 1969, Zienkiewicz, 1977
and Segerlind, 1984, therefore, a lumped formulation was used for the capacitance
matrix.The primary objective in this chapter is to briefly discuss the derivation of
element matrices for spherical and radial field problems.
4.1 The Spherical Field Problems

Misra and Young (1978) have derived the element matrices for transient heat
transfer in a sphere and have presented adequate details. To maintain continuity for the
reader, especially for the reader who does not have direct access to the reference, a

summary of the derivation is given here.

21

4.1.1 Governing equation
The differential equation for heat conduction in spherical coordinates, with heat

generation within the solid, is given by Carslaw and Jaeger (1959) as

 

 

2 2
Daquzaau+ Drp a[. at!)+ D, 3U 8U

=D— 4.1
rdrz r 3r rzsintpago rzsin2g03¢2+Q 'dt ( )

where D,, D}, D4,, and D. are the physical parameters and Q = Q(r,t) is the rate of heat

generation. The initial and surface conditions and physical properties of a sphere are such
that the isothermal surfaces are concentric spheres, therefore, the temperature is only a
function of the radius, r, and time, t. The problem can now be studied by rotating a one-
dimensional pin in three dimensions and can be integrated over the entire volume. The
radial distance can be divided into a finite number of elements. The three dimensional

equation, (4.1), reduces to

2
D,a(21+20’aU=D,§—t—j— forOSrSRandt>0 (4.2)
Er r Br 3:

 

where Dr is the radial thermal conductivity, U is the temperature and D, is the thermal
capacitance term which is the product of the density and the specific heat. The variables r
and t are the space and time variables respectively. According to Carslaw and Jaeger
(1959) the boundary conditions for (4.2) could be

U=U(r)atr=R;t>O

or a prescribed convection term at the surface

iDr—aaE+h(U —Um):0 at r=R, t>O
r

where h is the surface heat coefficient and U... is the temperature of the surrounding fluid.

The initial condition given by Calslaw and J ager is U = U (r) in 0 S r S R; t = O

22

4.1.2 Variational statement

The finite elements technique requires developing the element equations from
governing differential equations either by obtaining the variational or functional
statement of the physical problem or by directly transforming the governing equation
using Galerkin’s method when the functional statement is not readily available. Misra
and Young (1978) derived the element contributions by minimizing the functional
statement of a physical problem.

The method for obtaining a variational statement from the governing equation is
to rewrite the governing equation in the form of Euler-Lagrange equation, which for
several independent variables has been given by Schecter (1967). The equation Schecter
gives for spherical coordinates is

er a an 2 D aF , av
—————— —— on—-2 U 4.4
an Br Br) +rzsin2¢(d¢)+ ‘ at Q ( )

where F is the function to be determined. In a transient heat transfer problem, the
function F can be split into

F = F, + F3 (45)
where F1 is the function for internal heat conduction and F3 is the function for boundary
conditions.

The function F; in spherical coordinates is,

1 80 2D at] 2 D aU 2 av
F=—D — J— ——J’——— on——2 U 4.6
I 2[ '[ar] r2 (23¢) +rzsin2(p[d¢] + ' a: Q J ( )

which reduces to
F--l-D (>12 Ua—U- (47)
I 2 r pC at '

23

when U¢ f(¢,(p). The variable t} is BTU in (4.7).

r
The function F3 is found from

86:0

8U

 

For the boundary conditions given by (4.3), F3 is
1 2
F B = —2—h(U —U,,)

Misra and Young (1978), define the element stiffness matrix [km] as

7‘
[km]: [[0 am] BIN] )dV (4.8)

 

' Br Br

The element capacitance matrix [cm], is defined the same for all field problems and is

given in Segerlind (1984) as
[cm]: LD,[N]T[N] dV (4.9)

The incremental volume dV in (4.8) and (4.9) is dV=41trdr.

4.1.3 Element matrices

Figure (4.1) illustrates the node locations and elements for the axisymetric heat
transfer problems in a sphere. The interval between adjacent nodes is called an element
and a typical element ‘e’ is the interval between nodal points i and j. The temperature

within an element is assumed to vary linearly and is given by (Myers, 1971)
R. — r — .
U“) = cl") +c§"r = -——’ )Uf‘) + __r R’ U‘." (4.10)
L / L ’

24

The constants c1 and e; have a superscript ‘e’ because these are different for each element
and the superscript ‘e’ with U“) indicates the nodal value in an element. After the

necessary integration, the element stiffness matrix is given by

(e) 3_ 3 __
[K“’]=4D' (R’ R2")? 1] (4.11)
3(Rj-R,) -1 1

 

where R, is the radial distance to node i and RJ- the radial distance to node j.

Using the shape function matrix

Rj—r r—R'.
[N]-[ L L ]

 

in (4.9) and performing the matrix integration, the element capacitance matrix is given by

e (c)
60(Rj —R,.) c2, (‘22

where

cll = 2R3 —20R}R,.3 +301?ij 421:2,5

c,, = 3R3 — 5R3R, + 5R,R;‘ - 3R,5
c2, = c12

c22 =12Rf. —30R}R,. + 20RjR,2 - 2R,5

This capacitance matrix is for the consistent formulation. The lumped formulation is
obtained from the consistent formulation by placing the sum of each row on the diagonal

and placing zeros in the off diagonal positions. The lumped capacitance matrix is

e (e)
[C“’]=—fflo—c— d” O (4.13)
60(Rj—Ri) O d22

where d” = c] 1+ c1; and d22= c.2+ on. The final equations for d1 (and (122 are

25

d“ = 5R3 -5R;‘R, + zokfkf’ + 35R j R,“ —15R,.5
an =15Rj —35RjR,. + zorzj‘R,2 +51?ij —5R,:‘

The above mentioned element stiffness matrix, (4.11) and capacitance matrix,
(4.13) are used to build the global matrices [C] and [K] using the direct stiffness
procedure, Segerlind (1984).

4.2 Radial problem

The field equation in cylindrical coordinates (r,0, z) is

BZU D at] 0,32U azu aU
— '— —— D—— =D— 4.14
'ar2+r3r+r2862+ ZEizzdl-Q ( )

D i
at

 

where D,, De, Dz and D. are physical parameters, Q is the source term and U is the

unknown. If U is independent of 9 then (4.14) reduces to

BZU D Em BZU aU
'— D—— :1)
' 3r2 r 8r + 2 dz2 +Q

D ,—
a:

 

 

(4.15)

If the body is long in z-direction as compared to radius, the end effects are negligible, and

(4.15) reduces to

BZU D 8U 8U
D __+ r __+ _— D __ 4.16
' 3r2 r 8r Q ' a: ( )

 

This equation governs radial heat flow. Assuming that Dr is constant, (4.16) can be

written in the compact form

1 3 3(1) 3U
rl: 'ar rar)]+Q D'a: ( )

26

The boundary conditions associated with (4.17) are either U is constant or the convection

boundary condition
Baa—U: M—Ub +S (4.18)
r

4.2.1 Galerkin’s Finite Element Formulation
The weighted residual integral for left side of (4.17) is the volume integral given

by Segerlind (1984)
{R“’}=-L[NI g_a_ ra—¢ +Q]dV (4.19)
r Br Br /

The solution of field problems for cylindrical coordinates is discussed in several books.
The integral form of the element matrices for radial field problems can be obtained from
the cylindrical formulation by deleting all terms associated with the z coordinate. The
weighted residual integral associated with Galerkin’s finite element formulation for an

axisymmetric element is

{R“’}=[[[D,am 611v] +1) arm am] )dv)U“’}—[ QWW

 

Br 3r 1 Dz dz

—[[N] [0,? —cos sa+D a—Usin6)fl‘
Dz

The above is equation (13.21), Segerlind (1984). Deleting the terms associated with the z

(4.20)

coordinate direction and noting that the outside normal is always perpendicular to the

boundary, cos 0 = 1, the weighted residual integral for the radial element becomes

(e) _ BINIT 8[N] m} 3U
{R }—[L[D, Br Br )dV)U} —jQ[N]dV- L[N]T,[D 37]” (4.21)

The first integral in (4.21) multiplies the column vector of nodal values and defines the

 

element stiffness matrix [km]. The integral containing Q becomes {fm}, while the surface

27

integral is the inter—element requirement for interior element boundaries and the
derivative conditions for the element with a node on an internal or external boundary.
4.2.2 Element Stiffness Matrix

The element stiffness and capacitance matrices for the radial field problems are
not readily available in the literature. Most authors have chosen to write about the
axisymmetric problem, which is solved using two-dimensional elements. The one—

dimensional radial element has the equation as the spherical element

R. — _ .
U‘“ = Cf“ + Cf’r = {—114 r )1,“ + [—r LR' }/;" (4.22)

where R, is the radial distance to node i, Rj the radial distance to node j and L is the

element length. The row vector [N] is

 

 

 

R]. —r r—R.
' 4.2
[N]: [ L L :l ( 3)
and
3[N]=[‘_1 i] (4.24)
3r L L
while
T l
BIN] -2 425
8r 1 (. )
L

The element stiffness matrix [km] is defined by

[___k(e)] “Dre 31”] alNIJdV (4.26)

3,.

 

28

Substitution of (4.24) and (4.25), using dV=27trdr and integrating from Ri to RJ- produces

(c)_27r;Dr 1 —1
[k l- L [_1 1] (4.27)

after using the relationships

Rj+R,.

 

L=Rj—R,. and 7:

and noting that

2 2
R]. —R,.
2

 

=LF
4.2.3 Element Capacitance Matrix

The element capacitance matrix [cm], is defined the same for all field problems
[c“’]=[vD,[N]T[N] dV (4.28)
The matrix of shape functions, [N], for the radial problems is defined by (4.14). .

Substitution of (4.28), using dV=21trdr and again integrating from Ri to Rj produces the

consistent capacitance matrix. Add all of the coefficients in a single row and placing the

value on the main diagonal gives the lumped version of capacitance matrix which is

R. + 2’ 0
[C“’]= —2mD'L ‘ r — (4.29)
6 0 R J. + 2r
4.3 Closure

The element stiffness and element capacitance for the spherical and radial field
problems were developed in this chapter. There are several other details related to the

computer implementation that have not been discussed/included. The direct stiffness

29

procedure for constructing the global matrices and the incorporation of known and/or
derivative boundary conditions. These items are in most finite element books and

duplicating the information seemed unnecessary.

30

CHAPTER FIVE
METHODOLOGY

The research in this study is directed at developing an a priori time step estimate

that meets a required accuracy criteria. Stability and oscillations were given secondary

consideration. It is worthwhile to mention that no numerical algorithm can produce

valuable results unless it is stable and free of oscillations. This chapter discusses with the

methodology adopted to develop an a priori time step estimate.

5.1

Methodology

The methodology in this study is a refinement of the techniques used by Mohtar

(1994), Tan (1995) and Kwon (1998). The basic steps in this procedure are:

1.

Define the Physical Problem

This study was limited to the radial and spherical parabolic diffusion equations.
Each problem was converted to a system of ordinary differential equations by
using the finite element method in space and lumped formulation in time. The
investigation was limited to a cooling problem, which goes to equilibrium.

Define an Error Norm

An average error norm, L1 was defined

 

Ural — U
«i

— , refu-
e_22| (Um-U... (5.1)

 

 

and used in this study. The numerator is the difference between a set of calculated
values in space and time and a set of reference values for the same set of points in
space and time. The denominator contains the total number of sampling points, N,

and the largest difference between the initial condition and a final value. The

31

reference values were calculated using a 21-node grid in space and a time step
smaller than the time step values used to investigate other solutions.

The maximum error in each solution was also determined and retained as useful
data in developing the time step estimate. Mohtar (1994) investigated the use of

other error norms including L2 and L... norms. These norms provided less accurate

estimation of the desired time step.

Select a Defining Problem

The time step estimate is established using a physical problem that is considered
very difficult to solve numerically; a body with a constant temperature (or some
other variable) and the boundary temperature is changed instantaneously to
another value. This problem was selected because the analytical solution of this
problem has all of the frequency components and it has shortest time to steady
state. The initial and boundary conditions for both the radial and spherical
problems were

U(r,0) = 1, 0 S r S 1 (5.2)
with the boundary condition

U(1,t) =0, t>0 (5.3)
Both problems were solved using a radius of one and the physical parameters Dr

and D were also assigned a value of one.

32

Define the Sampling Points in Space

The sampling points in space have to be defined such that the grid being studied
has nodes at the same location as the reference grid. The sequence that satisfies
these requirements consists of 6, 11 or 21 nodes. The reference grid also has 21
nodes but smaller time steps. This set of grids is illustrated in Figure 5.1. The
calculated value at each internal node in the 6, 11 and 21 node grids was used in
the calculation of the error norm.

Define the Sampling Points in Time

The sampling points in time were defined at or very near 1/ 14, 2/14, 3/14,4/ 14,
5/14, 6/14, 7/14, 8/14, 9/14, 10/14, 12/14 and 14/14 of the 60 percent of the time

to steady state, tss, which was calculated using

I = —— (5.3)

where A. is the lowest eigenvalue of the system of ordinary differential equations

associated with the transient solution. The definition of tSS comes from, e4" , the
first term in the analytical solution of a parabolic diffusion equation. This
particular term lasts the longest. The value of 4 is used in (5.3) because when Mt
= 4, e'4 = 0.018 and over 98% of the transient has been completed.

Select the Solution Procedure in Time

There are several methods available for solving a system of ordinary differential
equations in time. Three single step methods were investigated in this study; the

forward difference, central difference and backward difference methods.

33

SAMPLING POINTS IN SPACE

 

 

 

 

 

 

 

1 2 3 4 5 6
C G C Q C O 6 Node Grid
11 ll Node Grrd
21 Node Grid
21
Figure 5.1

34

7. Define the Accuracy Level
A priori time step estimate must be defined relative to a level of accuracy. The
equations in this research were defined for an average error of one-percent and a
maximum error of approximately five-percent.
5.2 Computer Software
A finite element program, Segerlind (1987), was modified for use in the
numerical experimentation. Some of the results obtained from the Quickbasic®l program
were also verified using the commercially available mathematical program, MATLAB®2.
The Quickbasic programs solved the PDEs in space generating the ODEs and then
calculated the lowest and highest eigenvalues for the system of equations. The program is
capable of evaluating four different capacitance matrices including the lumped and
distributed (consistent) matrices. Only the lumped formulation was used during this
research. The computer programs were enhanced to minimize the amount of time

required for file management and graphing of the calculated values.

 

I Registered Trademark of Microsoft Corporation.

2 Registered Trademark of The Math Works, Inc.

35

CHAPTER SIX
RESULTS: UNIFORM GRID

A large number of experiments was conducted to determine the desired At and its
relationship with the lowest or largest eigenvalue of the system of ODEs associated with
the spherical and radial field problems. Each subdivision of this chapter covers the
experiments related to one of the field problems and the solution procedure in time; the
central difference, backward difference and forward difference single step numerical
schemes. Each time scheme was investigated using grids of 6, 11 and 21 nodes. Sectors
6.4 and 6.5 deal with deriving equations for the a priori time step.

6.1 Reference Values

Reference values were generated by solving a twenty one-node grid using the
central difference single-step numerical scheme. The time step (At) was kept equal to
0.6tss/ 112 for the unconditionally stable methods and At = 0.6tsJ8960 was used for
conditionally stable method. Numerical solutions with 31 and 41 nodes were abandoned
because there was minimal improvement in accuracy with a considerable increase in the
computational effort.

6.2 Spherical Shapes

For the unconditionally stable schemes, central difference and backward
difference, the 6, 11 and 21 node grids were solved in time using a At that varied from
0.064 (0.6tssll4) to 0.008 (0.6tsJ112). To maintain compatibility of the sampling points in

time domain, sampling points for At =0.6tss/ 14 were 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12 and 14;

36

similarly for At=0.6ts,/56 the points were 4, 8, 12, 16, 20, 24, 32, 36, 40, 48, 56. The
multiples of 14 used in the experiments were 28, 56, 70, 84, 98 and 112.

For the conditionally stable forward difference scheme 6, 11 and 21 node grids
were solved in time using a At that varied from (0.6tssl896) to (0.6tssl8960). All At were
kept in the multiples of 896 (896 is also a multiple of 14) to maintain compatibility of the
sampling points. For example, sampling points with At =0.6t,S/896 were 64, 128, 192,
256, 320, 384, 448, 512, 576, 640, 766 and 896, and so on. The multiples of 896 used in
the experiments were 896, 1792, 3584, 5376, 7168 and 8960. A larger number of time
steps was required for the forward difference scheme because of the stability criterion.

Every solution was compared with the reference in space as well as time domain
by calculating a L1 norm, the average error and deterrrrining the maximum error.
Summary of all the data used in experiments and the values obtained are attached as
Appendix A.

The results were plotted for every reading obtained during experiments. The time
step was kept on the horizontal axis and the average or maximum error was kept on
vertical axis. Details of data including variation of temperature at each sampling point in
space and time are lengthy and not attached with this document; however, a summary of
results, graphs and analysis for each experiment is presented.

6.2.1 Central Difference Method in Time
6.2.1.1 Six Node Grid

Figure 6.1 graphically represents the average error for a six-node grid. The six-
node grid did not produce accurate results. The average error remained above two-

percent even for smaller time steps. A more refined grid in space is recommended.

37

S semi

 

 

a... s8 3 aamaeafisz
: 3 en E 3 we a:
me"
E
m..."

 

 

In
N

 

 

 

I \

 

 

 

In
"2
N

\

 

 

me.~

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

fin

 

MUZmEm—nEHQ SHZMU QED SEOEZD A<U~mmmmm

10.1.13 afieraxv

38

6.2.1.1 Eleven Node Grid

The average error for an eleven-node grid is presented in Figure 6.2. The eleven-
node grid produced very accurate results with the average error remaining below one
percent for all time steps. The error increased sharply for time steps larger than 0.6tssl28

or At>0.0032.

6.2.1.2 Twenty-one Node Grid

Figure 6.3 represents the average error for the twenty-one node grid. The results
were very accurate, with the average error remaining below 0.2 percent for several time
steps. The error has a sudden increase at time steps larger than 0.6t,,/28 or At>0.0032 and
exceeds the error in the eleven-node grid.
6.2.1.3 Backward Difference Method in Time
6.2.1.4 Six Node Grid

The average error is presented in Figure 6.4. The six node grid did not produce
accurate results. Average error remained above two- percent even for the smallest of time
steps.
6.2.1.5 Eleven Node Grid

Figure 6.5 graphically represents the average error. The eleven-node grid

produced accurate results below At=0.016 where the average error remained below 1
percent. The error increased with the time steps larger than 0.6t,,/28 or At>0.0032.

6.2.2.3 Twenty-one Node Grid
The average error for the twenty-one node grid is presented in Figure 6.6. The
average error was below one percent when more than 56 time steps were used. The error

increased for the time steps larger than 0.6t55/28 or At>0.0032.

39

 

 

3.82 2 101

 

 

3 Sam

 

 

 

 

mmh. 95.3 3 macaw Ho .3252
new an wm ch vw «a Nam
i / \4...

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUZmamEmHQ AEHZWU DEG EMOBZD A<Uammmm

2.
N...
n...
v...
m6
ed
5.:
«6
a...

10.1.13 afimaav

4O

3

ma

3 23E

as as 3 9.3 .e 3.152

cm

an

um

mm

«3

 

 

\

\T

l+

 

\

 

 

 

metal 5 le

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUZmMmEEQ gHZmU QED EMOEZD A<Uammmm

‘9 “I (‘1
c c c
10.1.13 afieraav

“2
a

N4

41

we 65mm

a... s8 a 255 e 34.52
3 ma em 3 an we a:

 

 

 

 

 

 

 

motel eIOI

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUZmMmEnEQ QM<>>MU<m GEO EOEZD A<Uammmm

10.1.13 afieraav

42

3 23E

 

 

 

 

 

 

 

 

 

 

 

 

 

as as 3 3% .e .352
3 ma em 3. 3 ma N:
c
m...
\e
\T\ —
t m g
822 = 11 \
\\ N
\ m N
m
.\
Wm

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUZmEm—mEQ QM<>>MU<m QED EOEZD A<Uammmm

10113 afimaav

43

 

 

manual @ 1.1

 

 

he 21am

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

man. 6&3 3 38m he hon—~52
can Nah“ 3mm ehmm we: 83

c
m...

H
m4

N
m.N

11.111.11.111.
‘ m
\

:\ m n

.4

 

 

MUZmEmEEQ Qm<>>m0m AHMU EMOEZD A<Uammmm

10113 afimaav

6.2.3 Forward Difference Method in Time
The forward difference is a conditionally stable method in which stability and
oscillation limits depend upon the maximum eigenvalue, 714m, of the system of equations,

Segerlind (1984). The reference values were obtained by using a very small time step,

At=0.6tss/8960.

6.2.3.1 Six Node Grid

Figure 6.7 graphically represents the six node grid. Once again, it did not prove to
be accurate and average error remained above one-percent even for small time steps.
6.2.3.2 Eleven Node Grid

The average errors for the eleven node grid are plotted in Figure 6.8. The results
were accurate and the average error remained below one percent for smaller time steps
but suddenly increased for the time steps larger than 0.6tss/35 84.
6.2.3.3 » Twenty-one Node Grid

The average error for the twenty-one node grid is presented in Figure 6.9. The
results were very accurate. Average error remained below 1 percent but suddenly
increased for time steps larger than 0.6tssl3584.
6.3 Radial Shapes

The procedure and schemes used for radial shapes are similar to those used for the
spherical shapes. Results of eleven and twenty-one node grids have reasonable similarity;

therefore, they have been combined for the purpose of clarity and brevity.

45

 

m§Z¢+

 

 

 

he 2:3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

man. $3 3 Roam He .8852
wan «a: 3mm wumm we: 8%

c
m6

_
m4

N
Wm

little
\‘ m
\

e\ Wm

.4

 

 

 

 

 

 

 

 

mUZm—MmEEQ Qm<>>m0m Q20 320.323 A<Uammmm

10113 afimaav

46

we 2:3

a... s8 3 Ram 223.52
23 «a: 3mm 3% we: 8%

 

man...

8N6

 

3N6

 

wané

 

 

 

$662 3 lol

#\
8
N
c

 

 

 

put
"1
G

 

3N...

 

N35

 

 

 

 

 

 

 

 

 

 

 

 

 

 

m—Né

 

mUZm—mmmnflm Gym/3425:?“ QED EMOEZD A<Uammmm

10113 afieraav

47

we 8:5

 

 

 

 

 

 

 

 

 

 

as as 3 .85 he his:
wax flab“ 3mm 95mm weak. saw
a
\111111111fi
~66
\\

\ Q. a V
\ m
motel “N I m: A. 9%.
\ e
w
\ 2:. 1

\ as

.\
2.6

 

 

 

 

 

 

 

 

 

 

 

 

 

 

WUZm—MEEQ Qm<>9m0m QHMO EMOEZD A<Uammmm

48

6.3.1 Central Difference Method in Time
6.3.1 Six Node Grid

Figure 6.10 exhibits the average error for the six-node grid. The results were
inaccurate; average error remained above two-percent.
6.3.1.1 Eleven and Twenty-one Node Grid

Figure 6.11 displays the plot of average error for the eleven and twenty one-node
grids. The results for the time steps up to 0.6tsJ56 are accurate with the average error
remaining below 0.5 percent. There is, however, a sharp rise in the error for the time
steps larger than 0.6t,,/28.
6.3.2 Backward Difference Method in Time
6.3.2.1 Six node Grid

Figure 6.12 graphically displays the average error for the six node grid. The
results are inaccurate and suggest need for more refined grid in space.
6.3.2.2 Eleven and Twenty-one Node Grid

Figure 6.13 represents the average error for the eleven and twenty-one node grids.
The results were accurate with the average error remaining below 0.5 percent for several
time steps. The error increases suddenly at time steps larger than 0.6t35/56.
6.3.3 Forward Difference Method in Time
6.3.3.1 Six Node Grid

Figure 6.14 displays the average error for the six node grid. The results remain

inaccurate and further strengthen the idea of needing a more refined grid in space.

49

SeaAa

as s8 3 83m .2352
3 «a em 2 3 we a:

 

 

 

_7r \/

 

 

 

$602911

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUZm—MmEEQ 1395.750 QED EOEZD A/EDE

me.—

n.—

mu;

10113 afimaav

50

:ea&E

as .88 8 865 2235.2
3 ma em E S we a:

 

 

 

 

\Tillxfl lfi 1%

r r rillhil-

 

 

 

 

3.32 3 lxl
motel gala:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

muzmmmmmwfl 35750 QED EMOEZD A<HQ<M

In W t") N —-1 G
10113 afimaav

@

h

51

e e 65E

a... as 3 85 .e 3.252
z 3 em 2. a as a:

 

 

 

 

 

mgzell

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUZm—Mmmaa 23M0<m GEO $30523 Eng

9 M N v-1 6
10113 afimaav

In

\D

52

3

a... s8 6. Esme .352

cm

268%E

ch

3

ma

N:

 

 

K“

 

\

 

 

 

 

 

 

 

\

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUZMMmEEQ QM<>>MU<m QHMU EMOBZD A<HQ<M

m.N

m.m

10113 afimaxv

53

Se 8:3

as .88 8 8% e .352
ea 9: 5mm 3% we: 8%

 

 

 

 

 

metal 9 le

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

muzmmmmafl gag/Mai GEO EMOEZD imaé

mam.—

vwé

mew.“

mw.~

www.—

so
°°.
v-t

10113 afimaav

mew;

hm.“

mhwé

54

6.3.3.2 Eleven and Twenty-one Node Grid
Figure 6.15 represents the plot of the average results for the eleven and twenty

one node grids. Average error remained below 0.5 percent for the smaller time steps but

showed a sharp increase for the time steps larger than 0.6t55/3584.

6.3 Analysis
After considering the individual details of each experiment, I summarized the

results and developed the time step equations. The salient observations on the results are

as follows:

1. The most common feature among all the calculated results is the fact that the error
ratio decreases with decrease in the size of the time step, At. In other words, the
accuracy of a numerical solution improves with the decrease in At. There is a At
for some problems, however, where the error starts to increase.

2. In every scheme, there is a time step below which the improvement in accuracy is
not significant and the reduction in the step size becomes counterproductive due
to extra computational effort and possible round off errors.

3. There is definitely a region for At, where the solutions are within the desired
accuracy range and show little shift with the variation of the step size. This region
is between At>0.6ts,/84 and At<0.6ts,/28. This region was used to develop the a
priori time step estimates.

6.4 Derivation of the A Priori Time Step Equations
The a priori time step estimate equations were developed using the numerical

experiments conducted and analyzed in this chapter. Since the equations have been

55

Beeam

ea. .88 s 85 .e .3952
3» a: 32mm 2% 8: Saw

 

 

\

 

 

 

332 3le
3602 3 II

N

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUZmMm—mEQ QM<>>MOL QHMU EMOEZD grandma

1011a afimaav

56

derived from step change problems, the hypothesis is that they can be used for problems
with derivative boundary conditions.

Keeping in view the intricacies of the numerical schemes used, a separate
equation has been developed for each scheme.
6.4.1 Central Difference

A large number of numerical experiments were conducted and analyzed. The
analysis shows that accuracy level of less than one-percent has been obtained for all
analysis up to a time step as larger as 0.6t55/28. It has already been elaborated that tSS has
been obtained from the lowest eigenvalue it]. Based on the relationship and results
obtained during the numerical experiments the empirical equation deduced for the central
difference time scheme is as follows:
AM, = 0.025 (6.1)
6.4.2 Backward Difference

The analysis indicates a lot of similarity between the backward difference and the
central difference methods. However, it was observed that the average error remains
within the limit of less than one—percent in the range of a time step around 0.6tss/56.
Therefore the empirical equation obtained for the backward difference time scheme is
slightly different than the central difference:
At/l1 = 0.05 (6.2)
6.4.3 Forward Difference

Very small time steps were used in the forward difference scheme due to the
problem of stability and oscillations associated with this scheme. Therefore, the results of

forward difference method remain within the limit of desired accuracy. After deliberate

57

analysis it is has been found that the oscillation criteria for the forward difference method

is also the accuracy criteria. The empirical equation is as follows:

Ar=— .
2 (63)

These a priori time step estimation equations should provide a reasonable start
point in numerical solution of the parabolic diffusion equations that approach to
equilibrium. They should save considerable time by eliminating the present trial and error
forecasting of At and eliminating the computationally expensive trials with very small

time steps.

58

CHAPTER SEVEN
RESULTS: EQUAL VOLUME-GRID

This chapter discuses the solution of the spherical and the radial problems using
an equal volume grid instead of a uniform grid. Details about the reference values, the
error norm, and the solution procedures were identical to those discussed in Chapters
Five and Six. Therefore, only the results of this study are presented in this chapter. A
summary of the eigenvalues, time steps, average errors, maximum errors and other useful
data is presented in Appendix B.
7.1 Constant Volume Spherical Grid

7.1.1 Central Difference Method in Time
Figure 7.1 graphically represents the average error for the central difference

method and a six-node grid. The results were similar to the uniform grid. The average
error remained above 1.4 percent for smaller time steps and showed a sharp increase for
the step size larger than 0.6tsJ28. The average error for eleven and twenty-one node grids
is presented in Figure 7.2. The average error remained below one percent for the time
steps smaller than 0.6tss/56, and increased rapidly after that.

7.1.2 Backward Difference Method in Time

The average error for the backward difference method and the six node grid is
presented in Figure 7.3. The average error remains above 1.4 percent even for the
smallest of the time step value. The average error for the eleven and twenty-one node

grids is presented in Figure 7.4. Accurate results occur for At less than 0.6tsJ56 where the

average error remained below 1 percent.

59

E 8:5

EH. $00 9 macaw E 39:52
3 mm on on em S _

 

 

 

882 e 1T

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUmemmEQ
1333.750 QHMO mEDAO> A<DOm A<Ummmmmm

10.113 afimaxv

60

S 8&5

 

 

3H 68% 9 3on Co cone: Z
3 mm on 3. 4w NH _

 

\

 

 

882 a IT

882 : In:

 

 

 

 

\\

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUZmEmEEQ

A<MHZmU QED WEDAO> A<DOm A<Ummmmmm

10113 ofimoxv

61

2 Beam

as see 9 seem a
3 mm en 2 4w 2 _

 

O

 

"1
o

 

monoZ 0 110.1

 

 

 

 

 

 

1

10.113 aBeroxV

 

 

\

 

M

 

V?
m

 

 

 

 

 

 

 

 

 

 

 

 

 

 

v

MUZmMmEEQ QM<>>MU<m QED MEDAO> 1?.de éUEm—Emm

62

3

mm

3 63mm

as see a 88m 4
em on em

N:

 

 

 

 

 

\ 11111

 

882 a IT
882 : IT

 

 

\

 

 

S

V

 

K

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

V?
(\l

Wm

MUZm—mmmaa 23M0<m QHMO WEDAO> A<DOm invammmm

10013 oSeraAV

63

7.1.3 Forward Difference Method in Time

The time steps had to be kept very small due to an increase in Km”, for the system
of differential equations. The stability and oscillation criteria depend upon the Am“.

The average error for the six node grid is presented in Figure 7.5. The average
error remained above one-percent even for all time steps. The average errors for the
eleven and twenty one node grid are displayed in Figure 7.6. The results were accurate
and the average error remained below one percent for smaller time steps but increased
rapidly for time steps larger than 0.6tsJ3584.

7.3 Radial Shapes
7.3.1 Central Difference Method in Time

Figure 7.7 exhibits the average error for central difference method six node grid,
the error remained beyond two-percent. Figure 7.8 displays the plot of average error for
eleven and twenty one-node grids. The results for the time steps up to 0.6tss/56 are
accurate with the average error remaining below 0.5 percent. However, the error
suddenly went to 3 percent for the time steps larger than 0.6tsJ28.

7.3.2 Backward Difference Method in Time

Figure 7.9 graphically displays the average error for six node grid and. Figure
7.10 represents the average error for the eleven and twenty-one node grid. The results
were accurate with the average error remaining below 0.5 percent for several time steps.

The error had a sudden increase at time steps larger than 0.6t55/56.

cam

mm:

3.2. 688 9 £28m me 59:52

vwmm

2 8:5

05mm

we;

comm

 

 

 

 

 

 

movoz 0 lol

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mom.—
mom;
hm;

Km;
thé
mum.“
th;
mum;
chm;

Rm;

10113 08819AV

mUmem—mnzfl QM<>>MOm @220 WEDAO> A<DOm A<UHmmImm

65

e.» 8:5

 

 

 

 

 

 

 

 

 

 

 

 

 

 

we. see 9 seem? 352
gm 8: £2 8% we: 8%
la o
\llllll.
\ 8 o
3
\ . w
882 2 LT \ 2 c em
882 a IT \ N o m
\ m
\ mg
no
T At 11 L. ..
. moo

 

 

 

 

 

 

 

 

 

 

 

 

MUZm—mmmafl QM<>>MOL EMU m2340> 1?.de A<UHmmEmm

66

E 83E

as 38 3 8% .3 8:52
2 «a em 2 em a:

 

 

 

 

 

moceZelnl

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUZm—mmmmfin—
4.493.750 QED MEDAO> HZ<HmZOU A<HQ§

1-1

"2
p1

N

"2
N

10113 afimaav

67

w.» 8&3

as 38 3 332.. e352
2 an em E .3 a a:

 

 

 

 

 

s\\i\ . 11

 

\

 

 

3992 “NI
mucez 3 Inl

 

\

 

 

 

Q

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

WUmemnEHQ
139;va Q20 WEDAO> HZ<HmZOU A<HQ<M

G

"2
6

F4

"2
p1

N

"2
N

I")

"2
m

10113 afimaxv

68

as 8:33

as 38 3 anew Le 33.52
3 an em 2. 3 a:

 

 

 

 

 

 

 

 

 

m§Z©+

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUZm—MmEED
QM<>>MU<m QHMU MEDAO> HZ<HmZOU A<HQ<M

IQ an IQ N at) —1 IQ c
m N —1 a

V

101.13 aSeJaAV

69

OE 813E

as 38 3 83m 2.. .5352

 

 

 

2 ma em 2 em a:
\
\\.Hl/.\\\.
ill

 

 

 

3.52 “N Lal
3.52 3 II

\1.\\

 

 

\

 

 

\\

 

 

 

 

 

as

 

 

 

 

 

 

 

 

 

 

 

 

MUZm—mmnEHQ

QM<>>MU<m Q2230 mEDAO> HZ<HmZOU 132mg

F1

10113 afimaxv

70

7.3.3 Forward Difference Method in Time

Figure 7.11 displays the average error for forward difference method in six node

grid. The results remain inaccurate. Figure 7.12 represents the plot of the average results

for eleven and twenty one node grids. Average error remained below 0.5 percent for the

smaller time steps but shows a sharp increase for the time steps larger than 0.6tsJ3584.

7.4

Analysis

A large number of experiments were carried out on the problems already solved

by using the equal volume grid instead of the uniform grid. Radial length of the elements

was varied in such a manner that the volume of each element was kept equal. Review of

the analysis of these experiments is as follows:

1.

2.

There are no perceptible advantages associated with the equal volume grid.
However, the major disadvantage noticed during the study was a noteworthy
increase in the maximum eigenvalues of the system of equations generated by
using the equal volume grid. This phenomenon was observed in both, the
spherical as well as the radial problems.

Figure 7 .13 presents a graphical comparison of the maximum eigenvalues of the
equal volume grids with uniform grid for spherical problems. Similarly, Figure
7.14 gives the comparison of maximum eigenvalues associated with radial
problems.

The variance in the minimum eigenvalues of the system of equations using equal
volume grid as compared to uniform grid was not significant.

Figure 7.15 displays the comparison of the lowest eigenvalues of the equal
volume grid with the uniform grid for spherical problems. Figure 7.16 gives a

similar comparison for the lower eigenvalues associated with radial problems.

71

 

:s 2ng

as $8 3 Ram... 3252

wan Nah“ vamm 95mm we: Saw
mmé

 

RNA

 

«mud

 

mmaé

 

 

 

mote Z e lo! emu;

 

 

 

 

mmNA
\

\ meA

 

 

5mm;

 

 

 

 

 

 

 

 

 

 

 

 

 

meA

 

mUZmMm—mnzfl Qm<>>m0m Q36 szqO> A<DOm A<HQ<M

10.1.13 afimaAv

72

can

 

Q:

N: 2ng

as as 3 33222.52

vwmm

:5“ a:

E.

:3:

 

\

\Il

11111:

 

\=\\

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

N:.:
v:.:
::.:
::.:
fl:

NH:
3.:
3.:

mm:

mUZmEm—mam Qm<>wm0m mam: m=>540> A<DOm A<Hfl<m

mug] afieJaAv

73

m: 2:5

—N WHGCZ mu 9

:
llllllll.\l\4 s.
\ a:

\ 82

 

 

 

 

 

 

«Ea—3» =35— Iul \ :::N

are Eofiall \

 

 

 

 

 

\ 8mm
\ 8%

:::v

 

 

 

 

 

 

 

 

mmDA<>ZmUHm EDEUOQE qumOME éUEm—mmm

74

is 03$

5 $302 3 e

 

 

 

 

 

9:23» :35 lnl
2.5 5.5:: IT

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mmDA<>ZmDHm 23222 qumOME invammmm

ed

m.N

9N

EN

QN

ad

fin

NM

75

“N Wynn—oz fin w

 

 

 

 

 

gun—ON! 135nm .IDI
mum—U steam-.3 IOI

 

 

 

\

 

\

 

 

 

 

 

 

 

fl

 

 

mmDA<>ZmOHm EDEN/:2 qumOME imam

c9:

§N

8cm

8:?

86m

$8:

76

o: 2:5

 

 

 

 

 

 

“N = e
wm—GO” maxm
m.m
mm.m
v
«Ea—Pr 33M— Iul \ Wm
E6 .583. IT \\u mew

 

 

 

 

hm
llllll.\\ mnm

 

 

 

 

 

 

 

 

Wm

 

mm54<>ZmOHm ESE/mg qumOMQ 4.3%

77

CHAPTER EIGHT
EVALUATION OF THE TIME STEP ESTIMATES

This chapter recommends a procedure to apply the empirical time step equations
for the solution of spherical and radial field problems. This chapter also deals with the
verification of the time step estimates by applying them to the numerical solution of
problems different than what were used during the numerical experimentation.
8.1 A Priori Procedure

The recommended procedure for handling the numerical solution in the light of
the prepriorri time step equations is as follows:
1. Calculate M and A...“ for the system of ordinary differential equations developed

using the finite element method or finite difference method in space.
2. Calculate At for the intended numerical scheme using equations 6.1, 6.2 or 6.3.
3. Round At to a convenient value. Round the value down for step change boundary

conditions. The value can be rounded up for convection boundary conditions.
4. Solve the problem using the numerical scheme selected, printing the calculated

values as desired.
8.2 Evaluation

The time step estimates presented in Chapter Six were developed using the step
change problem on a solid sphere and solid cylinder of radius one and assuming
Dr=Dt=1. The ability of these equations to predict the time step for real problems was
evaluated using different materials and other boundary conditions.

The comparison problems were chosen with real values of the material

properties and realistic dimensions with the convection boundary conditions. The results

78

proved that the empirical a priori time step estimate equations produce reasonable results
and maintain the desired level of accuracy. A number of experiments were conducted to
verify the developed a priori equations. However, for the purpose of clarity results of a
radial problem with convection boundary conditions and actual physical properties are
displayed in the form of graphs and analyzed.

The Figure 8.1 presents the average error for eleven and twenty-one node central
difference method. It is clearly evident that the desired accuracy level of less than one
percent is available for the value of At lesser than 0.6tss/28. Similarly Figure 8.2 displays
the average error for eleven and twenty-one node grid using backward difference scheme
and verifies our a priori time step estimate equation.

Figures 8.3 and 8.4 present the plot of maximum error for central difference and
forward difference methods. It is indicated by the graphs that the maximum error remains

with in the specified range for both the schemes.

79

S mama

as $8 3 Eem *
3 ma em 2 an ”a a:

 

 

 

 

 

 

 

36:2 gulxl
36:2 3 Inl

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mUmemmmHQ
AEHZm—U qumOME ZOHH<UEE<> A<HQ<M

". ‘9 'Q V: "2 "i '1 °
: : : : : : :
10.1.13 aSmaAv

°°.
:

80

3

:N

3 05mm

as as. 3 £3 a

:m . :h

v:

a:

N:

 

 

 

 

k“

i“

 

 

 

 

mogz .nm '01.

 

 

 

 

 

R

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

QM<>>MU<m qumOME ZOHH<UEE<> 464nm;

mUmem—mafl

md

m.m

10113 afimaxv

81

mm 05mm

as .58 3 2:5 a

3 :m :b v: a: N:

 

 

 

 

 

 

 

 

motel “N lxl
motel S IOI

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

:N

\l . .

\ m
\\ .
\ mm
:N

mm

mUmemmEQ

AEHZm—U Emu—mafia ZOHH<UEE<> Jain—«am

HOHHH XVW

82

3

em

.2 25mm

a... $8 3 flew a.

:h

v:

a:

N:

 

 

lrll

 

R‘i‘k

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MUZm—MEEQ
QM<>>MU<m EmAmOME ZOHH<UEHM<> A<HD<M

:u

m—

:N

mm

:m

HOHHEI XVW

83

CHAPTER NINE
DISCUSSION AND CONCLUSION

The empirical equations for calculating the prepriorri time step estimate to solve
the system of ordinary differential equations for spherical and radial field problems have
been successfully accomplished and verified. The prepriorri time step estimates, meet all
accuracy criteria or the stability requirements (equations 6.1, 6.2 and 6.3).

The forward difference method is conditionally stable; therefore, both 7&1 and km”
must be evaluated, 2.1 is utilized to acquire the tSS and 2mm is used for estimating the
prepriorri time step. Time step value calculated for forward difference method should be
rounded down because of the stability and oscillation criteria. The time step values for
the central and backward difference schemes can be rounded up or down, to a suitable
value for easy division by integers and displaying the results.

The prepriorri time step estimate equations define a time step that will give
overall accurate results and not at any particular time. It has been observed that the error
distribution is not linear through out the time domain. The error is highest at small time
values and dies out as time increases. The distribution also changes with the integration
scheme. Some schemes are more accurate than others at the same point in the time
domain. The fact that the error is highest at small time values might explain why in the
central difference scheme increasing the time step did not reduce the accuracy. In the
later scheme, the inaccuracies present at small time values will not be included for large
values of the time step. On the other hand smaller time steps capture the numerical error

present at small time values.

84

10.1

Future Study

It is recommended that two and three-dimensional problems be studied in details
and validity of the existing empirical equations should be observed on them and
prospect of evolving new empirical equations may also be carefully thought
about.

Application of the empirical prepriorri time step equations to composite materials

and materials with multiple layers of insulation should also be studied.

85

:.~:_
0.8—

o.~:_
:.~:_
0.8—
0.8—
:.~:~
8.8—
8.8—

o.~:_
:.~:_
c.~:_
0.8—
8.8—
:.~:_

I83

omwn.n
cmwmd

cmwwd
:wwn.m
cmwmfi
:mwn.n
omwnw
owwmd
:mwm.m

cmwm.m
cmwmw
:mwm.m
:mwm.m
omwwd
omwmw

£2153 nonhuman—z

8w.m~_
9862

v.08:

~58:
nn~.w:
8m.w:
08.3—
:w:.:~_
~28:

v2.3
5.8.2
Sm:
cad
:8
ww.m

~80.—
vM~.~

con.—
n _ :4
new.—
mm.—
:2..—
8.—

em.—

:8.-
how.—
vow.—
8:.—
www.—
www.—

uehmu><

:9.— 88.:
83.38.:

8888.:
882. 8.:
88:8:
8:02.88
8888::
m~3 .88.:
82.38.:

8:38.:
$2.88.:
88m 88.:
n:-m:-m.o
m:-m-oo.c
m:-Mwmom.m

H. «:09

~80:
~80:

~80:
~8n.:
~8m.:
~8m.:
~8m.:
~8n.:
~80:

~80:
~80:
~8m.:
~8m.:
~80:
~80:

8
~:

3
m~
on
on
v:
:8
~_ _

0::

~82
8mm
85mm
8. N.
80%

mmh.

3.5 535:5

3:5

 

< EDmem<

86

20:;
5:;

5:;
30:;
5:;
20:3.
3:;
3:;
5:;

5:;
20:;
#0:;
5:3.
20:;

3:;

:.~:_
0.8.
o.~:_
o.~:_
05.

I83

o~m~um
o~mfin

a~mfin
m~mfim
o~m~un
o~mfim
o~mfim
o~mfim
o~m>.m

32.6
a~mfim
o~mfin
a~mfin
o~mfim
:~m\..m

cmwmw
cmwww
emwmfi
cmwmd
omwnfi

:92.
“he:

50::—
20.2
h _ w.—
3:
$0:
_wo.~
now.—

3.0:.
coca
v _ :. m
wo.~
hv~

v.—

wood:
5::
30mm—
no.9”—
mn _ .n~_

524:3 beam—:2

36.:
30:

:3”.—
m-.:
m:_.:
:3:
3 _:
50—:
~_ _.:

3::
$5.:
3::
5::
~:_.:
:2:

mend
~bm.m
: _ 0~
mov~
vv~.~

hob—H34».

mg.— 88.:
83.35::

:::::mm:.:
::::2. 5.:
:::now:::
8:: _ 5::
50:05::
o~3 38.:
53.35::

$625.:
$2.88.:
$3 50::
w:-Mo-m.:
n:-m-:o.c
3-93%:

5838.:
882.5:
:::wow::.:
8:: _ 58.:
3838.:

H. axon—

~50:
~50:

~50:
~50:
~50:
~50:
~50:
~50:
~50:

~50:
~50:
~50:
~50:
~50:
~50:

~50:
~50:
~50:
~50:
~50:

8h.

5
~:

3
w~
on
on
v:
m:
~:

can

~55
«.me
05mm
:2 N.
:03

3
w~
cm
on
em

2.8m

:

 

87

032
mamm—

0%:
0mm:
0%:
032
m.mmm_
mamm—
0%:

30:;
30:;
30:;
3:;
5:;

ES

8::
8:6

8::
corn:
8:6
8::
8::
8::
Siam

8::
00::
0:56
00::
8::
Stan

o~mfin
o~mfin
o~mfin
a~mfim
o~mfin

to:
who:

5. Kw~
v0;
h~v.~_
:0m
m:0~
vmo:
no:

5:
~_.m
e0~
mo.—
vow.—
n0—

2:.3
_:.m~
hcv.~_
m~:.:_
1::

52451— Sham—:2

 

vn0:
w~v:

~50—
.3:
m5:
2::
c _ :

:5:

mm:.:
3::
3::
~v:.:
v5:
moo:

m. _.m
5:.—
to:
hm»:
o~0:

uohmo><

a~3 58.:
8:35.:

5:238:
883 5.:
535::
88:5:
33:59:
a~3 38:
83.35::

amonmooo:
@0388:
$3 5:::
w:-m:-m.o
m:-m-oo.o
n:-mwmo0n

oooowmmo:
58:55:
:oomomoo:
8:: _ boo:
hocemmoo:

H. 5.09

~50:
~50:

~50:
~50:
~50:
~50:
~50:
~50:
~50:

~50:
~50:
~50:
~50:
~50:
~50:

~50:
~50:
~50:
~50:
~50:

3H.

w:
~:

#—
w~
on
:H.
V:
w:
~:

0::

~52
vwmm
ohmn
mg H.
:03

3
w~
cm
:H.

vw

macaw

-

 

88

hvwbm
www.cm
Smdm
Swen
Swdm
vadm
53.9”

www.cm
5%.me
hvwém
thdm
wadm
3:3

«:32
m.www_
m.wmm_
Ownn.
Wm?—

I33

cod
cod
cod
cod
cod
cad

OOOOOOO

cod

cod
cod
cod
cod
caN

OOOOOO

coN

Siam v: :m
0:5: wwwém
cows: 30.. _
Stun 9%:
82mm an$

52.53 8:832

Novd
94d
coin
as:
him

meN
:vd

fined
:3:
god
va.N
mth
Emu

50::
No;

Rm:
who:
omv:

notmo><

:ooooveo:
88:95:
88:08.:
88am 8:
$885.:
3938.:
:::::w:::

8828.:
::::m:::.:
:::wm:::.:
So: 50::
8mm 58.:
8858.:

8838.:
::::2. 5.:
:::mow::.:
8:: _ 58.:
$838.:

n. «:5

cow:
can:
:3:
can:
can:
can:
03.:

0%.:
com:
com:
com:
com:
:3:

Son:
Nam:
Sow:
Sow:
Sow:

8h.

3
mm
on
:5
v:
m:
N:

can

N2.—
vwmm
2mm
w: _ h
:00:

36 Egg:
aocofim

3
mm
cm
:5
vw

2.3m

 

89

VNHE
VNHE
was:
was:
3.53
3.5;
3.5..

3...:—
ems:
VNHE
v”;—
vmfiv—
VNSE

Sw.:m
$w.:m
hvw:m
hvwém
Sw.:m
hvw:m
5%.?”

ES

385
$86
mac:
386
mnood
nnood
nmood

38:
:25:
nmood
mm8.m
28:
30:5

:oN
:oN
:mN
:oN
:oN
:ad
:o.~

fizzle—~61— hOhhmmXfl—Z

www.mm
:86
«Ed
3N5
n:~.m
Eh:
:-.m

000000

880000

0

2.5.:
n a:
mom:
SN:
3N:
w _ m:
:N:

acm:
mom:
canon:
com:
_ 3.:
N _ N:

3.6

vm:.m
:3:
vmoN
mg:
mg:
EN

Sham—25V

:::::.v:::
:::::~m::
88:: 5.:
88mm 3.:
h::::: 5.:
:wmv _ 8::
:::::w:::

882 8.:
::::w:::.:
:::n~:::.:
5:: 58.:
:::~ 58.:
8838.:

:::::.v:::
:::::~m:.:
80:85.:
88:29:
h:::::_:.:
:mmzooo:
:::::w::.:

H. 3:5

:3:
:3:
:3:
:am:
:3:
:3:
:3:

:3:
£3:
:3:
:3:
:3:
:3:

:3:
:3:
:3:
:ow:
:ow:
:3:
:3:

mmh

3
”N
:n
:h
3
no
N. _

:3

«a:
mem
:3:
w: _ 5
::0w

3
wm
:w
on
V:
mo
N:

38:

:

 

9O

3:3V
mm.:ov
mm.::v
mm.:ov
mm.::v
3:3.
mm:av

2:3.
9:3»
2:3.
3:3
mm:ov
3:3

VNSE
VNHE
VNBE
VNSE
VNBE
VNHE
VNSE

ES

~_m_.m
~_m_.m
«SH—M
22:
~12:
~_m_.n
N_m_.m

N_m_.m
«.2:
~_m_.m
N_m_.m
~_m_.m
~_m_.m

mmood
mno:.m
nmood
mmocd
mm8.n
30:5
mmood

mmv. 3—
m8:—
ommé
:v:
mmm:
vowN— _
N.

000000

20.;
mvadm
02.2
mg:
mum:
mm:.w
5:5:

£2453 beam—.32

m _._
am:
:8:
5.:
:8:
N_v.nm

3::
:8:
3::
3::
8::
:8:

mmmd
new.—
:3.—
am:

min:
how:
3::

hohmo><

:::::v::.:
88:39:
8885:
88mm 5.:
h:::::_::
:wmv _ :8:
:::::w::.:

::8:. 8.:
::::m:::.:
:::w~:::.:
n:::_:::.:
8mm 58.:
8:: 58.:

:::::v:::
:::::~m:.:
88:: 5:
::::mm _ 0.:
h::::: _ :.:
:wmzaoo:
:::::w::.:

h. 5.5

:3:
:a:.:
:3:
:a:.:
:3:
:3:
:aw:

:am:
:3:
:3:
:3:
:ow:
:3:

:3:
:3:
:3:
:mm:
:3:
:3:
:3:

mmh.

3
mm
:m
on
g
m:
N:

:3

No:
wwmm
:3:
w: _ b
::$

1
mm
:m
:h
v:
no
N:

33m

NN

 

9|

N:
N:
N:
N:

N:
N:
N:
N:
N:

3:3
3:3.
32:3
3::
mm:ov
mm.:ov
3:3

I53

33—.
v22.
33..
33—.

v22.
33—.
3.2.
v2:-.
33—.
33—.

N_m_.m
.226
~_m_.m
~_m_.m
N_m_.m
N_m_.m
N_m_.m

52:83 seam—nu:

8::
mg:
m 3.:
:.::

000000

www.mv
wm.m~
:::.__
2m:
85.5
:0—
:3:

v:::
v:::
m::.:
v:::

32:
um _ .:
mm». _ .:
woo:
:8:
$8.:

3. _ .m
09:."
gm:
3 5.:
wow:
mvofim
mmv:

ucbmo><

5833.:
383:— :
2.33:— .:
8:533 .:

88:. 8.:
80:88.:
:::n~:::.:
5:: 58.:
83 5:0:
8:: 58.:

8893::
88:29:
8885.:
883 _ :.:
h:::::_::
:wmv _ :8:
:::::w::.:

a. SEQ

8::—
8::—
woo:—
woo:—

:3:
:0»:
:3:
:3:
:3:
:3:

:3:
:3:
:3:
:3:
:3:
:3:
:3:

8H.

:h
v:
m:
N:

:0:

no:
vwmm
:5:
w:_ n
:8:

:

$5 Eabzb
m cousc5>
3
mm
:m
on
v:
w:
2 _

2.3m

 

92

2.0.00
00060
00060
050.00

2.0.0.0
000.00
000.00
2.000
05060
00000

0.:
N:
N:

I53

«20—.
v:0_.
33—.
3.0..

v__0_.
v__0_.
v__0_.
v__0_.
v__0_.
v__0_.

36..
3.2.
.23..
36..
3.2.
33..
23..

33—.
«.20..
v__0_.

00— .0
.000
05.0

OOOOOO

m0mé~
0N0.N_
500
0:6
N0.»V
20.».
000.0

SEN—
90.0
000.0

521.55 Stags—2

000.0
30.0
000.0

OOOOOO

00.0

000..
30.0
000.0
N000
momd
010

00.0.0
00 _ .0
000.0

uohmo><

.00000NNO
50000000
2.0030— .0
002.0«3 .0

00000— 00.0
000000000
000000000
$005000
00mm 500.0
0000 _ 000.0

50000:.—
00002 000
0000330
5000030
50000000
2.0030— .0
002.003 .0

500003 .—
00000 _ 00 .0
00003030

.H 3.5

~000—
~000—
N000—
N000—

000.0
000.0
000.0
000.0
000.0
000.0

N000—
~000—
N000—
~000—
N000—
N000—
N000—

~000—
~000—
N000.

3.0

00.
Q0
00
N:

000

N0:
#09”
0000
00 _ h
0000

@—
0N
0m
0“.
V0
00
N:

3
0m
00

330

~N

 

93

000.00
000.00
000.00
000.00
000.00
000.00
000.00

000.00
000.00
000.00

83

v:0_. 03.00
v22. 000.:

”.20.. 30.0
33—. 000.».
3 _0_. 000.0
v__0_. 000.0
v22. ~00.0

v. _0_. 000.00
E _0_. 00¢.0
3 _0_. 000.—

EEIES batman—z

0 _ 0.0
000.—
000.0
30.0
V0.0

00v.0
00v.0

0 _ 0.0
00.0
V000

hobmo><

500003 ._
00000— 00.0
000000000
50000000
500000. .0
0000002 .0
0000003 .0

500003 .—
00000000
000000000

h. .500—

000.0—
000.0—
000.0—
000.0—
000.0—
000.0—
000.0—

000.0—
000.0—
000.2

8h.

v—
00
00
00
V0
00
0:

3
00
00

83m

 

94

00.0.0
00.0.0

00.0.0
00.0.0
00.0.0
00.0.0
00.0.0
00.0.0
00.0.0

00.0.0
00.0.0
00.0.0
00.0.0
00.0.0
00.0.0

83

000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0

00.
00.

.0000.
00.
00.
00.
00.
000.00.
00.

0.0
0.0
0..V
.000
00V
.0.V

53.1.51. 3:932

.0000
00V..

0000
000..
0..

0..

0..
0V0.00
.00..

V000.
0000..
000..
000..
V00.
000..

notme><

00V. .0000
00003000

000000 00.0
000000 .00
000000000
0000 . 000.0
000000000
00V. .0000
00003000

000000000
000000000
V000 .0000
00-m.0000.0
00-0.00000
00-m0000.0

H. 5.90—

0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

00
0..

V.
00
00
00
V0
00
0. .

000

000.
V000
0000
00. 0
0000

mm...

250$» 330m

.23.

 

0. 592.00%

95

0.00..
0.00..

0.00..
0.00..
0.00..
0.00..
0.00..
0.00..
0.00..

0.00..
0.00..
0.00..
0.00..
0.00..
0.00..

00.0.0
00.0.0
00.0.0
00.0.0
00.0.0

:3.

000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0

00.
000.00

000.0 . 0.
000.00V
V0000.
.0V.00
0V0.00
000..0V
0..0V

V0.0
0V0.0
V000
0V0.0
0V..0
V0.0

00.
00.
00.
00.
00.

5.21:3... 3:98.).

V0000
000.0

.000
.00..
000.0
000.0
00V0
000.00
VOV0

00.0
00.0
00.0
00.0
00.0
V0.0

000.0
VV0.0
V0. .

000. .
00V..

uctmo><

00V. .0000
00003000

000000000
000000.00
000000000
0000 . 000.0
000000000
00V. .0000
00003000

000000000
000000000
V000 .0000
00.m0000.0
00-m.0000.0
00-m0000 .0

000000000
000000.00
000000000
0000.0000
000000000

.0 3.00—

0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

mm...

00
0..

V.
00
00
00
V0
00
0..

000

000.
V000
0000
00.0
0000

V.
00
00
00

V0

2:5

:

 

96

..00.0
..00.0

..00.0
..00.0
..00.0
..00.0
..00.0
..00.0
..00.0

..00.0
..00.0
..00.0
..00.0
..00.0
..00.0

0.00..
0.00..
0.00..
0.00..
0.00..

I53

0000.0
0000.0

0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

000.0
000.0
000.0
000.0
000.0

52:53 batman—2

.0000
V00

000V. . .
000.000
000.00
0.0.00
000.00
V0.0V0
.000

OOOOOO

00V.V0_
000.00
V0000
000.00
000.00

.0000
00V.0

000.0
00.0

0V0.0
000.0
00.0
00.0

000.0
.000
0V0.0
000.0
0.0.0
V000

V000
000. .
000. .
000.0
000.0

uohmu><

00V. .0000
00003000

000000000
000000.00
000000000
0000. 000.0
000000000
00V. .0000
00003000

000000000
000000000
V000 .0000
00.0.0000
00-M0000.0
00-m.0000 .0

000000000
000000.00
000000000
0000 . 000.0
000000000

h. 5.09

0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

mm...

00
0..

V.
00
00
00
V0
00
0..

000

000.
V000
0000
00. 0
0000

V.
00
00
00

V0

2.30

00

 

97

00.00.
00.00.
00.00.
00.00.
00.00.
00.00.
00.00.

00.00.
00.00.
00.00.
00.00.
00.00.
00.00.

..00.0
..00.0
..00.0
..00.0

..00.0

83

0000.0
0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

0000.0
0000.0
0000.0
0000.0
0000.0

£2-83 Shana:

00.
00.
00.
00.
00.
00.
00.

OOOOOO

..0.V0.
.0000
..000
V00.V0
0V0.0

000.0
000..
0..

0..

0..
0V0.00
.00..

0000..
V00.
.00..
000..
000..
000..

000.0
00V..
000.0
0V0.0
00.0

uohmo><

00000V000
000000000
000000.00
000000.00
000000.00
000V.000.0
000000000

00000500
000000000
000000000
0000 .0000
0000 .0000
00005000

000000000
000000.00
000000000
0000.0000
000000000

h. 5.00—

000.0
000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0

0.00.0
0.00.0
0.00.0
0.00.0
0.00.0

3.0

V.
00
00
00
V0
00
0..

000

000.
V000
0000
00. 0
0000

250$» 330%
3002.00

V.
00
00
00
V0

2.80

 

98

00.0V0
00.0V0
00.0V0
00.0V0
00.0V0
00.0V0
00.0V0

00.00.
00.00.
00.00.
00.00.
00.00.
00.00.

53

0V00.0
0V00.0
0V00.0
0V00.0
0V00.0
0V00.0
0V00.0

0V00.0
0V00.0
0V00.0
0V00.0
0V00.0
0V00.0

0000.0
0000.0
0000.0
0000.0
0000.0
0000.0

Evan—51—

 

000.0 .0.
000.00V
V0000.
.0V.00
0V0.00
000..0V
0. .0V

OOOOOO

00.
00.
00.
00.
00.
00.

hgmlfl—Z

.000
.00..
000.0
000.0
00V.0
000.00
VOV.0

000.0
000.0
0000.0
00.0
.000
0000

39.930.

00000V000
000000000
000000.00
000000.00
000000.00
000V.000.0
000000000

00000.000
000000000
000000000
00005000
0000 .0000
0000 .0000

000000000
000000.00
000000.00
000000.00
000V.000.0
000000000

.0 5.00—

000.0
000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0

8,0

V.
00
00
00
V0
00
0..

000

000.
V000
0000
00 . 0
0000

00
00
00
V0
00
0. .

2.30

:

 

99

0.0V00
0.0V00
0.0V00
0.0V00
0.0V00
0.0V00

00.0V0
0.0V0 0
00.0V0
00.0V0
00.0V0
00.0V0
00.0V0
00.0V0

ES

0.0..0
0.0..0
0.0..0
0.0..0
0.0..0
0.0..0

0V00.0
0. 0. .0
0V00.0
0V00.0
0V00.0
0V00.0
0V00.0
0V00.0

5.2183 hob—MESS.

000.000
000.00
.000
000.00
V0.0V0
0

OOOOOO

00V.V0.
00.
000.00
V0000
000.00
000.00
00.
000.00

000. .
00V.0
000.0
00. .0
.V0.00

:00
000.0
000.0
. .000
000.0

V000

V000
000.0
000. .
000. .
000.0
000.0
V0000
000.0

uohmo><

000000000
000000.00
000000.00
000000.00
000V .0000
000000000

00000800
000000000
000000000
0000.0000
0000.0000
0000.0000

00009000
00003000
000000000
000000.00
000000.00
000000.00
000V.000.0
000000000

H. 5.0:

000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0

000.0
000.0
000.0
000.0
000.0
000.0
000.0
000.0

mm...

00
00
00
V0
00
0. .

000

000.
V00 0
0000
00. 0
0000

V.
V.
00
00
00
V0
00
0. .

2.30

00

 

lOO

BIBLIOGRAPHY

Alliare, PE. 1985. Basics of the Finite Element Methods: Solid Mechanics, Heat
Transfer, and Fluid Mechanics. W.C. Brown Publishers, Dubuque, Iowa.

Bathe, K. and EL. Wilson. 1976. Numerical Methods in Finite Element Analysis.
Pretence-Hall, Englewood Cliffs. New Jersey.

Boyce, E.W. and RC. Diprima. 1986. Elementary Differential Equations. Fourth Edition.
John Wiley & Sons, Inc. New York.

Carslaw, HS. and J.C.Jaeger. 1959. Conduction of Heat in Solids Second Edition.
Oxford University Press.

Churchill, RV. and J.W. Brown. 1987. Fourier Series and Boundary Value Problem.
Fourth Edition. McGraw-Hill Book Company, New York.

Cleland, AC. and R.L. Earle. 1984. Assessment of Freezing Time Prediction Method.
Journal of Food Science. Vol. 49(4): 1034-1042. Chicago, Illinois.

Fried, I. 1979. Numerical Solution of Differential Equations. Academic Press, New York.

Gear, G.W. 1971. Numerical Initial Value Problems in Ordinary Differential Equations.
Pretence-Hall, Englewood Cliffs. New Jersey.

Haghigi, K. and L.J. Segerlind. 1988. Modeling Simultaneous Heat and Mass Transfer in
an Isotropic Sphere-A Finite Element Approach. Transactions of the American
Society of Agricultural Engineers. Vol. 31(2):629-637.

Henrici, P. 1977. Error Propagation for Difference Methods. Krieger Pub. Co.
Huntington, New York.

Irudayaraj, J. 1991. Moving Evaporative Front Prediction Using The Finite Element
Method. Paper Presented at the 1991 International Winter Meeting of the
American Society of Agricultural Engineers, St. Joseph, MI.

Irudayaraj, J ., K. Haghigi, and R.L. Stroshine. 1990. Nonlinear Finite Element Analysis
of Coupled Heat and Mass Transfer Problems with an application to Timber
Drying. Drying Technology Vol. 8(4):731-749. New York.

Jaluria, Y. and K. Torrance. 1986. Computational Heat Transfer. Hemisphere Pub. Co.
New York.

101

Kreyszig, E. 1988. Advanced Engineering Mathematics. Sixth Edition. John Wiley &
Sons, Inc. New York.

Maadooliat, R. 1983. Element and Time Step Criteria for Solving Time Dependent Field
Problems Using the Finite Element Method. Ph.D. Dissertation, Michigan State
University.

Misra and Young. 1978. The Finite Element Approach for Solution of Transient Heat
Transfer in a Sphere. Paper Presented at 1978 International Meeting of the
American Society of Agricultural Engineers (No. 9107551). St Joseph, MI.

Mohtar, RH. 1994. Dynamic Time Step Estimates for Transient Field Problems. Ph.D.
Dissertation, Michigan State University.

Mohtar, RH. and L. J. Segerlind. 1992. Time Step Criteria for solving Unsteady
Engineering Field Problems. Paper No. 9107551. American Society of
Agricultural Engineers 0. St Joseph, MI.

Myers, GE. 1971. Analytical Methods in Conduction Heat Transfer. McGraw-Hill Book
Company, New York.

Myers, GE. 1977. The Critical Time Step for Finite Element Solution of Two
Dimensional Heat Conduction Transients. ASME Paper No. 77-WA/HT -14.
American Society of Mechanical Engineers, New York.

Narasimhan, TN. 1978. A Perspective in Numerical Analysis of the Diffusion Equation.
Advances in Water Resources. Vol. 1(3): 147-155.

Ortega, J. M. 1990. Numerical Analysis A Second Course. Reprinted by Society for
Industrial and Applied Mathematics. Part of the Classics in Applied Mathematics
Series (No. 3).

Ozisik, N. M. 1980. Heat Conduction. John Wiley & Sons, Inc. New York.

Patankar, S. V. 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishers
Company, New York.

Patankar, S. V. 1991. Computation of Conduction and Duet Flow Heat Transfer. A
Textbook Featuring the Computer Program CONDUCT. Innovative Research,
Inc. Maple Grove, Minnesota.

102

 

Peraire, J. J. Peiro, L. Forrnaggia, K. Morgan, and O. C. Zeinkiewicz. 1988. Finite
Element Euler Computation in Three Dimensions. International Journal for
Numerical Methods in Engineering Vol. 26:2135-2159.

Power D. 1987. Boundary Value Problems. Third edition. Harcourt Brace Javanavich,
Pub. Academic Press, New York.

Reddy, J. N. 1984. An Introduction to the Finite Element Method. McGraw-Hill Book
Company, New York.

Rushton, K. R. and L. M. Tomlinson. 1971. Digital Computer Solutions of Ground Water
Flow. Journal of Hydrology Vol. 12:339-362.

Segerlind, L. J. 1976. Applied Finite Element Analysis. John Wiley & Sons, Inc. New
York.

Segerlind, L. J. 1984. Applied Finite Element Analysis. Second Edition. John Wiley &
Sons, Inc. New York.

Segerlind, L. J. 1994-1998. Personal Communications.

Segerlind, L. J. and E. P. Scott, 1988. Selecting a Time Step Value for Numerical
Solution Involving Food Materials. Paper presented at 1998 International Summer
Meeting of the American Society of Agricultural Engineers (No. 88-6011). St
Joseph, MI.

Shih, T. M. 1984. Numerical Heat Transfer. Hemisphere Pub. Co. New York.

Smith, G. D. 1985. Numerical Solution of Partial Differential Equations: Finite
Difference Methods. Third Edition. Oxford Applied Mathematics and Computing
Science Series. Claredon Press, Oxford.

Stoer, J. and Bulivsch. 1980. Introduction to Numerical Analysis. Springer—Verlag, New
York.

Tan, W.T. 1995. Emperical Time Step Equations for Radial Field Problems. M.S.
Dissertation, Michigan State University.

Williams, W. E. 1980. Partial Differential Equations. Oxford Applied Mathematics and
Computing Science Series. Claredon Press, Oxford.

Zienkiewicz, O. C. 1971. The Finite Element Method in Engineering Sciences. McGraw-
Hill Book Co. London.

103

"Illlllllllllllllllll