‘ l .1 ?:i.:..r

. Lg: .
“Hr
..... .

 

. - 5. Q

.509.th n1... . ,

«I4 I\.!.. . ‘
thrills}?! ! . ‘ r.‘

,{IIIWUI “(Iv-v!
'. Ii..‘
'17:...- ‘I ll:

 

tut-4!:

Qalx till-5‘
it.
I}

w
m

\H

.3 trims; : :
Duh at . I}! swift-'10”!
RAVI;- Qn‘” .
. «H3: {3:1

4 if”? "NQ‘W

nmvg‘g”

awn-~99”

9'

3...?"

ntuwwunsvmnaom

H.
”J

Lil?! Li!"
#7.! :2 Il.O.:|I| I...

"M

‘ . . L I. LVA «(9“ ”.1150, II‘ V:
. g 1! {37.31.31an F413. . .. rm... I y .«w ‘ .

ti
‘ i
!
I

‘
\

 

 

 

 

Ti-ESIS

.93

o o

GAN STATE UENIV

IIIIIIIIIIIII lIIIIIIIIIIIIIIIIIIIIIIIIIII

1293 02061 1731

 

 

 

 

 

 

 

 

 

 

This is to certify that the

dissertation entitled

Time Step Criteria For An Extended Forward
Difference Method

presented by

Gihoun Kwon

has been accepted towards fulfillment
of the requirements for

 

 

PhD degree in Agricultural Engineering

 

Date November 22, 1999

 

MSU is an Affirmative Action/Equal Opportunity Institution 0-12771

 

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11/00 WM“

TIME STEP CRITERIA FOR AN EXTENDED
FORWARD DIFFERENCE METHOD
BY

Gihoun Kwon

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Department of Agricultural Engineering

Agricultural Engineering

1 999

ABSTRACT
TIME STEP CRITERIA FOR AN EXTENDED
FORWARD DIFFERENCE METHOD
By

Gihoun Kwon

Parabolic differential equations govern a variety of time dependent problems in
science and engineering. Many of the problems have physical parameters that
are a function of the primary variable. This makes the problem nonlinear and
necessitates a numerical solution. The objective of this study was to develop a
forward difference procedure usable with nonlinear problems and investigate its
time step requirements relative to stability and accuracy.

A forward difference method was developed assuming a finite element solution
in space and a system science approach in time. The solution procedure
requires a diagonal capacitance matrix [C] and a symmetric stiffness matrix [K]
that can be split into two triangular matrices, [K] = [K.,] + [K], where [K] = [K.,]'.
The solution procedure also uses one of Saulev's approximations to e". The

new forward difference method has the general form

[[C] +-;1IK.,I){¢}1., = [[C] - hiK]+-213[K.,]]{¢}j + hm,

where {do} is the vector of nodal values and h is the time step. This method is
explicit because the left hand side has an upper triangular form that can be

solved by back substitution. The stability criteria for this method is h < 40‘“,

which is twice the allowable time step associated with the traditional forward
difference method credited to Euler.

The time step requirements that produce an accurate solution were determined
by solving one— and two-dimensional problems using a range of physical
parameters and boundary conditions that produced a wide range of values for
Am. Reference values were calculated on a refined grid using small time steps
and a second order accurate central difference method. The one-dimensional
problems were solved on the unit length while the two-dimensional problems
were solved on the unit square. The time step criteria was defined using a one
percent error for the step change problem.

The a priori time step estimates for one-dimensional problems were:

11 node grid: h Am = 2.25

21 node grid: h Aw = 3.2

The a priori time step estimates for two-dimensional problems were:

11 by 11 grid: h Am = 1.33 Step change problem
h Am = 3 Derivative boundary conditions
21 by 21 grid: h Aw < 4 All problems

The a priori time step criteria were verified by solving unrelated one- and two-
dimensional problems. The one-dimensional problem was heat transfer through a
granite wall with heat input on one side and convection heat loss on the other side.
The two-dimensional problem was a solid circle with an eccentric circular hole that
contained a high temperature fluid. The a priori equations gave time step values

that produced numerical solutions with acceptable accuracy.

ACKNOWLEDGMENTS

I would like to give special thanks to Dr. Larry J. Segerlind who has helped
and encouraged me during the research and writing of this dissertation. My
major professor, Dr. Larry J. Segerlind. supplied the original motivation and
continued to generously give his time and valuable consultation during this work.
His guidance, support, and patience are essential for the completion of this
research. His knowledge. wisdom and smile had a great impact on me. Dr.
Segerlind was beyond the description. Thank you. Dr. Segeriind.

A great appreciation is extended to all members of my guidance
committee for their continual support and trust in me. Thanks Dr. Merva. Dr.
Hernandez and Dr. Martin. And i greatly appreciate Dr. Braits for his continual
support and suggestion, after leaving Michigan State University.

A special thanks goes to my family, my father. mother, and my wife,
Angela. for their continual understanding and encouragement throughout my

academic endeavor.

TABLE OF CONTENTS

Introduction ............................................... 1
Literature Review .......................................... 5
Theoretical Analysis ....................................... 16
3.1 General Integral Solution .............................. 19
3.2 Explicit Formulation ................................... 21
3.3 Numerical Stability and Oscillation Analysis ............... 24
3.4 Determination of Am ................................. 31

3.4.1 Forward Iteration Method ...................... 31

3.4.2 Fried's Method .............................. 33
Accuracy Criteria: One-Dimensional Problems ................... 34
4.1 Numerical Procedure ................................. 34

4.1.1 The Differential Equation ....................... 34

4.1.2 The Error Norm .............................. 38

4.1.3 Sampling Points .............................. 38
4.2 Experimental Results ................................. 42

4.2.1 Eleven Node Grid ............................. 42

4.2.2 Twenty-One Node Grid ........................ 45

4.2.3 Effect Of The Derivative Boundary Conditions ...... 45
4.3 A Priori TIme Step Equations ........................... 52
4.4 Validation Of One A Priori TIme Step Equation ............. 54
Accuracy Criteria: Two-Dimensional Problems ................... 59
5.1 Experimental Procedure ............................... 63
5.2 Experimental Results ................................. 65
5.3 A Priori Trme Step Estimates ........................... 70
5.4 Validation Of One A priori “i"Ime Step Equation .............. 71
Summary And Conclusions .................................. 77
References .............................................. 83

LIST OF TABLES

Table 4.1 The material coefficients with the corresponding values

of A, and km for 11 and 21 node grids ...................... 37

Table 4.2 Percentage average error for several values of a and

M, 11 node grid ........................................ 44
Table 4.3 Percentage average error for several values of a and

M, 21 node grid ........................................ 47
Table 4.4 Percentage average error for constant a and variable

derivative boundary condition M. 11 node grid ................ 49
Table 4.5 Percentage average error for constant or and variable

derivative boundary condition M, 21 node grid ................ 51
Table 4.6 Temperature distribution at four locations in a granite

wall, At = 0.02 hi‘ ....................................... 58
Table 5.1 The material coefficients with the corresponding values

of A, and A... used in the numerical experiments .............. 62
Table 5.2 Percentage average error for constant a and variable M.

11x11 node grid ........................................ 67
Table 5.3 Percentage average error for constant a and variable M.

21x21 node grid ........................................ 69
Table 5.4 Data for the 5x5 and 9x9 node grids ....................... 75 .

vi

Figure 3-1

Figure 4-1

Figure 4-2

Figure 4-3
Figure 4-4
Figure 4-5
Figure 4-6
Figure‘4-7

Figure 4-8

FIgure 4-9

Figure 4-10

Figure 5-1
Figure 5-2

Figure 5-3

Figure 5-4

LIST OF FIGURES

The stable and unstable regions for (3.50) ................. 30
Heat transfer in a one-dimensional insulated material

with heat convection on both boundaries .................. 36
Spatial and time domains .......... ' .................... 39
Coincident sampling point locations for the 11. 21. and

31 node grids ....................................... 40
Percentage average error as a function of time step,

11 node grid ........................................ 43
Percentage average error as a function of time step,

21 node grid ........................................ 46
Percentage average error or = 17.9 and variable M.

11 node grid ........................................ 48
Percentage average error or = 17.9 and variable M.

21 node grid ........................................ 50
Percentage average error for the step change problem ....... 53

Schematic representation of heat source and convection

heat loss boundary for a granite wall ..................... 55
Temperature distribution in a 9 cm thick granite slab

for numerical and analytical values. At = 0.02 hr ............ 57
Schematic representation of the two-dimensional region ...... 60

Sampling points for the 11x11 node grid for the unit square . . . 64

Percentage average error for a constant thermal diffusivity
of the material, or and variable M for 11x11 node grid ........ 66

Percentage average error for a constant thermal diffusivity
of the material. a and variable M for 21x21 node grid ........ 68

vii

Figure 5-5 The eccentric circles used for the verification calculations ..... 72

Figure 5-6 Location of the 15 sampling points used in the
verification calculations ................................ 74

viii

CHAPTER ONE

INTRODUCTION

A class of differential equations called parabolic equations governs a
variety of non-steady-state physical problems in the science and engineering
literature. Parabolic Equations are often referred to as diffusion equations and

have the general form

au
cEmV-(vu) (1.1)

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

the unknown variable: temperature, pressure head, moisture content, and so on.

Equation 1.1 is applicable to time dependent problems such as the diffusion of

neutrons or a gas into a solid, fluid flow and transport of solutes in a porous

media, and transient heat transfer. Examples of the diffusion process in
biological systems include but are not limited to:

1. The water flow equation for infiltration in the unsaturated soil region is
governed by parabolic type equations known as Richard's equations.

2. Salt movement and accumulation in the soil around micro-irrigation
emitters. This problem is important where the irrigation water contains salt
that moves through the soil and eventually accumulates at the fringe of the
wetted region. The same physical phenomena govern the movement of
solutes through a porous media and into the groundwater.

3. Transient heat transfer problems including the cooling of agricultural

products stored in bulk as well as the cooling of processed foods after
they have been placed in containers.

4. The drying process where heat is used to drive moisture out of a product.
This problem is governed by a pair of coupled partial differential equations
similar to (1.1) where each equation contains terms related to the
temperature and moisture.

Several engineering and mathematics books discuss the derivation of

(1.1) and give analytical solutions for simple shapes with mathematically "nice"

boundary conditions, Churchill (1987), Powers (1987), and Ozisik (1980). A

numerical solution of (1.1) is the only option for problems with irregular

boundaries and nonlinear material coefficients. The numerical solution of (1.1)

involves operating on the space terms using the finite difference method, the

finite element method, or a control volume approach. Each of these methods
converts the space-time problem into a system of ordinary differential equations

in time that has the general form
[01%LIKIIuI-IFHo} (1.2)

where [C] is a capacitance matrix coming from the transient term, [K] is a
stiffness matrix originating from the second order partial derivatives in space, and
{F} is a vector containing the known forcing terms such as boundary
temperatures and internal heat generation. The numerical procedure for
converting (1.1) into (1.2) is discussed in several books including Narasimhan
(1979), Patankar (1980), Segerlind (1984), and Smith (1985). The general

procedure for a finite difference solution of (1.1) is to simultaneously solve the

equation in space and time and skip the formulation defined by (1.2). Finite
difference solutions, however, can be written in a format identical to (1.2). The
finite element solution of (1.1) has the added advantage that [K] is symmetric for
many physical problems. The finite difference and control volume solutions to
these same problems often result in an unsymmetrical [K].

The most common procedure for solving (1.2) in time is to approximate
. . d{u} . . . .
the derivatIve term 71;— usrng finIte dIfferences and a srngle step procedure.

The general form of the single step procedure is (Segerlind, 1984)

(lC]+94t){U}. = (lCl—(1—9MthDIUl. +At((1-9){F}. +9{F}.) (1.3)
where 8 is a specified parameter. The three most common values for 8 are:

8 = 0: Euler's forward difference method

8 = The central difference method

Nl—x

8 = 1: The backward difference method
The subscripts "a" and "b" in (1.3) denote the values at the beginning and at the
end of the time step, respectively.

The forward difference procedure, 8 = 0, is ideal for solving nonlinear
problems because all of the required information is evaluated at the beginning of
the time step. The central difference scheme requires knowledge of the material
properties at both ends of the time step while the backward difference method
uses material properties at the end of the time step. The central and backward
difference methods require iteration calculations within a time step that increases

the computational time. The fonIvard difference method defined in (1.3) is known

to be conditionally stable and requires a small time step to obtain an accurate
solution, Smith (1985). Another difficulty with the fonIvard difference method is
that it is unclear whether the time step that satisfies the stability criterion also
gives an accurate solution.

Trujillo (1976) presented a procedure for solving (1.2) that involves
splitting a symmetrical [K] into upper and lower triangular matrices. Trujillo
obtained stable explicit algorithms similar to the central difference method. He
solved a single transient heat transfer problem with convection boundary
conditions to show that the new method gives results equivalent to the central
difference method.

This investigator has studied Trujillo's (1976) approach to solving (1.2) and
has found a forward difference method whose stability requirement allows a time
step twice the size of the traditional Euler's forward difference method, 8 = 0 in
(1.3). The objective of this dissertation is to present the mathematics of this
forward difference method along with an equation for selecting a time step value

that gives accurate results for one- and two-dimensional problems.

CHAPTER TWO

LITERATURE REVIEW

There are numerous published works on the numerical solution of partial
differential equations. Although several authors have published solution
procedures to linear and nonlinear time dependent problems, only a few authors
have discussed the determination of the time step value that produces accurate
results. The typical scenario is to present a numerical solution procedure and
compare it with an analytical solution of the partial differential equation using one
or more time step values. Those authors that have discussed the selection of a
time step value have done so in terms of a dimensionaless number that can not
be applied to composite bodies. This chapter contains a summary of the current
information on the solution of parabolic partial differential equations with the
primary emphasis on the determination of the time step value in (1.3) that
produces an accurate solution. The review of literature is limited primarily to
solutions using the finite element method because the theoretical development
given in Chapter 3 requires that [K] be symmetrical. Since this study is one part
of a continuing project started by Mohtar (1994), several of the references cited

here also appear in his dissertation.

2.1 Basic background information

The early finite element books contain brief discussions on using the finite
element method for the solution of time dependent field problems. Zienkiewicz
(1971) and Segerlind (1976) both discussed the numerical solution of the system
of ordinary differential equations but their discussions were confined to numerical
examples and they did not discuss any of the problems that can arise during the
solution process or compare different types of elements. Huebnes (1975) gave
the derivation of the capacitance matrix, [C] in (1.2), for transient heat transfer
but did not discuss the solution of the resulting system of differential equations.
Patankar (1980) 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 method in time due to its “friendliness" for all values of
time and solved the problem in space using the control volume method.

More recent books give the time dependent problem greater coverage.
Segerlind (1984) demonstrates some aspects of the physical reality and
numerical oscillation problems that can arise during the solution procedure.
Segerlind warns against the use of the quadratic elements because of physical
reality problems but he does not detail whether the numerical problems are long
term or short term or whether the errors are significant.

Shih (1984) discusses both the finite difference and finite element
numerical procedures and devotes a chapter to error analysis. Shih does not

discuss any estimates for At when solving time dependent problems. Shih also

devotes a chapter to the comparison of the finite difference and finite element

methods. His discussion covers smoothness of the basis 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."

Allaire (1985) discussed Euler's forward difference method for one-
dimensional problems in detail including the standard mathematical topics of
stability and numerical oscillations. Allaire does not discuss any solution in two-
or three-dimensions or make a comparison between linear and quadratic
elements in the one-dimensional case. Allaire used the Courant number as the
basis of comparison between solutions.

Jaluria and Torrence (1986) discussed the application of the three node
triangular element for solving heat transfer problems. They do not make any
comparisons with a two-dimensional finite difference solution or discuss any
other types of elements. The primary emphasis of their book was on the use of

the finite difference method to solve heat transfer problems.

2.2. Specification Of The Time Step Value

A major part of any transient analysis is to select a value for At in (1 .3) that
results in an accurate solution. A very limited amount of information appears to
be available on this topic.

Gear (1971) and Stoer and Bulirsch (1980) discuss the mathematical

approach used to estimate a time step value. They use the error between two

solutions obtained using time steps of At and 9% to determine an appropriate

step size. Their approach, however, does not give any information on how to

select the starting value for At. Another approach for the selection of a time step

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

Myers (1977) discussed the critical time step for two-dimensional transient

heat conduction problems. His discussion is based on using the maximum

eigenvalue, km, in the Euler's stability criterion of i or the nonoscillation

criterion .3)—- for the central difference method. Myers does not state whether

Am

these time step values gave accurate results. Dhatt and Touzot (1984) state that
the time step that eliminates stability and numerical oscillations may not produce
accurate calculations but they also do not give any suggestions of how to select
the time step value or show examples that supported their statement.

Shih (1984), Jaluria and Torrence (1986) and Segerlind (1984) avoid
explicit numerical evaluation of the time step. Each author discusses stability
and numerical oscillation problems but did not give a procedure for estimating the
time step as it relates to computational accuracy. Each author presented a

numerical solution procedure and compared it with an analytical solution using

one or more time steps. None of the authors commented on how they
determined the time step value that they used.

Rushton and Tomlinson (1971) used the alternation direction approach as a
numerical scheme. They studied stability and found that for different boundary
conditions, the Courant number, C, can be used to define accurate results. They
gave the following examples. 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 sudden change in discharge at a well, C should be less than
0.5. The authors indicated these were starting values and suggest a trial and
error procedure for selecting the optimal At value.

Henrici (1971) had an extensive discussion about the error propagation for
the difference methods in solving a partial differential equation. His theoretical
treatment did not include discussion of the time step for accurate results.

Wood and Lewis (1975) studied seven different finite difference time
marching schemes. They compared methods based on an accuracy criterion.
The authors related accuracy to oscillations and stability. They determined the
critical nonoscillatory time step for the, Crank-Nickolson method (C-N) using the
maximum eigenvalue. They showed, using numerical examples, that when the
time step value was increased beyond a critical value oscillations occurred.
According to their observation, inaccurate values appear in the backward
difference scheme for some time step values. They did not state whether
accuracy is a separate consideration in the numerical solution of parabolic

equations that needs needed to be addressed.

Reddy (1984) has a section on time dependent problems that is consistent
with much of the mathematics literature. He describes the stability in terms of the
roots of the characteristic equations and the eigenvalues of the global system.

Reddy includes a time step estimate for a structural dynamics problem,

At = 1m, where Tm", is the smallest period of natural vibration associated with
TI

the problem. Reddy also states that another estimate of 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

nondimensional form of (1 .1) in one dimension. Smith rearranged the difference

During the discussion of stability, Smith

 

equation and defined a term r =

At
(AxY'
stated that the explicit method is stable for r values less than 0.5. The implicit
Crank-Nickolson method has the advantage of being stable for all values of r.
Smith recommends r = 1 for an accurate solution for the Crank-Nickolson method
and discussed convergence and stability for some time stepping schemes.
Smith did not give a criterion for selecting the time step to obtain accurate results.
The r term defined by Smith was a dimensionless number but it has limited
application because the thermal diffusivity was assumed to be one.

Allaire (1985) called Smith’s r term the Courant number and included the
material properties that make up the thermal diffusivity. in addition to illustrating
stable and non-stable schemes, Allaire discussed an oscillatory stable scheme

that converges to the correct steady state values. Allaire showed the following

10

criteria to be true for the single step method.

0 < r s 0.25 No Oscillation
0.25 < r s 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 that reduces to the explicit
forward difference method and has a stability criterion of r < 0.5 when the
weighing coefficient is zero. He showed that the Crank-Nickolson method and
the fully implicit method 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 FC) < 0.5 for the
implicit method although lower values give better accuracy. One limitation when
using the Courant (Froude) number or the r term defined by Smith (1985) is that
the parameter does not include convection with the boundary conditions. Each
of the authors, Allaire (1985), Jaluria and Torrence (1986) and Smith (1985)
solved the heat equation using different boundary and recommended different
values of the Courant (Froude) numbers for accurate solutions.

Patankar (1991) presented a heat transfer computer program (CONDUCT).
This program uses the backward difference scheme to solve (1.3) in time, but
Patankar does not discuss the selection of the time step value and it is unclear
how or whether the selection is incorporated into the program.

Wood (1990) gives an extensive list of time-stepping schemes. He studied

stability, consistency, and oscillations and showed that these methods were

11

consistent with the analytical solutions. Wood advocated the use of time step
adjustment where the magnitude of the time step changes for a solution in time
but he didn’t give any formula for determining this time step.

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

Williams (1980) and Fried (1979) both studied the numerical solutions of a
partial differential equation 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 km (At)S 2

where kmax is maximum eigenvalue for the system of ordinary differential

equafions.

Haghighi and Segerlind (1988) solved the coupled heat and mass transfer
equations using the finite element method. They used Maadooliat’s (1983)
nonoscillation criteria to select a time step value and used elements that satisfied
the physical reality requirements that Segerlind (1984) discussed in his book.

Sarker and Kunze (1991) presented a finite element solution for
temperature distribution in a 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 the time step value in their paper.

lrudayaraj (1991) and lrudayaraj et. al. (1990) applied the finite element

method to the solution of a coupled heat and mass transfer problem. Both

12

papers used the stability criteria for selecting the time step, but there was no
check whether this time step ensured accurate results. The authors' calculated
results did not agree well with experimental data.

Liu et. al. (1984) used the stability criteria of a modified Runge-Kutta
method to solve a parabolic differential equation in time. Their work did not
discuss solution accuracy or whether the time step that satisfied the stability
requirement gave an accurate solution.

Peraire et. al. (1988) studied finite element solution of a fluid flow problem.

They used the Courant stability criteria of (At: k-rhgj to determine the time
u+

 

step. The variables in this equation are local sound speed, c, the average
element length, he, the velocity of fluid, u and a constant, k.

Alagusundaram, et. al. (1991) applied the 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,
but did not include numerical accuracy as one of them. They did not give the
time step value used to obtain the solution.

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 converge to a consistent value. They encountered a
stability problem and a physical reality violation that they call “jumping" at the
latent heat point. Although their solutions appear to be accurate, there is no

evaluation of a time step expression that could be translated to other problems.

13

(time to steady state)
100

 

Scott (1987) used an arbitrary accuracy criteria At =

to investigate the freezing and thawing of ice cream. Scott assumed that running
the problem for 100 time steps would be sufficient to obtain an accurate result.
Although this estimate might be a good starting point for some problems, there is
no mathematical justification for its use.

Maadoliat (1983) studied stability and physical reality oscillations of the
finite element numerical solution. He gave a set of conditions that must be
satisfied in order to avoid numerical problems such as numerical oscillations, but
he did not consider the accuracy criteria.

Logan (1992) used the maximum eigenvalue to define the time step. His

criteria is dependent on a value of B. When 8 < %, the largest time step, At, for

stability is based on

2

At= m (1.3)

where A"... is the largest eigenvalue. It is left to the user to choose the parameter
[3.

Bergan and Mollestad (1985) discussed algorithms for automatic
computation of time steps in discrete integration of dynamic problems. The
algorithm uses the current frequency, current period, and dynamic stiffness
parameter for automatic computation of time step. Other automatic time-step
control methods for dynamic systems have been described in the literature.

Zienkiewicz and Xie (1991) and Zeng et. al. (1992) proposed automatic selection

14

of the time step.

Orady (1992) discusses an automatic step length algorithm for stiff ordinary
differential equation systems. His algorithm adjusts an original time step by a
time variable factor that is a function of the truncation error and a safety factor. It
uses the Taylor expansion to estimate the truncation error. There was no
indication on how to choose the safety factor nor how to pick a starting time step.

Mohtar and Segerlind (1998) present an experimental approach to defining
a time step that integrate the two-dimensional field problem over a square region.
They describe their empirical equations as dynamic time-step estimates. The
calculation of the time step is based upon grid size and the smallest eigenvalue
of the system of ordinary differential equations, 9.1. Their time step estimates had
the general form

N”).,At = C (1.4)
where N was the number of nodes in the grid, M is the lowest eigenvalue in the
system of equations, and b and C were constants determined from experimental

studies.

15

CHAPTER THREE

THEORETICAL ANALYSIS

The objective of this research was to develop an explicit solution procedure
that can be used to solve nonlinear heat transfer and groundwater problems and
establish time step estimates that produce accurate results. The nonlinearity is
assumed to arise because the material coefficients are a function of the primary
variable. The idea for the procedure detailed here comes from the work of
Trujillo (1976).

The parabolic partial differential equation

_Iau

VZU — 3.1
a at ( )

with appropriate boundary conditions governs several physical problems
including conduction and convection heat transfer and groundwater flow. The
primary variable, temperature or pressure head, is denoted by U, and or
incorporates the physical properties of the problem. The quantity or can be a

function of U, or = f(U), making the problem nonlinear in time. The most general
form of VZU is

azu aZU azu
= + +

VZU
ax2 (By2 622

(3.2)

 

which is associated with the three-dimensional problem. The present discussion

is limited to one- and two- dimensions where (3.1) reduces to

16

__ = __ (3.3)

Of

62U a’U l 6U
6X2 +a—yT-E-aT (3.4)

 

respectively. Equations 3.3 and 3.4 assumes the material properties are not a
function of the coordinate directions.

Equations 3.3 and 3.4 can be converted into a system of ordinary
differential equations using either the finite difference method (OZISIK. 1994) or
the finite element method (Segerlind, 1984). The finite element method applies
the Galerkin method to the spatial coordinates at a fixed time. The element

contribution for Galerkin's method as given by Segerlind (1984) is
L 2 2
{RM} = -flW]T[a—q- ”iii/de (3.5)

The first term in the right side of (3.5) becomes the element stiffness matrix, [km].
The time dependent term in (3.5) can be evaluated using either of two different
formulations, the consistent formulation, or the lumped formulation. The
consistent approach uses a time derivative that varies linearly between nodal
points, while the lumped formulation uses a time derivative that is constant
between the midpoints of adjacent elements. The contribution from the time term

in (3.3) and (3.4) produces what is called the capacitance matrix which is

. _a_L 2 1
Ic"I— 6 I, 2] (3.6)

for the consistent formulation while the element contribution for the lumped

17

formulation is

m_g£1 0
kI—2I01I 07)

Summing the element matrices [k‘e’] and [cw] using the direct stiffness

procedure results in a system of first order differential equations defined by

[CIIUHKIIUI- {Fl= {o} (3.8)
Equation 3.8 can be premultiplied by [C]'1 and rearranged to obtain
IUI= lAllUl+ {F'} <39)
where
[A]: ——[Cl"lKl and {F'}= [Ci'lFl (3.10)
Equation 3.9 represents a system of linear differential equations that can be
solved several different ways. An integral solution patterned after Trujillo (1976)
is used here.
The mathematical and system science literature (Varga, 1962; Reid, 1983)

has many derivations where the authors substitute matrices for single variables

in mathematical operations. The only requirement is that the matrix operation be

. at . . at ,
defined. A common example IS e . The serIes expansron of e as a smgle

variable is

2 2 3 3
ea’=l+at+a£ +at + ...... (3-11)

 

 

The equivalent matrix operation ew‘ is

 

22 33
emufl+up+ugt+ML’+ ...... (an)

18

The following development uses this substitution approach. A single variable
solution is developed and the matrices are then substituted wherever the single

variable occurs.

3.1 General Integral Solution

A single equation equivalent to (3.9) is

dU
— = f -
t aU+ (313)

where a and f are coefficients, t is time and U is the primary variable. This
equation can be rearranged into

dU

(all—H) = dt (3.14)

Multiplying numerator and denominator by the coefficient 3 gives

(3% = dt (3.15)
A part of (3.15) has the form (ft—((3 where f(x) = g'(x). This allows the use of
integral relationship

[W W dX= ln((X) X) (3.16)

when f(x) = g'(x). Applying this integral relationship to (3.15) gives

adU '
S:I(___a U”) 1de (3.17)

Of

19

f

 

 

 

%[ ln(aU + r) lg” :11,“ (3.18)
Substituting the limits yields
éln[:go:ff]=(t—to) (3.19)
which can be rearranged to
IanLZTf) = a(t — 1,) (3.20)
and converted to an exponential form
aU + f : 94H") (3.21)
aUo + f
Another series of simple manipulations converts (3.21 ) into
U(t) = era' + 2(98’ — 1) (3.22)

asssuming to = 0.
Equation 3.22 is used to obtain the value of U(t) at a new time given the
initial value U0. Inserting (t+h) for t gives

U(t + h) = U(t)ea“""” + —;-(ea“*”’ — 1) (3.23)

Assuming f is constant in the interval (t, t+h) and starting the clock at zero for

each time step allows (3.23) to be rewritten as
ah f h
Um =Uje +;(e —1) (3.24)

where U,-.1 is the value at t = t + h and U,- is the value at time t. The series

expansion for ea, (3.11), allows (3.24) to be written as

20

 

 

2 2 2
UJ+I=UJ.[1+ah+a; + ...... ]+f[h+ah + ...... I (3.25)

The matrix equations equivalent to (3.24) and (3.25) for a system of

ordinary differential equations are
{1})”, = lube)” + Iakawh — 1) (3.26)

and

{u}, =[1+[A]h +%’32—+ ...... ){U}, +£h+ [AK-J. ...... ){F'} (3.27)
3.2 Explicit Formulation

Explicit solution methods in time are similar to Euler's method (Allaire,
1985) and require a diagonal capacitance matrix. The lumped formulation,
therefore, is used with this development. The finite element formulation for [K] is
used because this formulation always produces a symmetrical [K] even though
the physical problem may not have symmetry.

One way of deriving an explicit solution procedure is to split [A] into lower
and upper triangular matrices. The lower triangular matrix has the elements of
[A] that are below the main diagonal, zeros in the upper triangular region, and
each diagonal value is one-half of the original value. The upper triangular matrix
contains the second half of each element on the main diagonal, all the elements
above the main diagonal and zeros in the lower triangular region. This splitting
procedure produces a lower triangular matrix that is the transpose of the upper

triangular matrix and only one of the two matrices needs to be stored for

21

computer calculations.
The splitting of [A] gives
lAl=lLl+lUl (3.28)
where [L] is the lower triangular matrix and [U] is the upper triangular matrix. The
system of differential equations (3.9) becomes
{e} = (IL]+ lUIIl¢}+ {f} = [Llld>}+ lUll<bl+ {f} (3.29)
Since [U] is a common symbol for an upper triangular matrix, the symbol for the

primary or calculated variable is changed from U to (1). Define

g<t>=lUll¢l+lfl (3.30)
then

ld>1= [Ll{<v}+ 9(1) (3.31)

which has the same form as (3.9) with the general solution
{M0} = {1’0 )9“ + fiIeI‘I" -1) (3.32)

Solving (3.32) similar to (3.23) and (3.24) gives

I

{t}... = 6“" {¢}- + 2’46“)” -1) (3.33)
[L]
Inserting (3.30) into (3.33) produces

{I}... = efflil, .Mi‘flw —1) (3.34)

Using Pade’s approximation for e‘Llh (Varga, 1962; Trujillo,1976) as
-l
e’ILI" §[I+[i'2—yl) (l—E—z-h] for h>O (3.35)

An expansion of (3.35) produces

22

’3

‘-

[1+LEZIE)".II_I%I’1]=1_[L)I+(1LI’)2_(IL1”)3+ ......

Thus, the expansion agrees until the quadratic term for e‘L".

equation (3.34) with e'Lh results in

e-w(I...I=II,I+IUIILIII’II-e-WI

Substituting for e‘IL]h gives

2 2

{1_(1_lilfl).[l,fl

 

(“SI-(HEAR=1)» ,., III]...

Premultiplying (3.38) by [1 +LI'31E) PTOdUCGS
(I _%)(n,.. .-. (1 Alb—h... +h(lUl¢,- + ti.)
=I1+L95hjlilj+lUI hltl.+{f).h

(I —-[L2—]hII<I>},-n = (I +L112h+[U]h){¢}j + {fl/h

Premultiplying (3.39) by [C] yields

[IcI-Icltigjni... = [IcI+Ic1LIg-+Ic1uin]nr +[C]{f},-h

and using the definitions in (3.10) produces
[A] = -lC]"[K] = -[Cl“(th] + mm = -[C]"[KL] - [Cl"[Kul
Comparing the [A] in (3.43) with (3.28), gives
[L] = -lC]"[KLl and [U] = - [Cl"lKu]

Premultiplying (3.44) by [C], gives

23

-+{f})

(3.36)

Premultiplying

(3.37)

(3.38)

(3.39)

(3.40)

(3.41)

(3.42)

(3.43)

(3.44)

[CilL] = -ICIICI"IKu = IKLI (3.45a)
and
ICIIUI = -ICIICI"IKUI = [Ku] (3.45b)

Replacing (3.453, b) into (3.43),

[IcI+IK.I§I(¢,..I= [lei-thig-InivjinhIain (3.46)

A similar formula can be obtained by interchanging the [L] and [U] terms in

equations (3.30) and (3.31) resulting in
h h
IlC]+[K.]§]I¢,..I= [[Cl—lKulg—lKJhIlt-h{al.-h (34?)
Adding —[Ku]g+[Ku]g- to the right side, gives

IIcI+IKUIg]{I,..I-- [IcI-IK.I§-IK.h-IK.IQ+IK.I§]II,I+(th (3.48)

Collecting terms in the right side gives

[IcI+IK.I§]{¢,..}= [IcI-iK.I+IK.IIn+IK.I-’,1]{¢,}+{on (3.49)
And the final solution of (3.9) is

[IcI+IK.I§]{¢,.I-—- [IcI-IKh+iK.I-’,1]{¢,I+{qi,h (3.50)

Either of (3.46) or (3.50) can be used as a solution procedure in time.
3.3 Numerical Stability and Oscillation Analysis

The system of differential equations, (3.8), can be uncoupled using the

24

eigenvalues and eigenvectors of

(lKl— BlClXcD} = {0} (3.51)
to produce a new system of differential equations, Strang (1976),

{ZI-[AIZl— {F'}= {0} (3.53
where {Z} is the new variable and [7.] is a diagonal matrix where each diagonal

value is an eigenvalue. The eigenvalues of (3.8) are positive, real and unique

(Strang, 1976). An individual equation in (3.52) has the form

122+ 1.2.- _ f; = o (3.53)
dt

which is abbreviated here as

dZ .
—+7.Z-f =0 3.54
dt ( )

The remaining discussion centers around (3.54) rather than the more
complicated form given in (3.8).

The primary reason for uncoupling the system of ordinary differential
equations is because it simplifies the determination of the numerical stability
requirement for the solution procedure. If any one of the uncoupling first order
differential equations misbehaves, the entire system also misbehaves because
the final solution is a linear combination of all of the individual solutions.

The two problems that arise during a numerical solution of ordinary
differential equations are numerical stability and numerical oscillations. Stability
defines the behavior of the solution technique when a small error is introduced.
The error remains bounded for a stable technique. If the error grows and

becomes unbounded, the technique is said to be unstable and the results are

25

meaningless. Numerical oscillations occur in some solution techniques. The
calculated values on successive steps occur on opposite sides of the correct
solution when numerical oscillations exist.

The stability and oscillation characteristics of a solution technique can be
determined from the amplification factor. The uncoupled ordinary differential
equation has the standard form of (3.54) where it and f‘ are constants. The first

derivative term in (3.54) can be approximated by

d2 = 2,,l - z,

— 3.55
dt At ( )

where Z and Zn represent the values for Z on two successive time intervals.

Substitution of (3.55) into (3.53) and rearranging gives

[It I
or

Zi+l = (1’ Ami )zi 'I' (Aty
= A2, +f‘ (355)

where A = 1 - (A0)». is called the amplification factor.
Assuming the initial condition is Zn and writing a few steps of the solution
procedure in time gives
21 = AZ0 + f‘
and

Z2 =AZl +f’ =A(AZ0 +f‘)+f‘ =A2Z0 +(A+1)f‘

26

z, =Az2 +r' =A~iz0 +(A2 +A+1)f‘

Zn =A”Z0 +( "" +A"'3 + +A +1)f‘ (3.57)
Introduce a round-off error in Zo so that correct initial conditions is 2’0 =Zo +8.
The solution procedure now gives
2, =AZ, + f‘ =A£= Ag, + f‘

z, =AZl +f‘ = A2: +2120 +(A+1)f‘

Z, =AZ2 +f‘ =A3£+AZO +(A2 +A+1)'

Z, =AZ,_, +f' =A"E + AZo +( "" +A”’2 + +A +1)f‘ (3.58)
The error between the two solutions, (3.57) and (3.58) is
or, = Z, - 2,, = A": (3.59)
The amplification factor is to the nth power, therefore, it must satisfy
|A| <1 (3-60)
if the error is to remain bounded. When |A| > 1, the error grows and the

calculations become meaningless after a short period of time. Since the stability
of the system is not related to f', the forcing term is usually deleted during the
discussion of stability.

An amplification factor can be developed for (3.50) repeated here with {(1)}

replaced by {Z} and h is used to denote the step size At.

27

h h
(101+ $3.012)... = ([Cl—thl+ 51K. ]){Z}. + {al.h (3.61)
The stability analysis can be performed using a single uncoupled equation where

it
C = 1 and K = A, Ku = KL: —2—. The uncoupled equation for (3.61) is

OT

(1 + QEIIZIMF [I——3Z—"){z}, (3.62)

Rearranging (3.62),

 

 

—— {Z}, (3.63)
(. + f)
4
The amplification factor is
1- BL?»
_ 4
A m (3.64)
1+ —
4
which must satisfy -1< A < 1 or
1 -3123
4
-l < <1 (3.65)
I . I3:

The upper limit A < 1 is satisfied for all h > 0 since A is always positive. The

lower limit —1 < A is satisfied only if

28

h = AR; (3.66)

The smallest time step occurs when it = max, therefore, the stability criteria is

4
At< — (3.67)

IT) 8X

which is twice as large as the At obtained using Euler's forward difference

method (Allaire, 1985).
The numerical solution does not oscillate when the amplification factor, A, is

positive. This criteria is satisfied when the numerator of (3.65) is positive or
0 {titlija (3.68)

Solving this inequality gives

h < _4_ (3.69)
3;.

The largest eigenvalue also governs this criteria and the time step should satisfy

4
3).

max

At<

 

(3.70)

to avoid numerical oscillations. The stability and numerical oscillation criteria are

summarized in Figure 3.1.

29

 

Stable

 

No Oscillation

 

 

Oscillation

 

Unstable

 

Figure 3.1 The stable and unstable regions for (3.50)

 

3)»

max

30

>)I

3.4 Determination of km“

The time step for the explicit method developed in this chapter is related to
the largest eigenvalue, kmax, for the system of ordinary differential equations,
(3.8). The value of Am”, can be obtained two ways. It can be evaluated using the
vector iteration. method discussed by Bathe ,and leon (1976) or it can be
estimated using a theorem given by Fried (1979). The method outlined by Fried
(1979) is applicable only for the finite element formulation.

The evaluation of Max requires an iterative procedure that destories [K] and
[C]. These matrices must be rebuilt in order to solve the problem in time once
km... is known. The procedure given by Fried (1979) can be used while [K] and
[C] are being built but the estimated value for Am, exceeds the actual value.
Fried’s procedure, however, should reduce the computational time.

The vector iteration method given by Bathe and VIfilson (1976) was used in
this study because it was necessary to have an accurate value of max. A brief
review of the fowvard iteration method for determining km, is given here because
the procedure given by Bathe and Wilson is different than the procedure outlined

in most mathematics books.
3.4.1 FonIvard Iteration Method

Start with the system of differential equations in the form

[0161+ [Kl{¢} = {0} (3.71)

The forcing vector is deleted because it does not influence the value of 1...... The

31

system of differential equations (3.71) has roots 1., that satisfy

[K It») = iIdle} (3.72)
The capacitance matrix, [C], must be positive definite for the forward iteration
method. The method starts by assuming a unity eigenvector and iterates to a

solution. The procedure goes as follows

1. Define {x}1 as

 

 

2. Calculate

{YII : IKIIXII

3. Evaluate the following for n = 1, 2, .......

ICIIXIK.., = MK" (3.73)
[71% = NW... (374)

_ = {0.2.0}...
”Mm" {212.0}.

{YIKH
l (3.76)

(121;. in.)

(3.75)

 

 

{y}K+l =

The convergence proceeds as

{y}K+l ‘7 [KI‘IlbIn and p65“) —) Amax as n—)oo

provided {1)}:th at 0.

Bathe and Wilson (1976) give a set of calculations for the forward iteration

32

method using the matrices

r- - - m

 

 

5—410 2000
-4 6—41 0200
[K]: 1—4 6-4; [C]=0010
-01-45I [0001‘

 

 

The maximum eigenvalue is 10.63845 and the corresponding eigenvector is

1 0.10731“
0.25539
- 0.72827
_ 0.562274

I¢I

Ill

 

 

This example was used to check the software for evaluating kmax.

3.4.2. Fried’s Method

Fried (1979) proved that the largest eigenvalue for the system of ordinary
differential equations defined by (3.8) is less than the largest eigenvalue for the

individual elements. In equation form
AN 5 makai} 9 = 1,2, ..... (377)
where the element eigenvalue 15m is the largest eigenvalue of the element for
[K°]{<I>} = 1.7,[C°}{<b} (3.78)
Fried’s estimate for km. can be incorporated into finite element software and

1(8) and can be calculated as the individual elements are considered.

max

33

CHAPTER FOUR

ACCURACY CRITERIA: ONE-DIMENSIONAL PROBLEMS

The stability and oscillation criteria are important values in the
computational process but they do not necessarily define a time step that
produces accurate results (Dhatt and Touzot, 1984). A numerical procedure for
establishing the accuracy criteria for one-dimensional field problems solved using
(3.61) is discussed in the first part of this chapter followed by the definition of a
priori time step equations and a validation of one time step equation using a

unrelated problem.
4.1 Numerical Procedure

The investigation of accuracy using numerical experimentation includes a
partial differential equation, a solution procedure in time, an error norm, and an

error limit. These components are discussed herein.
4.1.1. The Differential Equation

The accuracy criteria developed in this chapter are for the one-

dimensional differential equation.

62¢ 69¢
__ = D _. 4.1
‘ ax2 ' at ( )

solved on the interval (0, H) with the initial condition
¢(x,0) =1 0 < x < H (4.2)

and the boundary conditions

At x = 0: :x—‘I = Me... —i>.> or ((0) = «I. (4.3)
_ . at .
At X ' H- _ a = M(‘I’bii - ¢r) 01' (NH) = ¢H (4-4)

The coefficients in (4.1), D)( and 0:, are material properties while M in (4.3) and
(4.4) incorporates a material property and the convection coefficient. The
quantities «to and on in (4.3) and (4.4) represent known values on the boundaries
while (boo and (lion in the derivative relationships are calculated values at the
boundaries. A schematic of the one-dimensional problem is given in Figure 4-1.
The partial differential equation (4.1) was converted to a system of
ordinary differential equations in time using the finite element method and a

lumped formulation. The system of differential equations was solved for several

different values of % (denoted by a) and M. The differential equation was also
t

solved for the step change problem with Ibo = din = 0 at the boundaries. The

values of % and M were varied to obtain a wide range of values for 1.1 and max.
t

The values for a and M and convection values and the corresponding values for
1.1 and kmax for 11 and 21 node grids are summarized in Table 4.1. The value M

= 0 corresponds to the step change problem. Several different M values were

used with or = 17.9 to investigate the influence of the derivative boundary

35

 

 

(

D

 

Insulated

al insulated material with convection

IITIBRSIOI'I
36

 

 

Figure 4-1 Heat transfer in a one-d
heat loss on both boundaries

Table 4.1 The material coefficients with the corresponding values of 1.1 and kmax

for 11 and 21 node grids

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11 nodes 21 nodes
or M A1
Kmax Kmax
17.9 0 0.55 21.8 87.7
1 7.9 5 0.29 23.6 90.2
1 7.9 1 8 0.45 34.2 105
17.9 25 0.47 41 .3 115
1 7.9 50 0.51 68.1 165
9.0 19 0.90 69.9 211
4.7 25 1.80 157 442
1.0 3.8 4.51 41 3 1600
0.5 3.8 9.02 826 3210
0.1 43 89.6 1 0800 27000

 

 

 

 

 

 

 

37

 

 

 

A)

CL)

 

conditions. All problems were solved using the initial conditions in (4.2).
4.1.2. The Error Norm

The error norm

il‘bcal - ¢refI
E, = l .. 100 (4'5)
N (¢max _ ¢min)

 

where N is the total number of sampling points was used in this investigation.
The value of N is five times the number of steps in time because the five
sampling points were used in space. Node 1 was not included in the evaluation
of the error even though it was a calculated value in several problems. The value
at node 1 was known for the step change problem. The error norm gives the

error as a percentage of the maximum difference in (j).

4.1.3 Sampling Points

The values calculated using (3.50) have to be compared with the reference
values at identical solution points in space and time, Figure 4-2. Grids with 11,
21 and 31 nodes were defined on the interval [0,1], Figure 4-3, to accomplish the
identical location requirement in space. Nodes 1, 4, 7, 10, 13 and 16 of the 31
node grid are coincident with nodes 1, 3, 5, 7, 9 and 11 of the 21 node grid and
with nodes 1, 2, 3, 4, 5 and 6 of the 11 node grid, Figure 4-3. The reference
values were obtained by solving (4.1) using a central difference scheme and the

31 node grid because this method is second order accurate (Gear, 1971). The

38

$4

 

Figure 42 Spatial and time domains

39

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 4 7 10 13 16
31 node grid

I 3 s 7 9 ll
21 node grid

I 2 3 4 s 6
11 node grid

0: Coincident sampling point

Figure 4-3. Coincident sampling point locations of 11, 21 and 31 nodes

40

values calculated using (3.61) were obtained for the 11 and 21 node grids.

The sampling points in time are a function of Xmax which is different for each
combination of a and M. Each solution had a unique set of sampling points
defined using the following procedure.

1. The values of M and AM were evaluated for the 11 or 21 node grid using the
designated values for or and M. These eigenvalues were obtained using a
separate computer program.

2. The time to steady state, T5, was calculated using

T... = {—9 (4.6)

1

This definition comes from e"" which is the first term in the analytical solution
of (4.1). Over 99% of the transient has been completed when A,t = 4.5 because
6'" = 0.0011. The sampling points in time were specifed in the first 60% of the
time to steady state. A study of several solutions indicated that the larger errors

occur early in the solution process.

3. A base time step was defined using

C
At:— (4.7)
1 4 . 4
where C = —, 1, —, 2, 3, and 4. The time step had to be less than —
3 3 km,

because this ratio defines the upper limit of a stable solution. The values of At

were rounded up or down to reasonable values.
4. The number of steps in time was calculated by dividing 60% of the time to

steady state by At and rounding the integer downward. The equation for the

41

number of steps in time was

0.60T 4 5 A 2.7A
N = 35 = 0.60 4 max = —-—m3‘- 4.8
At I4. I 1.. I CA. I ’

 

 

For example, when or - 1.9 and M = 0, 7.1 = 0.547, and in... = 21.8, the
number of steps in time for C = E:— is

_ 2.7(21.8) _
" (0.333)(0.547) “

 

(4.9)

which was rounded down to 320, a value dividable by two.
5. The reference values in time were calculated using the 31 node grid in space

and a time step one-half the value determined in step 3 above.
4.2 Experimental Results

The results of the numerical experiments are presented as a group of

tables and figures. The results are separated by the number of nodes in the grid.
4.2.1 Eleven Node Grid

The results for the 11 node grid for a mixture of or and M values are
presented in Figure 4-4. A summary of the values is also given in Table 4.2.

3

max

The largest error occurs in the step change problem up to At = The

 

A

average is below 1% for all combinations of or and M until the time step is

42

new moo: 3 d2..." 68: 65 Go 5:05: a ma cote manage oomEoEod .: 239....

 

 

000 6um NE
Hm Us. Hm Uh Uh Um
v m N v _ _
“ II o
\I I)

 

\ 1..

4
N

 

 

 

—.OId.nanIfl
m.cld.anI2IOl
.Id.wnI21.-.l
5.4-6.31211
o.old.o_IZI¢l .
9.218.”...le
Q: In .oIZIOl

4')

 

 

 

 

 

43

80883 SDWBAV $6

Table 4.2. Percentage average error for several values of a and M, 11 node grid

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

M 1. 1 1 4 2 3 < 4
or max _
31min lmax 311118.! lmax Amax lmax
At 0.015 0.045 0.060 0.090 0.13 0.18
17.9 0 21.8
°/oError 0.54 0.49 0.63 0.90 1.3 2.8
At 0.0097 0.029 0.038 0.058 0.088 0.11
1 7.9 18 34.2
%Error 0.14 0.30 0.39 0.63 1.3 2.9
At 0.0047 0.014 0.019 0.028 0.043 0.057
9 19 69.9
°/oError 0.14 0.29 0.39 0.61 1.3 4.3
At 0.0020 0.0060 0.0085 0.013 0.01 9 0.025
4.7 25 157
%Error 0.13 0.21 0.33 0.55 1.1 4.2
At 0.00080 0.0024 0.0032 0.0048 0.0073 0.0097
1 3.8 413
%Error 0.23 0.32 0.39 0.55 0.90 2.3
At 0.00040 0.0012 0.0016 0.0024 0.0036 0.0048
0.5 3.8 826
°/oError 0.20 0.32 0.39 0.55 0.90 2.3
At 0.000030 0.000090 0.00012 0.00018 0.00028 0.00037
0.1 43 10800
%Error 0.14 0.13 0.16 0.34 0.82 3.8

 

 

 

 

 

 

 

 

 

 

 

 

 

44

 

 

 

approximately 33—2—5.

4.2.2 Twenty-One Node Grid

The results for the 21 node grid for a mixture of a and M values are
presented in Figure 4-5. A summary of the numerical information is also given in
Table 4.3. The largest average percent error still occurs with the step change
problem although two other combinations of a and M have similar errors. The

average percent error is reduced relative to the 11 node grid and is below 1% for

all combinations of a and M until the time step is approximately 5—2.

4.2.3 Effect Of The Derivative Boundary Conditione

Several numerical solutions were performed keeping a, a constant (17.9)
and changing M. These numerical experiments were performed to determine
how the average error was related to M. The results for the 11 node grid is
presented in Figure 4-6 and summarized in Table 4.4. The results for the 21
node grid is presented in Figure 4-7 and summarized in Table 4.5. There are no
consistent trends within the two sets of results. The value of M = 5 produces a
long time to steady state and slowly changing 0 values but the error for this value
is larger than several of the other values of M. The error for the step change

problem is the largest or near the largest for both grids.

45

 

 

_.c|8.um|2|.l
WeldinIZIOl
_Io.unlz+
pendant—2+
93.6.2 uZlai
Q: Id.o— IEIII
0.: I6 .cIS—Iel

new one: _.N do.» as: 05 Co .6322 a mm 5:6 3805 83:69.8 .3. 9:9“.

 

 

 

 

 

 

 

 

 

uom 39W3AV '/o

46

Table 4.3. Percentage average error for several values of or and M, 21 node grid.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

M k 1 1 4 2 3 < 4
a max '—
31min max 31min max max max
At 0.0038 0.011 0.015 0.023 0.034 0.045
1 7.9 0 87.7
%Error 0.060 0.27 0.37 0.53 0.84 2.0
At 0.0030 0.0090 0.013 0.019 0.029 0.038
17.9 18 105
%Error 0.064 0.27 0.39 0.55 0.76 2.1
At 0.0015 0.0047 0.0063 0.0090 0.014 0.019
9 19 211
%Error 0.062 0.28 0.39 0.54 0.78 2.1
At 0.00070 0.0022 0.0030 0.0045 0.0068 0.0090
4.7 25 442
%Error 0.041 0.26 0.38 0.54 ' 0.76 2.1
At 0.00020 0.00060 0.00080 0.0012 0.0019 0.0025
1 3.8 1600
%Error 0.039 0.17 0.24 0.34 0.46 1 .0
At 0.00010 0.00030 0.00041 0.00060 0.00093 0.0012
0.5 3.8 3210
%Error 0.038 0.17 0.23 0.34 0.46 0.98
At 0.000010 0.000030 0.000049 0.000070 0.00011 0.00015
0.1 43 27000
%Error 0.014 0.14 0.30 0.44 0.65 1.7

 

 

 

 

 

47

 

 

Eco moo: E. .2 636:9 new Queue 6.. cote manage ommgcooaod .oIv 659...

 

 

 

 

 

 

 

 

 

ooofiwhmm—z:
ills.
v
o
4 w
.. x
N V
n
w
.5...
.\ m
A .. n 0
3|le 8
nuts—IT
fits—lei
«IS—Ill
cl—ZIOI “ . v
A
n

 

48

Table 4.4. Percentage average error for constant aand variable derivative
boundary condition M, 11 node grid

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

M A l l 4 2 3 < 4
(I max —
31min Amati 31min Amax [Imax Amax
At 0.015 0.045 0.060 0.090 0.13 0.18
1 7.9 0 21.8
% Error 0.54 0.49 0.63 0.90 1 .3 2.7
At 0.014 0.042 0.056 0.084 0.13 0.17
1 7.9 5 23.6
%Error 0.21 0.35 0.42 0.61 1.0 2.8
At 0.0097 0.029 0.039 0.058 0.087 0.1 2
17.9 18 34.2
%Error 0.14 0.30 0.40 0.63 1.2 4.1
At 0.0080 0.024 0.032 0.048 0.072 0.096
17.9 25 41.3
% Error 0.1 3 0.23 0.32 0.54 1 .1 3.9
At 0.0048 0.0015 0.020 0.029 0.044 0.058
1 7.9 50 68.1
%Error 0.30 0.28 0.1 5 0.29 0.73 2.9

 

 

49

 

 

 

Em ouoc _.w .5. osmtg new wine 5.. 22.6 389m 39:029. Niv 9:9“.

00. Cub.» NSF
........ i. U. U... H is
. a .s . 2
\\

 

IIelI

 

 

cnIZI-T
nuns—le
212141
«IS—Ill
al.—2+

 

 

 

 

 

HOWE! ZDWSAV '/o

60

Table 4.5. Percentage average error for constant or and variable derivative
boundary condition M, 21 node grid

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

M A I 1 4 2 3 < 4
O. max
31m” max 31min lmax lmax Amax
At 0.0038 0.011 0.015 0.023 0.034 0.046
17.9 0 87.7
°/o Error 0.060 0.26 0.37 0.53 0.85 2.3
At 0.0037 0.011 0.015 0.022 0.033 0.044
17.9 5 90.2
°/o Error 0.048 0.21 0.28 0.41 0.53 1.27
At 0.0031 0.0095 0.013 0.019 0.029 0.038
17.9 18 105
°/e Error 0.058 0.29 0.39 0.55 0.76 2.12
At 0.0028 0.0086 0.012 0.017 0.026 0.035
1 7.9 25 115
“/o Error 0.046 0.27 0.38 0.54 0.77 2.10
At 0.0020 0.0061 0.0081 0.012 0.018 0.024
17.9 50 165
%Error 0.013 0.18 0.27 0.43 0.62 1.60

 

 

51

 

 

 

4.3 A Priorl Time Step Equations

The hypothesis in this investigation was that an a priori time step estimate

has the general form

(Am... = C (4.10)
where C can be determined from the results of numerical experiments. The time

step chosen to determine C was the time step corresponding to the intercept for
an average error of one percent and the error for the step change problem. The

intercepts for the 11 and 21 node grids are presented in Figure 4-8. The

intercepts are $15- for the 11 node grid and X3_2_ for the 21 node grid,
max max
respectively. The a priori time step equations are:
11 node grid (Ammn = 2.25 (4-11)
21 node grid (404...... = 32 (4.12)

Both (4.11) and (4.12) give larger time step values than would be obtained
using the nonoscillation criterion, (3.70), of (A0)...” = 1.33 and also indicate that

values calculated using a time step near the stability criterion have more than a
1% error. A first analysis of (4.11) and (4.12) could lead one to conclude that a
larger time step can be used with the 21 node grid since 3.2 is greater than 2.25.

This is not true, however, because the value km... for the two grids are not
identical. Given a = 17.9 and the step change problem (M = 0), max = 21.8 for
the 11 node grid and 3...... = 87.7 for the 21 node grid. The respective time step

values are

52

EoBoa omcmco 69m 65 .8 5:6 manage Eooaon. .m-v 9:9“.

amhm m2:

XGE XQE

'&
'.<

 

 

 

 

 

 

 

mono: Fm mono: E.

 

 

80883 EOWBAV '/o

53

11 node r'd
9' At=3£§=0103
21.8

21 node 0
9” At = fl = 0.0365
87.7

A larger time step can be used with the 11 node grid.
4.4 Validation of One A Priori Time Step Equation

The a priori time step estimate (4.11) was used to solve a heat transfer
problem unrelated to the problems used to develop the equation. This problem is

illustrated schematically in Figure 4-9. It consists of 9 cm granite wall with a

surface heat flux of 100% on the left side, convection heat transfer of

100

 

m2 . on the right side, and an initial temperature is 20 °C. The thermal

conductivity of the wall is 1.6m—Vy6 where the therrnal diffusivity is 2.7*10‘3EV!.
- r

The problem defined in Figure 4-9 is presented and discussed by
Lienhard(1981). Lienhard solved the problem using the finite difference method

and compared the numerical calculations with the analytical solution.

 

.. 2Bi(k§ +Biz)cos[k.(§)1 -izii—‘i
T—T ——3—°— 1+Bi(1—3—]— L 3 “L2

’ ._ - 4.13
Initial L ":1 k3. (Bi + Biz + k:) ( )

L

where the kn are the roots of the equation kntankn = B. and B is the Biot number.
The analytical solution in (4.13) is credited to Lienhard(1981) since he does not

give a reference for it.

54

Granite Wall

 

 

 

—>
W
—> k=1.6—°—
q0=100—7 —>
m —>
-—> -3 W
0t=2.7x10 —
_: hr
1 10
I |
| I
0 3 6 9

q 0 = surface heat flux

hL = surface convection coefficient

hL =100

(I)...

w
m2 -°C

20 °C

 

Figure 4-9. Schematic representation of heat source and convection heat

loss boundary for a granite wall

55

The problem was converted to a system of ordinary differential equations
in time by apply the finite element method to a nine element, ten node grid. The
grid was uniform with each element 0.01 meter (1 cm) in length. The resulting
[C] and [K] matrices were analyzed to determine the minimum and maximum

eigenvalues which were 1.1 = 0.59 hr'1 and in... = 141 hr“, respectively. The
lowest eigenvalue was used to determine the time to steady state,

4.5 4.5
T = —- = = 7.63 hr
‘5 1., 0.59hr" (4°14)

while the value of 1...... was used in (4.11) to calculate a time step value

 

.25_ 2.25 _0016hr

7
At = " —
Am, 141m“l (4.15)

 

which was rounded up to 0.020 hr. This value is equivalent to 50 steps per hour.

The time step value of 0.020 hr is less than 1—4— = —4— : 0.028hrwhich is the

141

time step calculated using the stability criteria.

The problem was solved a second time using At = 0.009 hr. This time

step was below the oscillation criteria of At = fiL=000945hr A summagy of the

max

results at each surface, x = 0 and x = 0.090 m, and at two interval locations, x =
0.03 m and x = 0.06 m, are given in Table 4-6 and Table 4-7, respectively. The
results for the 11 nodes grid are displayed in Figure 4-10. The errors are
generally less than 0.10% and occur in the fourth digit. The a priori time step

equation provided accurate results for this example.

(It
0)

E mod u E .mo:_m>
.82sz new .motoEac Do.— nmfi 9:55 5:: £0 a m E 523566 2293th .34 9:9“.

2:9. .oEF

w m m m m NN YN EN m6 m6 N6 m.o 0.0 0.0

«r
I

 

mm

3,, 'BJmBJadwei

V
N

 

cog-om 66:95.2 IT _
5:28 623.92 I
II Ii II

 

  

 

 

 

mm

57

Table 4.6. Temperature distribution at four locations in a granite wall, At = 0.02 hr

Time,

Analytical Solution

Numerical Solution

 

 

 

 

 

 

 

 

 

 

 

 

 

 

hr X: 3 cm 6 cm 9 cm X= 3 cm 6 cm 9 cm
0 20.00 20.00 20.00 20.00 20.00 20.00 20.00 20.00
0.3 21.99 20.74 20.20 20.02 21.99 20.70 20.16 20.01
0.6 22.82 21.45 20.64 20.14 22.82 21.41 20.59 20.12
0.9 23.45 22.02 21.05 20.27 23.45 21.99 21.01 20.25
1.2 23.96 22.49 21.40 20.39 23.97 22.47 21.37 20.37
1.5 24.39 22.88 21.70 20.48 24.40 22.87 21.67 20.47
1.8 24.76 23.22 21.95 20.57 24.76 23.20 21.93 20.56
2.1 25.06 23.49 22.16 20.64 25.06 23.49 22.14 20.63
2.4 25.31 23.73 22.34 20.70 25.31 23.72 22.32 20.69
2.7 25.52 23.92 22.48 20.75 25.53 23.92 22.47 20.74
3.0 25.70 24.09 22.61 20.79 25.71 24.08 22.60 20.78
3.3 25.85 24.22 22.71 20.82 25.85 24.22 22.71 20.82
3.6 25.98 24.34 22.80 20.85 25.98 24.34 22.79 20.85

 

 

 

 

58

 

 

 

 

CHAPTER FIVE

ACCURACY CRITERIA: TWO-DIMENSIONAL PROBLEMS

The more common situation is to apply the forward difference procedure
developed in Chapter Three, (3.50), to two-dimensional problems. A priori time
step estimates for two-dimensional problems are defined and verified in this
chapter. The time step estimates are defined for the two-dimensional partial

equafion

2 2
D,§;§+Dy:7‘i’i=b,% (5.1)

where Dx, Dy, and D are material properties. These properties can be combined

into a single variable or 2 DB and (5.1) reduces to
i

when D)( = Dy = D. Equation 5.2 was used in this study along with the initial
condition

¢(x,y,0)=1 0<x<H, 0<y<H (5.3)
where H is the length of each side of a square region. A schematic of the two-

dimensional region is given in Figure 5-1 with the four sides labeled A, B, C, and

D.
The boundary conditions for (5.2) can be known values on any one of the

four sides or they can be derivative relationships. The possible boundary

59

 

A
Side C
H
8|: D / Side B
0 Side A H

 

Figure 5-1 Schematic representation of the two-dimensional region

60

conditions for (5.2) are

Side A: - 231 = Mot. - i.) or 30.0) = I. (54)
Side B: {E = Mii. — i.) or 901.10 = i. (55)
Side 0: 31: Mii. — t.) or ¢(x.HI= i. (5.6)
Side 0: “:9; = M0. —<I.> or 1(0.y)= i. (5.7)

The coefficient M in (5.4) through (5.7) incorporates material properties and a
convection coefficient. The quantity 0.. represents known values on the
boundaries that can be the same on each side. or different on each side or they
could be a function of the coordinate variable along each side. The quantity 0:, in
the derivative relationships are calculated values at the boundaries while (If is the
fluid temperature at a distance from the boundary. The calculated values for (iii,
are usually different at each node on a particular boundary.

The differential equation was solved for several different values of a and
M and for the step change problem with 0.. = 0 on each of the four sides. This
situation is denoted by M = 0. All problems were solved using the initial conditions
in (5.3). The values of or and M were varied to obtain a wide range of values for A,
and Am which are summarized in Table 5.1. The value of M = 0 corresponds to the

step change problem.

61

Table 5.1 The material coefficients with the corresponding values of 1.1 and max
used in the numerical experiments

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Grid M or A1 Amax Tss
0 0.0000005 0.0883 7.05 55.9
0.8 0.0000005 0.0225 7.16 200
2 0.0000005 0.0417 7.40 108
1 1X1 1
6 0.0000005 0.0654 9.70 68.8
10 0.0000005 0.0732 12.4 61.4
20 0.0000005 0.0802 19.4 56.1
0 0.0000005 0.0887 28.0 50.7
0.8 0.0000005 0.0225 28.5 200
2 0.0000005 0.0417 28.7 108
21 X21
6 0.0000005 0.0655 31 .5 68.7
10 0.0000005 0.0734 36.3 61.3
20 0.0000005 0.0805 49.7 55.9

 

 

 

 

 

 

62

 

5.1 Experimental Procedure

The accuracy criterion was developed by comparing solutions using (3.50)
to reference values at identical locations in space and time. Grids of 11x11,
21x21, and 31x31 nodes were defined on the unit square. The reference values
were obtained by solving (5.2) using the 31x31 grid and the second order
accurate central difference method (Gear, 1971). Equation 3.50 was used to
obtain nodal values for the 11x11 and 21x21 node grids. The error norm defined
in (4.9) was also used for the two-dimensional problems.

The space nodes were numbered from left to right starting along the x-axis
(y = 0). The sampling points were confined to one-eight of the grid as shown in
Figure 5-2. Figure 5-2 also shows the location of nodes 13, 14, 15, 16, 17, 25,
26, 27, 28, 37, 38, 39, 49, 50, and 61 for the 11x11 grid. Nodes 97, 100, 103,
106, 109, 193, 199, 202, 289, 292, 295, 385, 388, and 481 of the 31x31 node
grid were coincident with nodes of the 11x11 grid as well as nodes 45, 47, 49,
51, 53, 89, 91, 93, 95, 133, 135, 137, 177, 179, and 221 ofthe 21x21 grid.

The sampling points in time are a function of max that is different for each
combination of or and M. Each solution had a unique set of sampling points
defined using the following procedure. The value of 1.1 and kmax was evaluated
for the 11x11 or 21x21 node grid using the designated values for or and M.

These eigenvalues were obtained using a separate computer program that
generated [C] and [K] and used the forward and backward iteration procedures to

solve for M and Kmax (Bathe and Wilson, 1976). The time to steady, T55, was

63

 

 

. 2 Sampling points

Figure 52 Sampling points for the 11x11 node grid for the unit square

calculated using (4.5). The base time step was defined using

AIS.

 

(5.8)

max

where C =

 

, 1, g, 2, 3, and 4. The time step had to be less than

Cal—s

max

because this ratio defines the upper limit of a stable solution. The calculated
time step values, At, were rounded down to reasonable values. The

determination of the number of steps in time used in Chapter Four was also used

for the two-dimensional problem.

5.2 Experimental Results

The results for the 11x11 grid for a constant thermal diffusivity value of

a = 5x10.7 and a mixture of M values are presented in Figure 5-3. A summary of
the results are also given in Table 5.1. The largest error occurs for the step
change problem. The error for this problem appears to be significantly different
that the error for problems with derivative boundary conditions. The average

1.4

 

error is below one percent for the step change problem until At = k The
. 2.7

average error was below one percent for all nonzero values of M until At = T.
max

The results for the 21x21 grid are presented in Figure 5-4 along with a summary
in Table 5.2. The largest error also occurs with the step change problem but the
error is below one percent for all combinations of M including the step change

problem.

65

 

25 moo: Ex: .8 .2 oBmtm> new 6 6:89: 9: ho 33322—6 .mEEoE E8280 6 .8 .86 698on acceded .m-m 2:9“.

 

 

_o~ns_lel 9":le muEIXI ~u2+m.ou2|elou2+_

no...“ oEP

an XME

Is a
m m

rd E

. I
"<

 

 

L1

JOJJB BDBJBAV %

 

_ e8 . m n 8 2362.6 .8822» _

 

 

 

 

66

 

Table 5.2 Percentage average error for constant a and variable M.

11x11 node grid

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 1 4 2 3 4
a M 70..., 3.1 32 <

At 0.047 0.141 0.189 0.283 0.425 0.567
0.0000005 0 7.05

%Error 0.150 0.629 0.905 1.44 2.17 2.78

At 0.0450 0.135 0.180 0.270 0.405 0.540
0.0000005 0.8 7.16

%Error 0.00580 0.0648 0.124 0.286 0.579 0.900

At 0.0450 0.135 0.180 0.270 0.405 0.540
0.0000005 2 7.4

%Error 0.0131 0.137 0.251 0.542 1.04 1.53

At 0.0340 0.103 0.137 0.206 0.309 0.412
0.0000005 6 9.7

%Error 0.0189 0.181 0.313 0.642 1.19 1.72

At 0.0268 0.0800 0.107 0.161 0.241 0.322
0.0000005 10 12.4

%EI'TOT 0.015 0.155 0.263 0.532 0.989 1.46

At 0.0171 0.0510 0.0684 0.103 0.154 0.206
0.0000005 20 19.4

%Error 0.00596 0.100 0.167 0.327 0.615 0.937

 

 

67

 

 

 

 

 

 

Uta moo: merm .8 .2 03955 new a 5.99: 9: 20 Ezmztfi $855 E9800 m .6.— .otm mmmam>m Emoamn. .Ym 2:9“.

 

_ 8121 312+ mus—Ix! ~u2+m.ou2+on2+_

amfi we:

Kan: XNE rd:—

.2_« .22 222. «m A. «m

 

 

 

 

 

 

 

_ N.2 .. m u a 235.56 .959: _

 

 

 

 

Joua afieJaAv %

68

Table 5.3 Percentage average error for constant 01 and variable M.

21x21 node grid

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l l 4 2 3 4
a M A...“ <

31mm max Blmax lmax Amt): max
At 0.0119 0.0357 0.0476 0.0710 0.107 0.142

0.0000005 0 28
%Error 0.00440 0.0951 0.152 0.262 0.442 0.634
At 0.0116 0.0350 0.0467 0.0700 0.105 0.140

0.0000005 0.8 28.5
%Error 0.000160 0.00139 0.00451 0.0100 0.0640 0.140
At 0.0116 0.0340 0.0464 0.0690 0.104 ‘ 0.140

0.0000005 2 28.7
%Error 0.000150 0.00320 0.0101 0.0410 0.140 0.296
At 0.0105 0.0317 0.0422 0.0634 0.0950 0.127

0.0000005 6 31.5
%Error 0.000150 0.00574 0.0179 0.0723 0.232 0.472
At 0.00910 0.0275 0.0367 0.0550 0.0026 0.110

0.0000005 10 36.3
%Error 0.0000400 0.00466 0.0149 0.0625 0.209 0.420
At 0.00670 0.0201 0.0260 0.0400 0.0600 0.000

0.0000005 20 49.7
%Error 0.000640 0.00204 0.00679 0.0304 0.113 0.250

 

 

69

 

 

 

 

5.3 A Priorl Tlrne Step Equations

The a priori time step estimates are defined as the time step value that
produces a one percent error for the step change problem. The results are a
function of the grid size and the boundary conditions for the two-dimensional
problem. The following a priori equations can be applied to two-dimensional
problems.

11x11 node grid

Step change problem: (Athm z; (5.9)

Derivative boundary conditions: (A07.max = 2.7 (5.10)
21x21 node grid

All problems: (Ammax < 4 (5.11)

The a priori time step estimate for the 21x21 node grid indicates that
accurate results can be obtained using a time step value close to but below the
stability criterion. A time step value near the nonoscillation criteria should be
used when solving a step change problem with a 11x11 node grid. Field
problems with derivative boundary conditions take longer to go to equilibrium
then step change problems and a larger time step can be used when solving
these kind of problems. The time step for problems with derivative boundary
conditions can be approximately twice the size as the time step for the step
change problem. The time step used with a 21x21 node grid is actually less than

the time step used with a 11x11 node grid because km, increases as the number

of nodes increases.

70

5.4 Validation Of One A Priori Time Step Equation

The a priori time step estimate (5.9) was used to solve a heat transfer
problem unrelated to the problems used to develop the time step equation. This
problem is illustrated schematically in Figure 5-5. It consists of two eccentric
circles. A 2 cm diameter hole is enclosed in an 8 cm diameter circle. The hole is
centered on the vertical axis of symmetry and the bottom of the hole touches the
horizontal axis of symmetry of the larger circle. The inside of the hole is at 100

°C while the outside of the circle experiences convection heat loss with

W
cm2 -°

 

h = 1.5 and a fluid temperature of 0 °C. The thermal conductivity of the

cm2

SEC

The initial

 

solid is D: 0'753nw°—C while the thermal diffusivity was 56

condition is 100 °C.

Two grids were defined on the right side of the shape, Figure 5-6. The
first grid had five nodes along each of five radial lines extending from the center
of the hole. This grid had 25 nodes and is referred to as the 5x5 grid. The
second grid had nine nodes along each of nine radial lines. This grid had 81
nodes and is referred to as the 9x9 grid. Only the right side of the problem was
analyzed because of the vertical axis of symmetry.

The time-temperature history was calculated on the 5x5 grid using the
forward difference method defined by (3.50) while the reference values were
calculated using the 9x9 grid. Fifteen locations in space were used to compare

the forward difference method with the reference values. The fifteen locations

71

‘l

7

way

 

 

 

2
Thermal Diffusivity: a = 56%
SEC

'h=1.5—‘:—°
cm C

D = 0.75—W—
cm°C

0, =O°C

Boundary of the inner cylinder is at 100 °C

Diameter of the inner cylinder, 2 cm

Diameter of the outer cylinder, 8 cm

Figure 5-5 The eccentric circles used for the verification calculations

72

are defined by the shaded circles in Figure 5—6.
An eigenvalue analysis was performed for both grids and the results are
presented in Table 5.4. The a priori time step estimate (5.8) was used to estimate

the time step for the reference grid. Using (5.8) and the value for A“,

4 = 4 = 0.723 sec

3%, 3(1.84 sec'l)

At=

 

A time step value of 0.5 sec was used for the reference calculations which were
done using a second order accurate central difference method. Values were
calculated on the coarse grid using a time step value of one second. This value is
well below the stability criterion of 4A,", which gives 6.98 sec. The small value of
one second was used because of the need to compare calculated values at the
same locations in time.

The average error for the 15 sampling points in space and 773 sampling
points in time for the 9x9 gird and 1428 sampling points in time for the 5x5 grid was
1.45 percent. The calculated for the 5x5 grid are compared to the calculated results
for nodes, 9, 41, and 81 of the 9x9 grid in Figure 5-7. The two sets of calculated

values have a reasonable agreement.

73

 

 

 

 

 

 

 

Table 5.4 Data for the 5x5 and 9x9 node grids

Tu Am... A"...

 

714 sec 0.0058 sec" 0.573 sec"

 

0.0063 sec" 1.84 sec"

 

 

 

 

 

 

75

090:0 .0_:m0.._ 05 .o. 0.58 955800 00.5 .o. cow—5306 0.30.258. NA 0.39“.

 

a Sourcefiacm - I1. 850.2 85.0.20 0.20m“;

 

 

 

 

 

0000......
CV on em or o
o
. om
. ov
. om
. . . . . 8
£2 , .
- 0 q 1 l . o8

 

76

3.. ‘einieiadwe J.

CHAPTER SIX

SUMMARY AND CONCLUSIONS

The objective of this research was to (a) develop a forward difference
method that could be used to solve nonlinear problems governed by the

parabolic differential equation
cgthv-(w) (6.1)

and its boundary conditions and (b) establish a priori time step estimates for the
new fonlvard difference method that gives accurate results for one- and two-
dimensional problems.

Equation 6.1 was converted to a system of ordinary differential equations
in time using the finite element method with linear elements in space and a

lumped formulation in time. The resulting system of differential equations was

a
c—gh [Kl{¢}- {F} = {0} (6.2)
l
where [C] is a diagonal capacitance matrix, [K] is a symmetrical stiffness matrix,
{F} is the force vector, and {G is the vector of nodal values. Equation 6.2 was
solved in time using a procedure proposed by Trujillo (1976) that uses the

integral relationship

I; ’30: =ln(x) (6.3)

when f(x) = g'(x). one of Saulev's approximations for e“, and the fact that a

symmetric [K] can be split into a pair of triangular matrices, that is, [K] = [Ku] +

‘77

[K] where [KL] = [KU]T. Trujillo's procedure also requires that [C] be a diagonal
matrix.

The new forward difference procedure is
(0092114010,...=[1c1-2121+-‘-;£1K.1)10 +2110, (6.4)

where At is the time step. Equation 6.4 is an explicit solution procedure because

the left hand side has an upper triangular form that can be solved by back

substitution. A stability analysis of (6.4) shows that the time step must satisfy

 

4
At 5 —
1w (6.5)
for a stable solution and
4
At 5
31m (6.6)

to eliminate numerical oscillations. The stability requirement (6.5) allows a time
step value that is twice the size of the time step associated with Euler‘s fonlvard
difference method. The time step for the nonoscillation criterion is increased by
33 percent over the value for Euler's method.

The establishment of accuracy criteria using numerical experimentation
requires a partial differential equation with its boundary conditions, values
calculated using a solution procedure in time, reference values, an error norm,
sampling points, and the definition of an acceptable error level. The a priori time
step estimates were determined for the partial differential equation (6.1) solved

numerically in time using (6.4). The error norm was defined as

78

1 Zl‘i’cal - ¢refI

an
E =—
f N (¢max-¢min)

=1OO

 

which gives the error as a percentage of the maximum difference in 0. The a
priori time step estimates were assumed to have the general form

(mhm=c $0)
where C is determined using the intercept of a one percent error and the error
curve for the step change problem plotted as a function of Am... The step change

problem produced the maximum error or an error close to the maximum for all of

the one- and two-dimensional problems that were solved.

6.1 One-Dimensional Problems

Grids of 11, 21, 31 nodes were defined on a unit length along the x-axis.
The nodal locations for the 11 and the 21 node grids matched nodal locations in
the 31 node grid. The 11 and 21 node grids were used with (6.4) while the
reference values were calculated using the 31 node grid, a the second order
accurate central difference method, and a time step value one-half the size used
to solve (6.4). The material properties and the numerical values in the derivative
boundary conditions were varied to produce a range of Am. values from 21.8 to
27000. The sampling points in time were unique for each set of material
properties and boundary conditions because they were a function of Am”. The
sampling points in time were determined using a procedure that incorporated the

time to steady state.

79

The a priori time step estimates for the one-dimensional problem were

11 node grid (AGAM = 2.25 (6.9)

21 node grid (At)2,,,,, = 3.2 (6.10)
These equations indicate that the time step must be below the stability criteria of
(A070.1am = 4 in order to obtain an accurate solution but that the time step can be
greater than the value calculated using the nonoscillation criteria of (A025,... =
1.33. Equation 6.10 has a larger value for C but this does not mean that a larger
time step can be used with the 21 node grid because A“, increases as the grid is
refined.

The a priori estimate (6.9) was used to solve a time dependent heat
transfer problem involving in a granite wall that had a heat source on one side
and convection heat loss on the other side. The values calculated using (6.4)
agreed very well with values calculated using the analytical solution given by
Lienhard (1981). The two solution procedures produced values that matched in

the first three significant digits.

6.2 Two-Dimensional Grids

The a priori time step equations for two-dimensional problems were
developed using 11x11 and 21x21 node grids defined on the unit square. The
reference values were obtained using the second order accurate central
difference method and a 31x31 node grid. The sampling points in space were

limited to one-eight of the region because of the symmetry of the unit square.

80

 

The material properties and the numerical values in the derivative boundary
conditions were varied to produce a wide range of 3...... as was done with the
oneodimensional problems. The sampling points in time were a function of A..."
and were determined using the same procedure that was used for the one-
dimensional problems.

The largest error occurred for the step change problem for all
combinations of material properties and boundary conditions that were
considered. The step change problem had significantly larger error values for
the 11x11 node grid than those obtained for any derivative boundary condition,
therefore, two a priori time step estimates were defined for the 11x11 grid. The
errors for the 21x21 node grid were more closely grouped and a single a priori
equation was adequate for this grid. The a priori time step estimates for the two-

dimensional problems were:

1 1x11 node grid
Step change problem: (All... =1.33 (6.11)

Derivative boundary conditions: (Atym = 2.7 (6.12)

21x21 node grid
All problems: (Amm < 4 (6.13)

The a priori time step estimate for the 21x21 node grid indicates that
accurate results can be obtained using a time step value close to but below the
stability criterion. A time step value near the nonoscillation criteria should be
used when solving a step change problem with a 11x11 node grid. Field
problems with derivative boundary conditions take longer to go to equilibrium

then step change problems and a larger time step can be used when solving

81

these types of problems using an 11x11 grid. The time step for problems with
derivative boundary conditions can be approximately twice the size of the time
step for the step change problem.

The a priori time step estimate (6.12) was used to solve a two-dimensional
problem consisting of an eccentric circular hole contained inside a circular shape.
There were convection boundary conditions inside the hole and on the outside of
the circular shape. The finite element solution and (6.4) was solved using a 25
node grid while reference points calculated using a 81 node grid. Excellent
agreement was obtained between the two solutions for nodes on the surface of
the hole and for nodes on the outer surface. The a priori time step estimates
developed in this study provide excellent starting values for At when solving a
parabolic differential equation for one- and two-dimensional problems using the

fonlvard difference procedure defined in (6.4).

 

REFERENCES

Alagusundaram, K, D.S. Jayas, and WE. Muir. 1991. A Finite Element Model of
Three-Dimensional CO2 Diffusion in Grain Bins. ASAE Paper No. 91-6558.
American Society Of Agricultural Engineers, St. Joseph, MI.

Allaire, P. E. 1985. Basics of the Finite element Method. Solid Mechanics, Heat Transfer,
and Fluid Mechanics. Wm C. Brown Publishers, Dubuque, Iowa.

Bathe, K. and E. Wilson. 1976. Numerical Methods in Finite Element Analysis.
Prentice-Hall, Englewood Cliffs, NJ.

Bergan, G. and E. Mollestad. 1992. An automatic time-stepping algorimm for
dynamic problems. Comput. Methods Appl. Mech. Engng., 48:299-318.

Churchill, Ruel V. and James Ward Brown. 1987. Fourier Series andBoundary Value
Problems. Fourth edition. McGraw-Hill. New York.

Cleland, AC. and R.L. Earle. 1984. Assessment of Freezing Time Prediction
Methods. Journal of Food Science. 49: 1 034-1 042.

Dhatt, Gouri and Gilbert Touzot. 1984. The Finite Element Method Displayed. John
Wiley 8 Sons, London.

' Fried, Isaac. 1979. Numerical Solution of Difi’erential Equations. Academic Press,
New York.

Gear, G. Williams. 1971. Numerical Initial Value Problems in Ordinary Difl'erential
Equations. Prentice-Hall, Englewood Cliffs, NJ.

Haghighi, Kamyar and Larry Segerlind. 1988. Modeling Simultaneous Heat and
Mass Transfer in an Isotropic Sphere-A Finite Element Approach. Trans of
ASAE. 31 (2):629-637.

Henrici. Peter. 1977. Error Propagation for Difference Methods. Krieger Publishing
Company, Huntington, NY.

Huebner, Kenneth H. 1975. The Finite Element Method for Engineers. Wiley
lnterscience, New York.

lrudayaraj, J. 1991. Moving Evaporative Front Prediction Using the Finite Element
Method. ASAE Paper No. 91-7553. American Society of Agricultural

Engineering, St. Joseph, MI.

 

 

 

lrudayaraj, J., K Haghighi, and Fl.L. Stroshine. 1990. Nonlinear finite element
analysis of coupled heat and mass transfer problems with an application to
timber drying. Drying Technology, 8(4):731-749.

Jaluria, Y. and K Torrance. 1986. Computational Heat Transfer. Hemisphere
Publishing Company, New York.

Leinhard, John H. 1981. A Heat Transfer Textbook. Prentice-Hall, Englewood Cliffs,
NJ.

Liu, Kan W., T. Belytschko, and Y. Fei Shang. 1984. Partitioned rational Runge
Kutta for parabolic systems. Int. Journal for Numerical Methods in Engineering,
20:1581-1597.

Logan, D.L. 1992. A First Course in the Finite ElementMethod. PW S-KENT Publishing
Company, Boston, MA

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

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

Mohtar, RH. and L.J. Segerlind. 1998. Dynamic Time-Step Estimates for Two-
Dimensional Heat-Conduction Transients. Int. Journal for NumericalMethods

in Engineering, 42: 1 ~14.

Myers, GE. 1977. The Critical Time Step for Finite Element Solutions to Two
Dimensional Heat-Conduction Transients. ASME Paper 77-WA/HT—14. Am.
Soc. of Mechanical Engineers, New York.

Narasimhan, TN. 1978. A perspective on numerical analysis of the diffusion
equation. Advances in Water Resources, 1(3):147-155.

Orady, EA 1992. An automatic steplength control algorithm for stiff ODE's
systems. Comput. Struct., 43 (5):935-939.

Ortega, J. M. 1 990. NumericalAnaIysis: A Second Course. Reprinted by the Society for
Industrial and Applied Mathematics (SIAM), Philadelphia, PA.

Ozisik, Necati M. 1980. Heat Conduction. John Wiley 8. Sons, New York.

Ozisik, Necati M. 1994. Finite Difierence Methods in Heat Transfer. CRC Press. Boca
Flaton, FL.

Patankar, Suhas V. 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere
Publishing Company, New York.

Patankar, Suhas V. 1991 . Computation of Conduction and Duct F low Heat Transfer. A
textbook featuring the computer program CONDUCT. Innovative Research, lnc.,
Maple Grove, MN.

Peraire, J., J. Peiro, L. Formaggia, K Morgan and QC. Zienkiewicz. 1988. Finite
element Euler computations in three dimensions. Int. Journal for Numerical
Methods in Engineering, 26:2135-2159.

Powers, David. 1987. Boundary Value Problems, Third edition. Academic Press, New
York.

Reid, J. Gary. 1983. Linear Systems Fundamentals : Continuous and Discrete, Classic
and Modern. McGraw—Hill, NY.

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

Rushton, K F1. and L. M. Tomlinson. 1971 . Digital computer solutions of groundwater
flow. Journal of Hydrolog, 12:339-362.

Sarker, Nripendra N. and Otto Kunze. 1991. Finite Element Prediction of Grain
Temperature Changes in Storage Bins. ASAE Paper No. 91-6560. Amen'can
Society of Agricultural Engineers, St. Joseph, Ml.

Scott, E. P. 1 987. Simulation of Temperature and Quality Profiles in Frozen Foods Subject
to Step Changes in Storage Conditions. Unpublished Ph.D. Dissertation,
Michigan State University, East Lansing, MI.

Segerlind, L.J. 1 976. Applied Finite ElementAnalysis. John Wiley 8 Sons, New York.

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

Shih, Tlen-Mo. 1984. Numerical Heat Transfer. Hemisphere Publishing, New York.

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

Stoer, J. and R. Bulirsch. 1980. Introduction to NumericalAnalysis. Springer-Verlag,
New York.

 

Strang, Gilbert. 1976. LinearA lgebra and It 'sApplication. Academic Press, New York.

Trujillo, OM. 1976. An Unconditionally Stable, Explicit Algorithm for Finite Element
Heat Conduction Analysis. Nuclear Engineering and Design, 41:175-180.

Varga, RS. 1962. Matrix Iterative Analysis. Prentice-Hall, Englewood Cliffs, NJ.

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

Wilson, Edward L. and Robert E. Nickell. 1966. Application of the Finite Element
Method to Heat Conduction Analysis. Nuclear Engineering and Design, 4:276-
286.

Wood, W.L. and RW. Lewis. 1975. A comparison oftime marching schemes for the
transient heat conduction equation. lnternationalJournalforNumericalMethods
in Engineering, 9:679-689. '

Wood, W. L. 1990. Practical Time Stepping Schemes. Oxford Applied Mathematics and
Computing Science Series. Clarendon Press, Oxford.

Zeng. F., N.E. Wiberg, X.D. Li, and Y.M. Xie. 1992. A Posteriori Local Error
Estimation And Adaptive Time-Stepping For Newmark Integration In
Dynamic Analysis. Earthquake Engng, Struct. Dynamitx, 22:338-351.

Zienkiewicz, 0.0. and Y.M. Xie. 1991. A simple error estimator and adaptive time
stepping procedure for dynamic analysis. Earthquake Engng. Struct. Dynamics,
20:871-887.

Zienkiewicz, 0.0. 1971. The Finite Element Method in Engineering Science. McGraw-
Hill, New York.

86