LIBRARY
Mlchlgan State
Universlty
PLACE ll REI'URN Boxmmwombdnckomm yourrocotd.
TO AVON) FINES mum on orbdonddoduo.
r—————‘—_—————————————W
DATE DUE DATE DUE DATE DUE
[—T—T—j
MSUIIMWWVEWOWIW
DYNAMIC TIME STEP ESTIMATES FOR
TRANSIENT FIELD PROBLEMS
By
Rabi Hassan Mohtar
A DISSERTATION
submitted to
MICHIGAN STATE UNIVERSITY
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Agricultural Technology and Systems Management
Department of Agricultural Engineering
1994
ABSTRACT
DYNAMIC TIME STEP ESTIMATES FOR
TRANSIENT FIELD PROBLEMS
By
Rabi Hassan Mohtar
Parabolic equations govern a variety of time dependent problems in science and
engineering. Space integration of the ﬁeld equation produces a system of ordinary
differential equations. A common problem during the numerical solution of these
equations is specifying a time step that is small enough for accurate and stable results and
still reasonable for economic computations. This is a challenging mathematical problem
that is yet to be solved.
This study presents an experimental approach to estimate the time step that
integrates the ﬁeld equation within 5 % of the exact solution. Time step estimates were
determined and evaluated for one-dimensional and two-dimensional linear square ele-
ments. A section is also included on the effect of shape and size of two-dimensional
quadrilateral elements on numerical stability.
Time step estimates were determined for four numerical schemes for one-dimen-
sional problems and three ﬁnite element, and three ﬁnite difference schemes for two-
dimensional problems. Comparisons between ﬁnite element and ﬁnite differences as well
as correlations with the Froude number were conducted. The time step estimates are
functions of grid size and the smallest eigenvalue M. For some problem an initial mesh
is required to evaluate M.
The results indicates that the ﬁnite element and ﬁnite difference methods are
identical in one-dimensional problems and have similar time step estimates. These two
methods had similar accuracy two-dimensional problems. The central difference schemes
were superior to the other schemes as far as the ﬂexibility in allowing a larger time step
while maintaining good accuracy. Backward difference and forward difference schemes
were very close in their accuracy; the difference between the two schemes was attributed
to the stability requirement of the forward difference scheme.
The dynamic time step estimate for the central difference method was:
A, A N"18 = 1.13 for the one-dimensional problems, and A, A N055 = 1.6 for two—
dimensional problems, where N is the number of nodes in the domain.
The Froude number equivalence is deﬁned as the Froude number that gives a time
step equal to the one computed using the presented regression equations. For the range
of problems presented in this work and for the central difference scheme the Froude
number equivalence ranges from 0.5 to 2.7 in the one-dimensional problems and 0.46
to 9.13 in two-dimensional problems.
Copyright by
Rabilkwsm1hdohmu
1993
DEDICATION
This work is dedicated to a person that gave and was never given, to my mother who is
still around me praying for my success and savior.
To the soul of
WADAD SAID MOHTAR
To her memory I write this arabic poem
azalﬁtuuan.
éED‘DJHJaémJS‘WD.‘
Jagchlluﬂldamlul
éﬂkéésoéhwc
Wdeiié‘ach
wmajuéz
3e)!“ W J 3ft” 11-51.)?
gal u an. 3J5...» “J t,
1.1., 02.3.33 L. any“;
gauze-2 £13:- 31 L. J55
9-3; uh“; 4,313.49 'wuu' ' ‘1
This is also dedicated to my beloved father Hassan and the rest of my dear family
Sami, Samia, Asia, Salma, Nadia, Ramzi, Rifaat, Jihan, Lina, Hazar, David, Yasmeen,
and Samer.
Last but not least, a dedication to all who are suffering from oppression all over
the world particularly in the Arab nation.
ACKNOWLEDGMENTS
A special thanks to all of the smiling faces at MSU. A special appreciation to
those who help in ways they do not know about. To all of those I say keep it up folks.
A very special thank to Dr. Larry Segerlind for being a guide and a friend. Dr.
Segerlind was more than that. His exceptional cheerful personality is infectious to those
who work with him. His knowledge, wisdom and smile had a great impact on me.
Thank you Dr. Segerlind.
A great appreciation to all members of my guidance committee for their continual
support and trust in me. Thanks goes to Dr. Bralts for the time he spent advising me,
for his ﬁnancial support during the ﬁrst three years of my studies at MSU and for his
suggestions to ground my work to the agricultural areas. Thanks to Dr. Merva for
believing in me and for his teachings and the time spent correcting the dissertation. A
special thanks to Dr. Kunze for the long hours he spent with me teaching and researching
soil and water problems and for his in-depth comments on my research.
Appreciation to all friends, faculty and staff at the Agricultural Engineering
department especially Dr. von Bemuth, Lilian, Clara and all GSAG folks.
I would like to show gratitude to Gladis for her friendship and her successful
effort in keeping my faith in God healthy.
A special thank you to the I-Iarriri Foundation for providing me with ﬁnancial and
moral support to come to the USA to continue my graduate work. Without their support
this was not possible.
A list of those that need to be thanked ﬁlls many pages, to all of you friends and
family a warm thank you, you are all in my heart.
vi
TABLE OF CONTENTS
CHAPTER ONE (INTRODUCTION)
Introduction
Application to the Agricultural Problems
CHAPTER Two (JUSTIFICATION)
Justiﬁcation
CHAPTER THREE (OBJECTIVES)
Objectives
CHAPTER FOUR (ONE-DIMENSIONAL PROBLEMS)
Analytical Solution
Numerical Solution
Methodology
Time Step Estimate
Evaluation of the Time Step Estimates
Effect of Element Size on Accuracy
Comparison With Froude Number
New Scheme
CHAPTER FIVE (TWO-DIMENSIONAL SQUARE GRID)
Analytical solution
Numerical solution
Methodology
Finite Element and Finite Difference
Accuracy Ratio
Time Step Estimation Process
Results
Discussion
Finite Element and Finite Difference
Sensitivity of the Error to Time Step
Effect of Element Size on Accuracy
Comparison With Froude Number
Evaluation of the Time Step Estimates
CHAPTER SIX (THE EFFECT OF SHAPE AND SIZE ON STABILITY ...)
Mathematical Background
Methodology
Results and Discussion
CHAPTER SEVEN (SUMMARY AND CONCLUSION)
Recommendations
Future Research
REFERENCES
vii
page 1
page 1
page 3
page 3
page 5
page 16
page 16
page 20
page 22
page 24
page 25
page 29
page 38
page 46
page 49
page 53
page 56
page 57
page 5 8
page 59
page 60
page 62
page 68
page 69
page 79
page 79
page 86
page 87
page 87
page 92
page 106
page 106
pagelll
page 115
page 128
page 134
page 134
page 135
LIST OF TABLES
Table 4-1 Froude number equivalence for various A, values. page 54
Table 5-1 Number of eigenvalues and degrees of freedom for the various grid sizes
used in the derivative boundary condition problem page 97
viii
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-11
Figure 4-10
Figure 5-1
Figure 5-2
Figure 5-3
LIST OF FIGURES
Graphical representation of boundary and initial condition for the step
change problem. page 27
Variation in the error ratio as a function of the time step and the number
of nodes for the central difference scheme. page 31
Actual and regression curves for the data points that gave accurate numeri-
cal results as good as ANODE for the central difference
scheme. page 33
Actual curves for the data points that gave accurate numerical results as
good as ANODE all schemes. page 35
A logarithmic plot representation of the data shown in
Figure 4-4. page 36
Performance of the estimated time step for the CD scheme and N=16 grid
problem for three different problems. page 44
Performance of the estimated time step for the derivative boundary condi-
tion 2 problem for the various schemes and N =16 grid. page 45
Effect of element size on accuracy represented by the variation in error
ratio as a function of number of nodes. page 48
Comparison between estimated time step and Froude Number for the
Galean scheme. page 51
Froude number equivalence changes with A. page 52
Calculated time step performance for the step change problem and theta
= 0.333. page 55
Changes of the error as a function of time step for GCD scheme and a
grid size corresponding to 1089 total number of nodes. page 70
Changes of the error as a function of time step for various mesh sizes
under GFD scheme. page 71
Changes of the error as a function of time step for various mesh sizes
under GCD scheme. page 72
ix
Figure 5-4
Figure 5-5
Figure 5-6
Figure 5-7
Figure 5-8
Figure 5-9
Figure 5-10
Figure 5-11
Figure 5-12
Figure 5-13
Figure 5-14
Figure 5-15
Figure 5-16
Figure 5-17
Changes of the error as a function of time step for various mesh sizes
under GBD scheme. page 73
Changes of the error as a function of time step for various mesh sizes
under FFD scheme. page 74
Changes of the error as a function of time step for various mesh sizes
under CFD scheme. page 75
Changes of the error as a function of time step for various mesh sizes
under BFD scheme. page 76
Computed optimal time step as a function of N for the six
schemes. page 77
Computed optimal time step changes as a function of the number of nodes
for GFD and FFD schemes. page 81
Computed optimal time Step changes as a function of the number of nodes
for GCD and CFD schemes. page 82
Computed optimal time step changes as a function of the number of nodes
for GED and BFD schemes. page 83
Computed time step variations with the number of nodes for the three
ﬁnite element schemes. page 84
Computed time step variations with the number of nodes for the three
ﬁnite difference schemes. page 85
Computed time step variations with the number of nodes for GFD and
GBD schemes. page 88
Computed time step variations with the number of nodes for FFD and
BFD schemes. page 89
Finding a Froude number equivalence through best ﬁtting the GCD time
step estimate equation and the Fr. time step estimate for a ﬁxed boundary
condition and material properties i.e. ﬁxed A]. page 91
Froude number equivalence changes as a function of A, for the six
schemes. page 93
Figure 5-18
Figure 5-19
Figure 5-20
Figure 5-21
Figure 5-22
Figure 5-23
Figure 5-24
Figure 6-1
Figure 6-2
Figure 6-3
Figure 6-4
Figure 6-5
Figure 6-6
Figure 6—7
Figure 6-8
Figure 6-9
Figure 6-10
Froude number equivalence changes for the different schemes and
different values of A,. page 94
Error propagation with the time step for GFD scheme. Time steps are
reported as ratios of actual time step to the estimated one for the speciﬁc
scheme and grid used. page 98
Error propagation with the time step for FFD scheme. page 99
Error propagation with the time step for GCD scheme. page 100
Error propagation with the time step for CFD scheme. page 101
Error propagation with the time step for GBD scheme. page 102
Error propagation with the time step for BFD scheme. page 103
Shape parameter deﬁnition. page 113
Shape parameters values for some selected quadrilaterals. page 114
Effect of area on the largest eigenvalue for a square element. page 116
Effect of area on the largest eigenvalue for a trapezoidal
element. page 117
Effect of shape on the largest eigenvalue for a unit area. page 119
Effect of shape on the largest eigenvalue for an area of four. page 120
Comparative chart for some selected shapes as to their maximum
eigenvalue. page 121
Practical example 1 : dividing a rectangle into squares and
triangles. page 123
Practical example 2 : dividing a trapezoid into triangles. page 124
Comparison between a quadrilateral and a triangular element. page 127
xi
CHAPTER ONE
INTRODUCTION
A variety of time-dependent problems in the science and engineering ﬁelds are
governed by a class of equations called Parabolic Equations. In the engineering
literature they are often referred to as Diﬂ‘icsion Equations. These equations have the
general form:
6U
_ = V. (1-1)
cm k (VU)
where c is the capacitance coefﬁcient, k is the conductivity coefﬁcient and U is the
unknown variable: temperature, moisture content, pressure head, etc.
Equation 1-1 applies to transient heat conduction in solids, gas diffusion, ﬂow of ﬂuids
and transport of solutes in porous media, to name a few. The derivation of this equation
is included in almost every engineering and mathematics book dealing with the solution
of the above problems. Powers (1987), Ozisik ( 1980), Patankar (1980), and Churchill
(1987) are a few examples. Applying numerical methods, ﬁnite elements (FE) or ﬁnite
differences (FD) to (1-1), will change the time and space-dependent partial differential
equation (PDE) into a time-dependent system of ordinary differential equations (ODE’s).
This conversion is discussed in many books dealing with numerical solution of PDE;
Segerlind (1984), Smith (1985), and Narasimhan (1978). The ODE has the general form
[dig—l?- +[K]{U}—{F}={0} (14)
where [C] is the capacitance matrix coming from the transient term in the PDE, [K] is
the stiffness matrix coming from the second partial derivative with respect to space in the
PDE, and {F} is the forcing function. Since there is no source term in the PDE this
vector is zero before the boundary conditions are incorporated.
Finite element or ﬁnite difference methods are used to solve (1-2) in the time
domain. Although the FE methods shows clear advantages over FD methods in the space
domain for solving (1-1), this advantage does not extend to the time domain.
There are many schemes available in the literature for solving (1-2) in the time
domain. Many criteria are used to test these schemes including stability, oscillation, and
accuracy. Despite the fact that many authors have presented and discussed solution
procedures for the system of ODE’s in (1-2), there is a lot of art and experience involved
in selecting the proper scheme and the time step that is needed to reach an accurate and
stable solution particularly in two- and three-dimensional problems.
3
APPLICATION TO AGRICULTURAL PROBLEMS
The process of diffusion is well grounded in the biosystem processes. This
section will brieﬂy introduce some of the processes that are governed by the diffusion
convection ﬁeld (1-1). These problems include but are not limited to:
Salt movement and accumulation around micro-irrigation emitters This is important
where the irrigation water contains salt that moves through the soil and eventually
accumulates at the fringe of the wetted region. The problem escalates where the
water used to ﬂush the soil is limited. The same physical phenomena governs the
movement of solutes through the porus media to and through the ground water.
Water infiltration through the soil The fundamental water ﬂow equation in the
unsaturated soil region is governed by a parabolic type equation known as
Richard’s equations.
Grain drying is also governed by diffusion The system here is governed by heat as
well as moisture transfer. Accurate numerical schemes are critical in the optimal
design of grain dryers. Parameters such as the time needed to dry the grain and
the rate of drying are critical in determining the dryer speciﬁcation. The same
can be said about drying manure and other organic products.
Odor transport from animal housing bins The diffusion of odor through air is
parabolic. This is becoming more important as urbanization is expanding in to
existing farm land. The concentration of sulfur in the air as affected by structure-
induced turbulence in the ﬂow is one good example.
4
DOCUMENT LAYOUT
This dissertation consists of discrete chapters. Each chapter is a separate entity
containing introduction, methodology, results, and discussion. Literature citations for
all chapters are at the end of the dissertation. A general literature review (justiﬁcation)
and conclusion for each chapter is included in the separate chapters.
CHAPTER TWO
JUSTIFICATION
There are more names than there are things
Anonymous
A complete literature search of published work on numerical solution of PDE and
their implementation is a huge task. The aim in this section is to show that among the
extensive published resources on the subject, there is a question that the authors have not
answered. This question is ”What is the time step value used in their numerical solution
to ensure stability, eliminate numerical oscillation, and ensure an accurate solution.
Although numerical solution of (1-1) and (1-2) are presented in numerous books
and technical articles, none of these publications explicitly answer the question stated
above. This section discusses some of the literature over the last thirty years. Citations
are concentrated on the ﬁnite element and ﬁnite difference literature.
The solution of the time dependent ﬁeld problems using the ﬁnite element method
was discussed only brieﬂy in early ﬁnite element books. Heubner (1975), discusses the
derivation of the capacitance matrix [C], for transient heat transfer but never discusses
the solution of the resulting system of ODE’s. Zienkiewicz (1971) and Segerlind (1976)
discuss the numerical solution of the system of ODE’s but do not discuss any of the
problems that can arise during the solution process and they do not compare the different
types of elements.
6
More recent books give the time dependent problem greater coverage but could
be misleading to the inexperienced analyst. Allaire (1985), for example, concentrates
most of his discussion on Eulers’s single step explicit method with one-dimensional
problems because the calculations are easily detailed. This method is known to be
unstable (which Allaire discusses) and is not the most accurate of the single step
methods. Allaire dose not discuss any solution in two or three dimensions and makes
no comparison between linear and quadratic elements in the one-dimensional case.
Segerlind (1984) discusses some aspects of the numerical oscillations and physical
reality problems that can arise during the solution procedure. He warns the reader to
avoid using the quadratic elements because of physical reality problems but does not
detail whether the numerical problems are long term or short term and whether the errors
are signiﬁcant.
Until recently, most application oriented books in heat transfer and groundwater
ﬂow centered their discussion of numerical solutions on ﬁnite difference methods. Base
and length of presentations are in favor of the later method. Jaluria and Torrence (1986)
discuss the three node triangular ﬁnite element for solving heat transfer problems but do
not make any comparisons with a two dimensional ﬁnite difference solution. These
authors do not discuss any of the other types of two-dimensional elements. Their
discussion could lead one into thinking that the three node triangle is the most appropriate
for a numerical scheme. Segerlind (1984) indicated that the four node quadrilateral
element is superior to the three node triangle. The presentation in Jaluria and Torrence
(1986) can be contrasted with Patankar (1980) who limits the discussion of the ﬁnite
7
element method to two pages and does not give any equations for the method. Patankar
also recommends the use of the backward difference scheme due to its ”friendliness" for
all values of time and grid size.
Shih (1984) has a chapter on accuracy and error bounds. Most of his discussion
analyzes the error bounds for different orders of the ﬁnite element method. Shih does
not discuss any estimate for At when solving time dependent problems. He also has a
chapter on the comparison of the ﬁnite difference and ﬁnite element methods. He 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”. Shih
(1984), Jaluria and Torrence (1986) and Segerlind (1984) avoid explicit numerical
evaluation of the time step. They discuss stability and numerical oscillation problems but
none give a procedure for estimating the time step as it relates to the accuracy of the
computation. The typical scenario is to present a numerical solution procedure and
compare it with an analytical solution of the PDE using one or more time step values.
The authors, however, never said how they determine what numerical value of time step
to use. Dhat and Touzot (1984) comment that the time step value that eliminates stability
and numerical oscillations may not produce accurate calculations. They also do not give
any suggestions of how to select the time step value.
Gear (1971) and Seer and Bulirsch (1980) discuss the mathematical approaches
to determine a time step value. They deﬁne an error as being the difference between two
8
solutions with time steps of At and At/2 to determine an appropriate step size. This
approach however, does not give much information on how to select a starting value for
At.
Myers (1977) discusses the critical time step for two dimensional heat conduction
transient problems. His discussion, however, centers on estimating the maximum
eigenvalue for use in the Euler stability criterion or the Crank-Nickolson oscillation
criterion of (At 5 2/maximum eigenvalue). Myers does not discuss the determination
of At as it relates to the accuracy of the integration.
Another approach for selecting At is to limit the maximum change in any nodal
value to a certain percent of its previous value. This approach is used in some
commercial ﬁnite element software when solving nonlinear problems. This method
suffers from the need to repeat the calculations if the time step is too large and also does
not give any information on how to select a starting value for At.
Reddy ( 1984) has a section on time dependent problems that is consistent with
much of the mathematics literature. Reddy describes the stability in terms of the roots
of the characteristic equations and the eigenvalues of the global system. Roots of that
equation should be bounded by one to avoid oscillations. For a structural dynamics
problem, he included an estimate for At that gives accurate solution: At = Tum/1', where
Tmill is the smallest period of natural vibration associated with the approximate problem.
According to Reddy, another estimate to At can be obtained from the condition that the
smallest eigenvalue of the characteristic equation be less than one.
Smith (1985) discusses the explicit Euler’s method for solving the non-
9
dimensional form of (1-1). Smith rearranged the difference equation and deﬁned a term
r = cit/(bx)? During the discussion of stability Smith stated that the explicit method is
stable for r values less than 0.5. The implicit Crank-Nickolson 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 and gave time step expression that satisfy both criteria. There was no criteria
for selecting time step based on the accuracy. The term r deﬁned by smith does not
include material properties since cp/k term was deﬁned as one.
Allaire (1985) called Smith’s r term the Courant number. Allaire’s variable
included the material properties. In addition to illustrating stable and non-stable
schemes, Allaire deﬁned an oscillatory stable scheme as having spatial oscillation that
eventually dies out with the solution converging to the correct steady state values.
Allaire showed the following criteria to be true for the single step methods:
0 < r S 0.25 No oscillation
0.25 < r 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 the weighted explicit implicit scheme. For Theta = 0, the scheme
reduces to the explicit method and has stability criterion of r < 0.5. 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) deﬁned Allaire’s (1985) Courant number as the
10
Froude number, Fo. These authors suggested using values of F0 < 0.5 for the implicit
method although lower values gave better accuracy. They never give any example of
what the lower values should be.
One limitation to the use of Courant (Froude) number or the r term deﬁned by
Smith ( 1985) is that the parameter does not change with the boundary conditions. Each
of the authors, Allaire (1985), Jaluria and Torrence (1986) and smith (1985) solved the
heat equation for a different boundary and initial conditions. They recommended
different values of the Courant (Froude) numbers.
Patankar (1991) presented a heat transfer computer program (CONDUCT). This
program uses the backward difference scheme to integrate, but Patankar never discusses
selecting a time step.
Wood (1990) gives an extensive list of time stepping schemes. His list included
most of the known schemes and some new one’s. He studied stability, consistency, and
oscillations where the term "time step" was mentioned at several places. For many of
these schemes numerical results were tabulated using various time steps and the
corresponding error was presented. The author showed that these methods were
consistent with the analytical solution. Wood also refers to the use of time step
adjustment where the time step changes every time step but never give any formula for
determining time step.
The mathematical literature introduces the concept of consistency as a prerequisite
condition to be satisﬁed by any scheme in order to be mathematically sound. Brieﬂy,
this states that as At approaches zero, the error between the actual and simulated
11
solutions should go to zero. What they do not tell us is how small At has to be in order
to satisfy a certain accuracy criteria.
Ortega (1990) deﬁned and discussed three types of errors that are all associated
with the time Step: a) discretization (global) error, b) convergence error, and c)
rounding error. He did not indicate how to deﬁne the numerical value for the time step
that will minimize these errors.
Rushton and Tomlinson (1971) used the alternating direction approach as a
numerical scheme. They studied stability and found that for different boundary
conditions the Courant number, C, that generate accurate time steps changes. For a
sudden change of pressure head on the boundary, C should be less than 1.0. For a draw
down at a well C should be less than 0.05. For sudden change in discharge at a well,
C should be less than 0.5. The authors though, suggest a trial and error procedure for
selecting the optimal At value.
Henrici (1977) had an extensive discussion about the error propagation for the
difference methods in solving the PDE. His theoretical treatment did not include
discussion of the time step for accurate results.
Williams (1980) and Fried (1979) both studied the numerical solutions of PDE
and used the time step criteria that satisﬁed 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 (At 5 2/maximum eigenvalue).
Haghighi and Segerlind (1988) solved the coupled heat and mass transfer
equations using the ﬁnite element method. They used Maadooliat’s (1983) non
12
oscillation criteria as well as the physical reality conditions that Segerlind (1984)
discussed in his book.
Nripendra and Kunze (1991) presented a ﬁnite 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 time step in their paper
Irudayaraj (1991) and Irudayaraj et. a1. (1990) applied the ﬁnite element method
to the solution of a coupled heat and mass transfer problem. Both papers used the
stability criteria for selecting the time step. There was no check whether this time step
ensured accurate results. The authors’ result did not agree with experimental data in the
literature. The stability criteria of was used by Liu et. a1. (1984). These authors used
a modiﬁed Runga Kutta method to solve the parabolic system. Their work did not
discuss solution accuracy.
Peraire et. a1. (1988) studied ﬁnite element solution of ﬂuid ﬂow. They used the
Courant stability criteria of (At 5 K*he/u+c), where c is local sound speed, he is the
average element length, u is the velocity of ﬂuid, and K is a constant.
Alagusundaram et. a1. (1991) applied the ﬁnite element method to model the
diffusion of carbon dioxide in grain bins. Their results were off compared to the
measured values. They listed several reasons for this discrepancy. They did not mention
how they determined the time step. They did not even state what time step value they
use. They did not state whether the size of the time step might be one reason for the
inaccuracy of their calculation.
13
Wood and Lewis (1975) studied seven different ﬁnite difference time marching
schemes. They compared methods based on an accuracy criteria. The authors related
accuracy to oscillations and stability. They determined the critical non oscillatory time
step for the Crank-Nickolson (C-N) based on the maximum eigenvalue. They showed
numerically that when increasing the time step beyond the critical time step oscillations
occurred. Wood and Lewis observed inaccurate values in backward difference scheme
for some time step values. They did not state that accuracy is a separate consideration
in numerical solution of parabolic equations and needed to be addressed and that the time
step needed to be adjusted accordingly.
Cleland and Earle (1984) studied the freezing time of food material using six
ﬁnite 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 called "jumping" of the latent heat peak.
Although there is evidence of accuracy in their solutions, there is no evaluation of a time
step expression that could be translated to other problems.
Abdalla and Singh ( 1985) simulated the thawing of food using the ﬁnite element
method. They presented comparisons between analytical and predicted values. However
they did not state what time step value they used.
Segerlind and Scott (1988) were among the ﬁrst to deal with the time step
estimates from the accuracy perspective. They presented a time step estimate for one-
and two-dimensional problems that produce accurate results. They did not give any
derivation for their estimate and stated that much of it is based on their experience.
14
They did not Show any evidence that their time step estimate really works. However,
they have stated an important observation that the oscillation criteria is too conservative.
The time step requirement is excwded by a factor of two before oscillation were
observed.
Ne-Zheng Sun (1989) studied numerical solutions for the coupled groundwater
ﬂow and advection-dispersion equation. He applied a variation of the linear ﬁnite
element method and compared his solution with analytical ones. No indication was given
as to what time step was used in his analysis.
Scientists reporting new time stepping schemes seem to discuss stability, and
oscillations only. Accuracy is an important criteria in numerical solutions and was not
given enough, if any, attention. Yu and Heinrich (1987), Segal and Praagman (1986),
Fong and Mulkey (1990), Riga] (1990), and Schreyer (1981) used the stability time step
requirement (At < C h2/2), where C is the thermal capacitance when performing a
numerical solution for the heat conduction equation.
Shu-Tung Chu and Hustrulid (1968), and DeBaerdemaeker et. a1. (1977) did not
deﬁne a time step estimate when they discussed numerical solution of the diffusion
equation. Scott (1987) uses the following arbitrary accuracy criteria At=(time to steady
state)/100. In other words Scott assumes that running the problem for 100 time steps
should be sufﬁcient to ensure stability and accuracy. Although this estimate might be a
good starting point for some problems, there is no justiﬁcation for its use nor any
experimental work that back this claim was given.
Maadoliat (1983) studied stability and physical reality oscillations of the ﬁnite
15
element numerical solution. He concluded with a set of conditions that must be satisﬁed
in order to avoid both numerical problems and recommends a time step estimate
accordingly. He did not consider the accuracy criteria.
CHAPTER THREE
OBJECTIVES
Based on a review of the literature, there is a compelling need for establishing a
sound estimate of the time step value that satisfy stability, non oscillations, as well as
accuracy criteria. This should eliminate the trial and error procedure in attaining a
reliable numerical solution.
The most known about the theoretical error estimate is the order of its power
relation which is method dependent. In other words, the error is the product C h”, where
C is a constant, p is the order of accuracy of the integration scheme, and h is the grid
size. Theoretically and for single step methods, p is two for central difference and one
for the other schemes. Practically, these values of p might be far from those theoretical
values and they change during the course of computations (Baker 1993). Without an
analytical expression that quantiﬁes C and its variation with time step this relation is
useless. A complete analytical approach to ﬁnding time step estimate that is practical and
applicable to a variety of problems is a tedious and may be an impossible task. Some
of the approaches discussed in chapter two focusing on evaluating a theoretical error term
could be a challenging and probably an amusing mathematical exercises but are not
functional engineering tools to evaluate a time step. This research uses an experimental
route that may not be at the same level of sophistication as the mathematical approaches
but sufﬁcient to yield applicable and far more useful estimate for a time step that
16
17
accurately integrates the PDE in (1-1). This ability to estimate an accurate time step
becomes extremely important for a non-experienced user that has limited physical
understanding of the problem being solved.
The parameters and/or factors that will affect the numerical solutions of (1-1) are
numerous and all need to be considered. Some of these factors and/or parameters are:
- the dimension of the problem: one, two or three.
- the level of time stepping scheme/s: two, three, or four.
- the time stepping schemes for a particular level.
- the type of capacitance matrix: lumped or consistence.
- the space and time discretization: ﬁnite element or ﬁnite difference.
- type of element: linear, quadratic,
- problem selection including appropriate initial and boundary conditions.
The other important parameters that affect the dynamics of the transient solution
are space and time steps, number of nodes, number of elements, minimum and maximum
eigenvalues, time to steady state, shape of element for two or more dimensions.
The number of options are numerous and need to be narrowed down in order to
manage the problem. Some authors have shown that two level ﬁnite difference time
marching schemes are superior to three or four level time schemes. Wood and Lewis
( 1975) and Shayya et. a1. (1991), Segerlind (1987) and (1992) are examples.
Segerlind (1987) performed a preliminary study on selecting a time step value.
His ﬁndings are the basis for redeﬁning the objectives of this research. Segerlind studied
the accuracy of determining the global system of ODE’S eigenvalues. He studied two
18
types of oscillations; ampliﬁcation factor oscillations, and physical reality oscillations.
Segerlind then studied accuracy and gave a rough estimate for the time step. The main
concepts learned from his study are:
- accuracy of eigenvalues does not imply a solution accuracy. This observation put less
value on the need to accurately estimate the full set of eigenvalues.
- when solving time dependent problems, the linear element is preferred to a quadratic
element.
- lumped capacitance matrix is superior over the consistent matrix formulation.
- ﬁnite difference time discretization is superior to the Galean ﬁnite element time
discretization.
The following restrictions were incorporated in this study:
- use ﬁnite difference two level (single step) time marching schemes: Euler, Crank-
Nickolson, Galerkin, and backward difference.
- use lumped capacitance matrix.
- use linear elements in one-dimensional problems.
- use square bi-linear elements in two-dimensional problems.
The following initial and boundary conditions were considered:
Initial conditions: sine wave, linear variation, and uniform value with a step change at
each end at time zero.
Boundary conditions: Dirichlet type and derivative boundary condition.
19
The overall objective of this research was to develop a time step estimate for
solving (1-2) that satisfies an accuracy criteria, satisfies stability requirements and
to recommend a numerical scheme to be used for this equation.
The specific objectives for the study are:
1. Develop a time-step estimate for each numerical scheme that satisﬁes an accuracy
criteria for one-dimensional linear elements and two-dimesional square elements.
2. Develop recommendations on the solution of the parabolic partial differential equation
based on the ﬁrst objective.
3. Compare and analyze the time step estimate of objective one with a commonly used
time step criteria such as the Courant or Froude numbers.
4. Compare the accuracy of the ﬁnite element and ﬁnite difference space formulations
for two-dimensional problems.
5- Study the effect of element size and shape on the maximum eigenvalue for two-
dimensional four nodes quadrilateral elements. Comparisons with the three node
triangular elements and recommend which elements to use in transient two-
dimensional ﬁeld problems.
CHAPTER FOUR
ONE-DIMENSIONAL FIELD PROBLEMS
Parabolic differential equations governs a variety of time-dependent problems in
science and engineering. One-dimensional ﬁeld problems have the following general
form
122 =1) 1‘]. (4.1)
" 6X2 ‘ at
Equation (4-1) applies to a transient heat conduction in solids, gas diffusion in air,
ﬂow of ﬂuids and transport of solutes in porous media, to name a few. U is the
temperature in a heat conduction problem, pressure in a ﬂow problem, or concentration
in a solute transport problem. The material properties D, and Dt could be lumped into
one variable (D, / D.) commonly referred to as the thermal diffusivity in heat transfer
problems, 0:.
Applying numerical procedures such as the ﬁnite element (FE) or the ﬁnite
difference methods (FD) to (4-1) with its boundary and initial conditions changes the time
and space-dependent partial differential equation (PDE) into a time-dependent system of
ordinary differential equation (ODE’s) as in Equation (4-2).
where {U} is the vector set of the unknown nodal values, [K] is the stiffness matrix, and
20
21
mag” + [K].{U} - {F} = {0} (4,2)
comes from the second partial of U with respect to x, [C] is the capacitance matrix, and
comes from the time dependent term in the PDE, while {F} is the force vector, comes
from the source term in the PDE as well as the boundary conditions. All matrices are
global and built from element matrices using the standard stiffness procedure, Segerlind
( 1984).
The global matrices [C] and [K] are positive deﬁnite symmetric matrices. [K] is
singular before boundary conditions are imposed. The matrices [C] and [K] are built
from element contributions whose coefﬁcients depend on the type of interpolation
function used to solve the problem. Two types of capacitance matrices, [c‘°’], and the
stiffness matrix, [k‘°’], for the linear one-dimensional element are summarized in (4-3).
The lumped capacitance matrix comes from using a weighing coefﬁcient in the Galean
derivation of [c‘°’] different from the one used to develop [k‘°’]. The same weighing
functions are used to obtain both [c‘°’] and [k‘"] in the consistent formulation. Although
many investigators indicated that both capacitance matrices have comparable accuracy,
there are some disadvantages associated with the consistent formulations pertaining to the
stability of the numerical solution, Visser (1965), Wilson and Nickell (1966), Brocci
(1969), and Zienkiewicz (1977). Other capacitance matrices have been proposed, the
average and the optimum, Segerlind (1987).
22
Linear Element Matricies
Capacitance
Consistent - [cm] _ D, h 2 1
' 6 1 2
_D' h 1 0 (4-3)
Lumped : [C (e)] =
O 1
Stiﬁ‘ness
D 1 -1
(e) = 1
[k 1 T [-1 1]
Analytical Solution of the System of Ordinary Differential Equations
The system in (4-2) can be solved analytically using modal analysis, Cheney and
Kincaid (1985) and Potter and Goldberg (1987), a summary of which is presented here.
The method was programmed and used to compare the analytical solution of PDE and
the numerical solution of the system of ODE’s.
Multiplying (4-2) by [C]'1 gives:
3f}? + [K']{U} - {F'} = {0} (*4)
where:
[K'] = [C]'1 [K]
{F'} = [C]" {F}
The solution of (4-4) is obtained by assuming a solution of the type:
{Um} = {a} e'“ (*5)
and substituting it into (4-4). The result is the standard eigenproblem
23
([K'] - A [11) {U} = {0} “'6’
where A is the eigenvalue and {U} the associated eigenvector. Since [C] and [K] are
symmetric and positive deﬁnite all the eigenvalues are real, positive and distinct and the
associated eigenvectors are also real and linearly independent.
If we deﬁne [E] as a matrix whose columns are the eigenvectors of
([K‘]-A[I]){U} ={0}, then [B]1 [K’][E] = [A] , a diagonal matrix whose values are the
eigenvalues of [K‘], Strang (1976). The [E] matrix is used to uncouple the system of
differential equations.
Deﬁne a new vector of nodal values {Z} as {U} = [E]{Z}, then 3{U}/6t = [E]
6{Z}/6t. Substitution into (4-4) and multiplying by [E]1 yields
{2'} + [A] {z} - {F‘}={0} (*7)
Since [A] is a diagonal matrix, the system of linear differential equations has been
uncoupled in the new variable {Z}. An individual equation in (4-7) has the form
Zr+xin-fr=0 (4.8)
with a solution of
z,(:) = b; ‘V + {:1 (4-9)
1
where bi is determined using the transformed initial conditions, {2(0)} = [E]" {U(0)}.
The ﬁnal solution is a linear combination of the Z,(t) solutions
24
(1,0) = Ele(t) + 15,2220) +... + Eﬁnzno) (4—10)
the modal analysis approach is seldom used because of its extensive computations. It is
used in this study to evaluate the analytical solution of the system of ODE’S.
Numerical Solution of the System of Ordinary Differential Equation
Numerically, ﬁnite element or ﬁnite difference methods can be used to solve (4-2)
in the time domain. Although the ﬁnite element method shows clear advantages over the
ﬁnite difference method in the space domain this advantage does not extend to the time
domain, Segerlind (1984). There are many time domain schemes available in the
literature for each method. Several criteria are used to test these schemes including
stability, oscillation, and accuracy.
The numerical solution for the ODE’s is usually conﬁned to the single step
schemes deﬁned by
([Cl+04t[K]){U}2 = (IO-<1—0)Ath1){U}.+At(0{d.+<1-0){F}.) (4'11)
where
0 = O is the Euler scheme
0 = 0.5 is the central difference scheme
0 = 0.67 is the Galean scheme
0 = 1.0 is the backward difference scheme
0 = 0, 0.5, and 1 are associated with the ﬁnite difference method, while 0 = 0.6667 is
associated with the ﬁnite element method, Segerlind (1984).
25
The central difference (CD) or the Crank-Nickolson scheme is second order
accurate in time compared to the other single step schemes that are only ﬁrst order
accurate in time, Gear ( 1971). The central difference scheme has been shown by many
authors to be an accurate and stable scheme for solving (4-12).
The Objectives for the One-Dimensional Study are:
1. Develop a time-step estimate for each numerical scheme, of (4-1 1), that satisﬁes an
accuracy criteria.
2. Develop recommendations on the solution of the parabolic partial differential
equation based on the ﬁrst objective.
3. Compare and analyze the time step estimate of objective 1 with a commonly used
time step criteria such as Courant or Froude numbers.
METHODOLOGY
The main criteria for the numerical schemes used in this research was accuracy.
Stability and oscillations were not a direct primary part of the study. The accuracy
criteria must satisfy any stability requirement for the numerical scheme under
consideration. This section will present the methodology for developing a time step
estimate equations that will integrate the ODE accurately.
The analytical and numerical solutions of the system of ODE’s were compared
to the analytical solution of the PDE for the unit step change problem. The step change
problem was selected because its analytical solution contains all the frequency
26
components and has the shortest time to steady state.
The initial conditions for the step change problem are
u(x,0) = 1 0 S. x s 1 (4-12)
with the end points satisfying the boundary conditions
u(0,t) = u (l,t) =0 t > 0 (4-13)
These conditions are shown in Figure 4-1.
The analytical solution for the step change problem is
on
u(x,t) = [42 1 (sin(2n+1)wx)e“(2’“”"")]
1.0 (MI)
(4.14)
The comparisons between the numerical and analytical solutions were made using
an accuracy error ratio, e, deﬁned by
)3 2 INoDE, -A PDEU|
e = j-l M (4.15)
2 )3 IAODEg-APDEyI
j-l i-l
where NODE is the numerical solution for the system of ODE’s, APDE is the analytical
solution for the PDE, AODE is the analytical solution for the system of ODE’S, n is the
sampling point in the space domain, m is the sampling point in the time domain. The
accuracy ratio is based on the assumption that the numerical solution of the system of
ODE’s must be as good as analytical solution of the ODE’s.
27
U(x,0) = 1, 0 <= x <= 1
U(0,t) = U(1,t) = 0, t > 0
Initial Values 1
Figure 4-1. Graphical representation of boundary and initial condition for the step
change problem.
28
The accuracy ratio in (4-15) is a ratio of two L1 norms where each norm is
evaluated for a set of points in the space-time domain. Other deﬁnitions of the error
ratio were studied including the L2 and the Lou norm. These ratios resulted in a less
accurate estimation of the time step through out the process of deﬁning a time step
estimate. The inﬁnity norm was less conservative than the L1 norm in some cases, but,
underestimated the time step in other cases.
The sampling point locations in the space dimension were at each node. The
sampling point locations in the time dimension were at or very near 2.5, 5, 10, 20, 40,
and 80 percent of the time to steady state (t,,) where t,, was deﬁned as
t = (416)
:3.
rs A,
where Al is the lowest eigenvalue for the system. This eigenvalue was evaluated
numerically (J ACOBI method) or analytically when an analytical solution was available.
The deﬁnition in (4-16) comes from e’“, the ﬁrst term in the analytical solution
which lasts the longest as t increases; when Alt = 4, e4 = 0.018, and over 98% of
the transient has been completed.
An existing one dimensional ﬁnite element program was modiﬁed for the use in
this investigation, Segerlind (1987). The computer program evaluated the analytical and
numerical solutions for the PDE and the ODE’S as well as the eigenvalues for the
eigensystem of (4-6). The program evaluates four different capacitance matrices
including the lumped and the distributed (consistent) matrices. The lumped formulation
and a linear element was used in the current investigation. The analytical solutions of
29
the PDE (APDE) as well as the numerical solution of the ODE (NODE) were checked
against similar problems found in the literature Smith (1985) as well as hand calculations.
All tests indicated that the computer programs used in this investigation were
programmed correctly.
Time Step Estimates
The speciﬁcation of the time step to overcome oscillations and physical reality
problems has been addressed by many authors, Segerlind and Scott ( 1988), Ortega
(1990), and Ortiz and Nour-Omid (1986). The time step speciﬁcation is complicated by
its dependency on the grid, boundary and initial conditions. It also changes with the time
stepping scheme used. In order to estimate the time step nwded to numerically integrate
the PDE accurately, a numerical experiment was conducted. In this experiment the step
change problem with the linear lumped capacitance matrix and D, = D, = 1, were used.
Eight different uniform grids, each with a different number of nodes (N) in the region
[0,1] were analyzed. For each N, the time step was allowed to vary generating a series
of numerical scenarios. Each scenario was characterized by time step, N, initial and
boundary conditions, material properties, and the numerical scheme. The L1 norm for
the difference between the APDE and NODE for all nodes at each sampling point was
computed. A similar L1 norm for the difference between APDE and AODE for all
nodes at each sampling point was also computed. Each L1 norm was the sum of 6 by
N values of the difference between APDE and either AODE or NODE. The value of
six comes from the fact that there are six sampling points in the time domain, and N is
30
the number of sampling points in the space domain. The error ratio deﬁned earlier in
(4-15) and used in this investigation is the ratio of these Ll norms. Each case, deﬁned
earlier in this section, produced one error ratio value. Varying N and time step for a
certain case creates a set of points that ﬁts a parabolic distribution curve. A sample
curve for the CD scheme is shown in Figure 4-2. Similar graphs were generated for the
other three single step schemes of (4-11). A common feature among these results is that
as the size of the time step was increased for a particular N, the error ratio increased.
In other words, the accuracy of any scenario improves as the time step gets smaller.
Nevertheless, there is a time step below which the accuracy improvement for any
scenario does not justify the extra computation generated by reducing the time step. The
time step at which the computation is considered at an acceptable accuracy was set at 5 %
error level. This means that the NODE solution is allowed to deviate 5% from the
AODE. The time step corresponding to this accuracy point was found from direct
numerical experimentation values and could be recognized from a blown up version of
Figure 4-2. This ﬁgure shows that as N increases the time step should decrease to stay
accurate. It also shows that the relation between time step and accuracy for a certain
mesh is a non-linear power function.
The time step value that gives numerical results that are as good as the AODE
within the 5 % error tolerance was recorded for each case. The data sets corresponding
to N and the time step were found to be ﬁtting a power function of the form At=aN",
where a and b are parameters to be determined.
31
Changes in Error Ratio As A function of Time Step And
Number Of Nodes For Theta = 0.5
4 v e
.2 3 ..
:‘2
h 2 ..
o /‘
“E l .. Om—o‘Q/e/O “A 's‘ -<.‘ ’43 I
0 3 3 3 f c 4:
0 0.005 0.01 0.015 0.02 0.025 0.03
Time Step in seconds
O
"""'—N=5 +N=6 N=7 —°_N=9
—‘—N=ll +N=13 N=15—0—N=17
Figure 4-2. Variation in the error ratio as a function of the time step and the number
of nodes for the central difference scheme.
32
The power relation can be linearized by taking the natural logarithm of both sides. The
variables N and At were then regressed to ﬁnd a best estimate for the time step as a
function of N. The linear regression equation was then transformed into its non linear
power form by taking the exponential of both sides of the regression equation.
Since the time to steady state is related to the smallest eigenvalue, this eigenvalue
was incorporated as part of the regression equation. This means that instead of
regressing the time step against N, the product of At and A, was regressed against N.
This allows the time step estimate to vary with the boundary conditions as well as the
material properties. Including such parameters into the At estimates makes this estimate
dynamic.
Figure 4-3 shows the set of data points, corresponding to the time steps that gave
numerical results as accurate as the AODE, along with the time step prediction regression
equation (4-18) for the central difference scheme. It shows that for N > 7 the relation
between N and time step can be considered linear. There are two regions separated by
the N = 7 that can be regressed linearly and independently. This means separate
equations for each region. Two independent time step prediction equations for each
scheme depending on the mesh is not appealing from the application point of view. In
addition, since the objective was to obtain a single estimate two separate equations was
not desirable. A similar regression equation was developed for each scheme. The
division at N=7 was noticed for all schemes. Furthermore, the AODE solution was not
accurate for N < 7. Therefore, the time step estimate was developed only for N > 7.
33
Time Step That Give Numerical Results As Good As
ANODE Theta = 0.5
(125 -
T
0.2 ~- '
OJ -'
Time Step ‘ Lambda 1
(105 "
0 2 4 6 8 10 l 2 l 4 16 l 8
Number Of Nodes
—'— NUMERICAL —D— REGRESION
Figure 4-3. Actual and regression curves for the data points that gave accurate numerical
results as good as AODE for the central difference scheme.
34
Figure 4-4 shows a plot of the set of data points corresponding to the time steps
that gave numerical results as accurate as the AODE for the four schemes. The central
difference scheme was able to accommodate a larger time step for a particular problem
and still maintain high accuracy. This is not true for N < 7 where the Galerkin schemes
is superior. Due to scaling factors the other three schemes appear to be lumped together.
A closer look of the ﬁgure shows that the Galerkin scheme is next best in allowing larger
time steps for a certain problems and while maintaining a certain accuracy. The other
two schemes, backward difference and Euler show comparable behaviors.
Figure 4-5 shows a natural logarithmic scale of the data in Figure 4-4. It was
included to show the parallel trends in the 0 > 0.5 and 0 < 0.5 schemes. Euler and
central difference scheme have similar powers that makes them run parallel lines on a
log-log scale. They have different coefﬁcients. Similarly, the backward difference and
Galerkin schemes run parallel but with different coefﬁcients. This behavior may or may
not have any signiﬁcance. 1
The time step prediction equations that were developed as described above for
the step change problem using the four single step integration schemes are shown in (4-
17) through (440).
35
Time Step Values For Numerical Results That Are As
Accurate AS Analytical ODE
0.025 ~-
0.02 «»
33‘ 0015 ~-
U:
E 0.01 .-
E—
0.005 1- \
I .\
0 ¢ : W , 3 , L w 3.___,
0 2 4 6 8 10 12 14 16 18
Number of Nodes
——-— theta=0 —D— theta=.5 . theta=.6667 —<>—theta=l
Figure 4-4. Actual curves for the data points that gave accurate numerical results as
good as AODE all schemes.
36
Time Step Values For Numerical Results That Are As Accurate As
Analytical ODE
Ln—Ln Scale
0 : 1
2; -1 CD 0.5 l 3
E -2 ~-
5 -3 ..
«- -4 J.
D.-
.2 -5 ..
ﬁ -6 ,-
E -7 -_
-8 ..
Number of Eigenvalues
—"—‘"‘ theta=0 —G— theta=.5 ’ theta=.6667 ——°— theta=1
Figure 4-5. A logarithmic plot representation of the data shown in'Figure 4-4.
37
Euler Scheme _ .27N-1.6
At - x, (447)
Central Difference _ 1.13N'1-‘8
A’ ‘ ———>\l (4-18)
Galerkin Scheme _ 70N '3-79
A’ ‘ A, (4.19)
Backward Difference _ 30.6N'3-9‘
A’ ‘ —>\1 (4.20)
These prediction equations are the basis for further investigations and satisfy the
main objective in this study.
Unlike the Froude, Courant, and r values reviewed earlier, the above equations
are dynamic. The time step estimate changes with the boundary conditions, and material
properties. These changes are captured by the smallest eigenvalue. In addition, the time
step estimate varies with the mesh size.
The above time step estimates, (4-17) through (4-20), have been developed based
on a lumped error deﬁned in (4-15). These estimate deﬁne a time step that will give
overall accurate results and not at any particular time. It was shown that error
distribution is not linear through out the time domain. The error is highest at small time
values and dies out as time increases. The error distribution also changes with the
integration scheme. Some schemes are more accurate than others at the same point in
the time domain. The fact that the error is highest at small time values might explain
why in the central difference scheme increasing the time step did not reduce the
38
accuracy. In the later scheme, the inaccuracies present at small time values will not be
included for large values of the time step. On the other hand smaller time steps capture
the numerical error present at small time values.
EVALUATION OF THE TINIE STEP ESTIMATES
After developing the time step estimates using the unit step change problem, the
equations were used to estimate the time step needed to solve other problems with
different initial and boundary conditions. The main objective of this section is to show
that the time step estimate equations, (4-17) through (4-20) developed in the previous
section, are applicable for a different set of problems within the range of meshing used
to develop them, 5 < N < 17, in the domain [0,1]. These equations were developed
using the most rigorous problem, a step change. Applying the time step estimate to other
problems should provide a conservative estimate.
The four evaluation problems were
a) Sine Wave
u(x,0) = Sin 1rx, 0 S x s l
u(0,t) = u(1,t) = 0, t > 0
where the analytical solution, u, consists of a single term, is given by the following
u(x,t) = e("1l'10(81n1rX) (4'21)
Even though the complete set of eigenvalues exist for this problem, An = nzrrz, all of the
39
coefﬁcients of the analytical solution except the ﬁrst are zero. It is important to observe
that the eigenvalues, A“, expression is shared by problem b below, as well as the step
change problem in the previous section. This is because all of these problems has the
same zero boundary conditions.
b) Ling Vao'ag'on
u(x,0) = 2x, 0 s x S 1/2
u(x,0) = 2(1-x), 1/2 5 x S 1
u(0,t) = u(l,t) = 0, t > 0
where the analytical solution, u, is given by the following
u(x,t) = f: i2(sinﬂ)(sinmrx)e“""”’ 14-22)
.3.
«Ina. n 2
Note that since sin nw/2 = 0 for n = even, only odd numbered eigenvalues enter
into the solution. The discontinuity at x = 1/2 produces the largest error at that point, this
error will decrease with time as long as the boundary conditions remain constant, Smith
( 1985). The coefﬁcient l/n2 goes to zero faster than l/(2n+ 1) in the step change
problem, thus faster convergence occurs with the analytical solution.
c) Envativo Boundary Conditions 1
u(x,0)=1, 05x51
au(0,t)/ax = u, t > 0
6u(l,t)/6x = - u, t > 0
40
This is equivalent to
6u(0,t)/8x
Mu - S, t > 0
6u(l,t)/8x = - Mu - S, t > 0
where M = l, S = 0 are the parameters that are used in the ﬁnite element formulation,
Segerlind (1984).
The analytical solution is
u(x,t) = 42 secan
cos2an(x -0.5)e‘"‘“"" (4-23)
”-1 (3 +403)
where the eigenvalues 011,, are the solutions of a tan a = 0.5, Smith (1985).
d) W
u(x,0)=1, 05x51
6u(0,t)/6x = 0, t > 0
6u(l,t)/ax = -H2 u, t > 0
This is equivalent to
6u(0,t)/6x = 0, t > 0
au(l,t)/8x = - Mu - S, t > 0
where M = H2, S = 0
The analytical solution is
°° e ’“ﬂ-I‘H, 00ng
= 2 (4-24
u(x’t) 2u°,,., “139322sz cosB L )
where Bm’s are positive roots of 6,, tan BmL = H2, Ozisik (1980).
41
The ﬁrst two problems a) and b), have the same A, as the step change problem,
of 9.87 and the time to steady state of 0.4053 sec. The same sampling points time
domain were used for these three problems and were 0.01, 0.02, 0.04, 0.08, 0.16, and
0.32 seconds. Problem c) has a A, = 1.707 and a time to steady state of 2.343 sec. The
time domain sampling points were .06, .12, .24, .48, .96, and 1.92 sec. Problem d) has
A, = 2.04 and the time to steady state is 0.98 sec. The time sampling points were .025,
.05, .1, .2, .4, and .8 sec. Time sampling points for all of the four problems were at
or near 2.5, 5, 10, 20, 40, and 80 percent of time to steady state.
The end goal in any numerical computation is to achieve accurate results that are
as close to the analytical PDE as possible. For that reason a different accuracy ratio was
deﬁned as:
)3 NODE,
e = L.— (4-25)
2 APDE,
M
where the sum in both the numerator and denominator is over the number of nodes, N
and the sampling points in the time domain. The time step was evaluated for the ﬁve
problems deﬁned earlier, including the step change, for two mesh sizes, using the four
single step scheme time step estimate equations (4-17) through (4-20). Multiples of this
time step, one-half, twice, three times, were computed to deﬁne a series of numerical
runs whose list constitute the complete combination of the following parameters:
- numerical scheme where 0 takes values of 0, 0.5, .6667, or 1.
- Number of nodes, N, was: 8 or 16 in the interval [0,1].
42
- Problem initial and boundary conditions were as follows: sine wave, linear variation,
and step change initial conditions. Two different derivative boundary conditions
where M = 1 = H2, S=0 for both conditions.
- Time step At: estimated from (4-17) through (4-20) and its multiples.
For each of the above scenarios, the accuracy ratio ,e ,in (4-25) was computed
at each of the six time domain sampling points. The ratio, e, was less or greater than
one and any deviation from one was considered a deviation from accuracy. Values of
e > 1 means that the sum of the numerical nodal values of U at any time domain
sampling point is greater than the APDE. Conversely, e < 1 means that the sum of the
numerical nodal values of U at any time domain sampling point is less than the APDE.
For that reason the sum of the absolute deviation from one for all the six sampling points
in the time domain was computed. This sum was deﬁned as the sum of the deviation.
The average of this sum across the six sampling points was also computed. This average
deviation is in a sense the error in the computation for that problem. Thus each scenario
produced a single number representing the numerical error or deviation from the APDE
solution. Starting from our development of the time step estimates, this error should be
less than 5% for the computed time step. The error should increase beyond 5% for
some multiples of that computed time step. The half computed time step value solution
is expected to have an accuracy comparable to the solution at the computed time step.
The number of calculated time steps is a new parameter that was introduced to study the
performance of the estimated time step. It is a non-dimensional parameter and
represented the actual time step used in the computation scaled by the estimated time step
43
for that problem. A value of two, for example, represents a different time value for
different problems. This parameter was only used as an evaluation criteria for the
estimated time step for each combination. It is used to prove that the criteria works for
various problems and not for a comparison between problems.
Figure 4-6 shows the performance of the central difference scheme for N= 16 and
three different initial conditions. Figure 4-7 shows the calculated time step performance
for the second derivative boundary condition problem for the four different single step
schemes. Similar curves were generated for other scenarios. The x-axis in these Figures
was labeled as the number of calculated time step which is the ratio of the time step and
the estimated time step from (4-17) - (4-20). The y-axis is the average deviation which
is the decimal error of the NODE compared to the APDE. Both ﬁgures Show that as the
number of calculated time steps increased beyond a value ranging from 1.0 to 4.0 the
calculations become inaccurate.
Figure 4-6 shows that the time step estimate is accurate for the step change
problem and conservative for the sine wave and the linear variation problem. As the
time step used in numerical computation exceeds twice the estimated time step, solution
for step change problems become less accurate. The other initial condition problems
accommodate larger time steps before they show inaccuracy. For all problems the ﬁgure
shows that reducing the time step below the calculated time step did not improve
accuracy.
44
Calculated Time Step Performance Theta = 0.5, N = 16
0.09 ‘-
0.08 1' 0
0.07 '-
0.00 0
0.05 1-
0.04 0
0.03 “r
0.02 "
0.01
1011
vrat
Average De
0 l 2 3 4 5 6 7 8 9 10
Number of Calculated Time Steps
O
—'— Sine Wave —0— Lin Variation Step Change
Figure 4-6. Performance of the estimated time step for the CD scheme and N =16
grid problem for three different problems.
45
Calculated Time Step Performance for The Derivative Boundary
Condition Problem 11:10, S=O at X=1 , N=16
7.1 0.03 ~-
g 0.025 "
Q)
a
a 0.02 -»
.O’IT:
e ° 0.015 «~
9 0.01 v
Q)
ED
1‘3 0.005 "
2
< 0 < : t c : : : :
0 2 4 6 8 10 12 14
Number of Calculated Time Step
—'— theta=0 —D— theta=0.5
O
theta=.667 —°— theta=l
Figure 4-7. Performance of the estimated time step for the derivative boundary condition
2 problem for the various schemes and N = 16 grid.
46
Figure 4-7 shows that for 0=0.6667 and 1, time steps larger than the estimated
value produce inaccuracies. For 0=0 time steps larger than the calculated were unstable.
For the central difference scheme it took a time step of four times the estimated value
to produce inaccuracies. Time steps below the calculated value did not improve the
accuracy for the central difference scheme.
In conclusion, (4-17) through (4-20) produce time step values that produce
accurate numerical results when compared with analytical solutions. Although the time
step estimate was somewhat conservative for a few problems, it gave good results for all
the problems analyzed. It is of great importance to note that the time step estimate not
only ensures accuracy but was also below the stability and oscillation requirements which
conﬁrm the belief that If oscillations occur in the calculated values, the time step is too
large for accurate results.
Effect of Element Size on Accuracy
The mesh size 6x and the time step 6t are two interdependent variables. What has
been studied so far is to how to estimate the time step that produces accurate results for
a certain mesh. The mesh size was ﬁxed and the time step was optimized. Another
approach to accuracy is to select the proper mesh that produce accurate results for a ﬁxed
time step. This approach may not be feasible, but, the impact of the grid size on the
numerical accuracy was studied hoping to increase the understanding of this variable.
The effect of element size on the solution accuracy was investigated by solving
the step change problem deﬁned earlier with a constant time step of 0.005 seconds while
47
changing the number of elements on the [0,1] interval from 3 to 21. The linear lumped
element was used with D, = D, = 1.
The same error ratio criteria deﬁned in (4-15) was used in this investigation. For
each scheme the L1 norm was computed, as the sum of the deviation from APDE as well
as AODE (at 5t=0.005 for various N values in [0,1]), for all nodes and each sampling
point. The error ratio was then determined as the ratio of the two norms. Note that
AODE solution is independent of the time step, therefore it needed to be evaluated only
once for each grid size.
The error ratio as related to the number of elements for the four schemes is
shown in Figure 4-8. The results show that all schemes become less accurate relative
to the AODE as the element size decreases. The ﬁgure also shows that the central
difference scheme is more tolerant than the rest of the schemes to decreasing element size
while maintaining a constant (St.
The Galerkin scheme was second best followed by the forward and backward
difference schemes which have about the same accuracy. The important conclusion here
is that to maintain a high level of accuracy, the time step should decrease as grid is
reﬁned.
48
Changes In Error Ratio As A Function 0f Element Size For The Four
Methods; Time Step Was Fixed To 0.005; Step Change Problem
5 ' e
4.5 _
4 ..
~2- 3.5 ~- -
“E 3 -~ /.
,5 2.5 -~
2 0 /
1.5 ~- .6 /
I ' 0. /ﬁ _.-;———— w : t
0 5 10 15 2O 25
Number Of Nodes
—'— Theta = 0 —D— Theta = 0.5 ’ Theta = 0.6667 —°— Theta = l
Figure 4-8. Effect of element size on accuracy represented by the variation in error ratio
as a function of number of nodes.
49
Comparisons With Froude or Courant Number
The time step estimates of (4-17) through (4-20) are compared with the Froude
or Courant Numbers. Both numbers are deﬁned as
Froude N0. = kAt 1426)
pch2
where k, c, and p are given parameters in the PDE. The three variables are often taken
as unity for simplicity.
The Froude number is a function of grid size, time step, and material property.
It does not change with the boundary condition. Therefore the time step criteria
developed using the Froude number does not change if the boundary conditions of the
problem are changed. The boundary conditions have an impact on the dynamics of the
problem. The speed at which the problem moves to steady state is a function of the
boundary conditions. The time to steady state is a function of the eigenvalues that
depend on the boundary conditions.
The comparison of the time step estimate (4-17) through (4-20) and the Froude
number has to be based on a ﬁxed value of A,. Using a cartesian coordinate with N
being the x-axis, and the time step as the y-axis, the Froude number equation plots as a
series of curves similar in Shape to the time step (4-17) through (4-20). A regression
ﬁt was performed for each scheme and the Optimal Froude number was determined for
a ﬁxed value of A,. This optimal Froude number reduces the difference between the time
step estimate using (4-26) and the time step estimates using (4-17) through (4-20) to zero.
The corresponding Froude number was deﬁned as the Froude number equivalence.
50
Results of this regression for A, =9.8 are shown in Figure 4-9 for the central difference
scheme. Similar results were produced for the other three schemes and various A,
values. Each scheme yielded different and distinct Froude number equivalences which
are given in Table 1.
The Froude number equivalence of the backward difference and Euler schemes
were the same. This suggests that the two schemes have comparable accuracy. This
result is not a surprise from the theoretical point of View since both schemes stand at
equal distances and opposite sides from the center of the approximation domain. The
values also conﬁrm an earlier observation that the Euler and backward difference
schemes run with parallel accuracy and there is no real advantage in using the backward
difference scheme. It is inaccurate in the region where Euler scheme is unstable.
Except for the central difference scheme, the optimal Froude numbers are quite
different from the values given in the literature.
Figure 4—10 shows the Froude number equivalence as a function of 0 (scheme) for
three different values of A,. The distribution for each curve is symmetric about 0=0.5
(central difference). Although there is an inﬁnite series of such curves the ﬁgure Shows
important features and characteristic of the comparison between the two time step
estimate approaches . The ﬁgure shows that the central difference scheme is superior
to the other schemes. The difference between central difference scheme and the other
schemes is signiﬁcant. The other important feature is that the Froude number
equivalence is not a ﬁxed parameter and is not sufﬁcient to describe the dynamic
behavior of the time dependent solution of the parabolic equation.
51
Predicted Time Step Compared to Froude Number
0.016 .. I
0.014 r-
0.012 *-
0.01 ~r
0.008 -
' 0.006 ,-
0.004 -
0.002 *-
Tlme Step
; I“¢!=I=I_ :—
1,
0 2 4 6 8 10 12 14 16 18 20
Number Of Nodes
—'— Theta=0.6667 “—0—- Fo/N~2
Figure 4-9. Comparison between estimated time step and Froude Number for the
Galerkin scheme.
52
Froude Numer Equilelance Changes With the Smallest
Eigenvalue
F roude Number
O
_'— lambda 1:9.8 —0— Lambda 1 =2.707 Lambda 1 = 1.7
Figure 4-10. Froude number equivalence changes with A.
53
In addition, The difference between Froude number equivalence for each A, value is
signiﬁcant.
The Single Step Scheme 0 = 0.33
Given a function U(t) and the interval [a,b], the mean value theorem for
differentiation can be used to develop an equation for U(t). This theorem states that
there is a value of t, call it C, such that ¢(b)-4>(a)=(b-a) d6/dt(¢), where c=a
corresponds to the Euler’s scheme, c=b corresponds to backward difference scheme.
In both cases the mean value theorem, the basic principle behind all single step schemes,
approximates the value of U(t) at the center of the domain [a,b] or at At/2 by its value
at either a (for Euler) or at b (for backward difference). Since the function is symmetric,
the error associated by either schemes is the same. Although (4-17) and (4-20)
corresponding to 0=0 and 1 respectively give different time step values, the calculated
values are not signiﬁcantly different. The mean value theorem is explained in Cheney
and Kincaid (1985), Potter and Goldberg (1987), and Kreyszic (1962).
From the discussion of the mean value theorem, if we deﬁne another single step
scheme with 0:.333, this scheme should have accuracy comparable to the Galerkin
scheme. To test our hypothesis a numerical experiment was conducted. The experiment
uses the step change problem deﬁned in the methodology section, two different mesh
sizes N =8, and 16, 0:.333, and time step estimates computed by the Galerkin scheme,
(4-19).
54
The objective is to show that the same time step estimate could be used for both
schemes 0=.6667 and 0.333. The same procedure of evaluating the average deviation
was followed as in the evaluation section above. The results are seen in Figure 4-11.
This ﬁgure shows that the time step estimate of 0:0.6667, (4-19) can be used when
0=0.333. For the case of N = 16, the time step estimate was somehow conservative, it
become inaccurate only after 10 multiples of the estimated time step was exceeded. For
N=8, time step values greater than twice the estimated time step yield inaccurate results.
Similar to all the cases analyzed, the ﬁner mesh (N = 16) resulted in improved accuracy
(compared to N=8).
The optimal Froude number equivalence values for each scheme (0: 0, 0.333,
0.5, 0.6667, and 1) and for three different values of A, are listed in Table 4-1.
Table 4-1. Froude number equivalence for various A, values:
THETA x, = 9.8 x, = 2.707 x, = 1.7
0 0.05 0.2 0.3
0.333 0.2 0.5 0.7
0.5 0.5 2.0 2.7
0.6667 0.2 0.5 0.7
1 0.05 0.2 0.3
55
Calculated Time Step Performance Theta = 0.3
for the Step Change Problem
Avera e Devratron
Error)
.0 .0
c3 9) C)
o o —-
0n _. or
\ r. :
0 5 10 15 20 25 30
Number of Calculated Time Step
—'—N=8 —0—-N=16
Figure 4-11. Calculated time step performance for the step change problem and
0:0.333
CHAPTER FIVE
TWO-DIMENSIONAL SQUARE GRIDS
Parabolic differential equations governs a variety of time-dependent problems in
science and engineering. The two-dimensional ﬁeld equation have the following general
form:
2 2
Dﬂ+Dﬂ+GU+Q=DaU (5-1)
" 3X2 ’ 3Y2 ‘ E
where the second partial of U with respect to X or Y are the diffusive term, Q is the sink
or source term and the partial of U with respect to time is the transient term.
Equation 5-1 can be applied to a transient heat conduction in solids, gas diffusion
in air, ﬂow of ﬂuids and transport of solutes in porous media, to name a few. U is the
temperature in a heat conduction problem, pressure in a ﬂow problem, or concentration
in a solute transport problem. For isotropic media, the material properties D,, Dy and
D, could be lumped into one variable (D/ D,) commonly referred to the thermal
diffusivity in thermal problems, a, where D,=D,=D.
Integrating (5-1) with its boundary and initial conditions over space coordinates
using the ﬁnite element (FE) or ﬁnite difference (FD) methods yields the time dependent
system of ordinary differential equations (ODE’s):
56
57
am}
[Cl—at
” (”(1. + [Kly 1‘ [KlG) {U} - {F} = {0} (52)
where {U} is the set of the unknown nodal values, [K], is the stiffness matrix, arising
from the second partial of U with respect to x, [K]y is the stiffness matrix, arising from
the second partial of U with respect to y, [K]G is the stiffness matrix, arising from the
G term in the PDE (5-1), [C] is the capacitance matrix, arising from the time dependent
term in the PDE and {F} is the force vector, arising from the source term in the PDE
as well as the boundary conditions. The term [K] will be used to denote the matrix sum
[K],+[K],+[K]G. The matrices [C] and [K] are positive deﬁnite symmetric matrices.
[K] is singular before the boundary conditions are imposed. The matrices [C], [K], and
{F} are global and built from element contributions, using standard stiffness procedure
described by Segerlind (1984), whose coefﬁcients depend on the type of interpolation
function used to solve the problem, for example, linear or quadratic.
Analytical Solution of the System of Ordinary Differential Equation
The system in (5-2) can be solved analytically using modal analysis described in
Chapter 4. The technique is identical to the one-dimensional case and is not discussed
here. In one-dimensional problems, the analytical solution of the system of ODE
(AODE) was used to test the accuracy of the numerical scheme. Using this technique
beyond the developmental phase was not recommended. It requires large memory
storage space and intensive computations. In addition, the technique is not suitable for
problems where the initial condition is zero. Which are common in the heat transfer
58
literature. During the computations the initial condition end up in the denominator and
a division by zero causes the computations to cease. For these reasons, the AODE
solution was not included in the accuracy criteria used in two-dimensional problem
methodology below.
Numerical Solution of the System of Ordinary Differential Equation
Numerically, ﬁnite element or ﬁnite difference methods can be used to solve (5-2)
in the time domain. Although the ﬁnite element method shows clear advantages over
the ﬁnite difference method in the space domain this advantage does not extend to the
time domain, Segerlind (1984). There are many time domain schemes available in the
literature for each method. Many of these were introduced in Chapter two. Several
criteria are used to test these schemes including stability and oscillation. The current
study will explore ﬁnite element and ﬁnite difference space discretization and three time
discretization schemes corresponding to 0 values of: 0, 0.5, and 1; a total of six schemes:
1. Galerkin Forward Difference GFD: Galerkin ﬁnite element in space, forward ﬁnite
difference in time.
2. Galerkin Central Difference GCD: Galerkin ﬁnite element in space, central ﬁnite
difference in time.
3. Galerkin Backward Difference GBD: Galerkin ﬁnite element in space, backward
ﬁnite difference in time.
4. Forward Finite Difference FFD: Explicit ﬁnite difference in space, forward ﬁnite
difference in time.
59
5. Central Finite Difference CFD: Explicit ﬁnite difference in space, central ﬁnite
difference in time.
6. Backward Finite Difference BFD: Explicit ﬁnite difference in space, backward ﬁnite
difference in time.
The Objectives of the Two-Dimensional Study are:
1. Develop a time-step estimate for six numerical schemes given above that satisﬁes an
accuracy and stability criteria.
2. Develop recommendations to the solution scheme for the parabolic partial differential
equation based on the ﬁrst objective.
3. Compare and analyze the time step estimate of objective one with a commonly used
time step criteria such as Courant or Froude numbers.
4. Compare the accuracy of the two methods of space discretization, ﬁnite element and
ﬁnite difference.
METHODOLOGY
The methodology is brieﬂy described in the paragraph below and will be detailed
later in the section. An engineering problem governed by (5-1) and its boundary and
initial conditions is solved using one of the six schemes and several time steps. For each
time step the error between exact and approximate solution is computed. For this series
of solutions a response line of time step against error is generated. As the value of time
60
step increases the error is expected to grow. There exists a time step that will accurately
integrate the problem within a certain accuracy level. Time step below this value will
produce results that are below the deﬁned accuracy limits but too small for efﬁcient
computations. Time steps above this step produce error too large for accurate integration
and could induce oscillations (in case of Euler’s scheme). The numerical problem above
is solved using several grid sizes generating a relation between grid size and optimal time
step for the particular problem being solved. Results are summarized in regression
equations that are used for time step estimate. The procedure is repeated for the other
six schemes. The time step estimate equations were veriﬁed using a different set of
initial and boundary conditions to make sure the results are valid for a variety of
problems.
Finite Element and Finite Difference Comparison:
There is a major difference between the ﬁnite element and the ﬁnite difference
methods. Each method uses a different approximation or integration scheme to develop
the element matrices. The ﬁnite element method uses the Galerkin’s weighted residual
formulation and a set of Shape functions in the integration process. As a result of this
space integration the initial PDE is transformed to a system of ODE’s. A similar
integration is carried out in the time domain. The ﬁnite difference on the other hand
approximates the slope of the function between two nodes using the Taylor’s expansion
series. In ﬁnite difference, the space and time integration is done simultaneously so the
ODE’S intermediate step is seldom presented. The end result is that the element stiffness
61
matrix for each method is different.
In one-dimensional problems the ﬁnite element and ﬁnite difference methods have
the same element stiffness matrices. For this reason a discussion of the differences
between the two methods was not done. In two dimensions, the frnite element derivation
is given in many ﬁnite element books, Segerlind (1984), Reddy (1984). The derivation
of the ﬁnite difference element matrix is not a common presentation in the literature.
The element capacitance matrix for the ﬁnite element method over a rectangular
element is
I
I
Consistent
U
is
N
HN-liN
(c) l
[c l 36 1
L
(5-3)
J
Ohwr-‘N
1),/1 o
[c “’] — Lumped
ll
0
OO—C
CHOON#Nr—I
[HOG
where A is the element area and the element stiffness matrix for the ﬁnite element
method over a rectangular element is
2-2-11' '21-1-2‘
[k(‘)] = D: a '2 2 1 "'1 + Dy b 1 2 '2 “'1 (54)
6b -112-2 6a -1-221
1‘1‘22. -‘2'112.
where a and b are the side lengths of the rectangle. A square element is obtained by
allowing a=b and adding the two matrices in (5-4).
The ﬁnite difference stiffness matrix was derived by Segerlind (1989) and is
p
2 -10 -1'
[km] =2 ’1 2 ’1 0 (55)
2 0-12-1
_-10 -12j
The ﬁnite difference capacitance matrix is lumped and is identical to (5-3). The force
vector for both methods is zero when Q is zero otherwise the QL2 term is evenly
distributed over the four nodes.
Accuracy Ratio
The accuracy ratio is a measure that was deﬁned to quantify how well the
numerical solution approximates the exact solution. The accuracy measure used with the
two-dimensional problems was different than the measure used in the one-dimensional
problems. In the former analysis the hypothesis was that the numerical solution for the
system of ODE’s (NODE) could be as accurate as the analytical solution of the system
of ODE’s (AODE). The one-dimensional analysis proved that this hypothesis is good
when the number of nodes exceeds 7. For a grid that is coarser than 7 nodes in [0,1],
63
both the AODE as well as the NODE were inaccurate for all time step values. A 5%
error or less could be misleading because the solution could still be far from the exact
solution. The difference between NODE and APDE could be very large in the case of
a coarse grid irrespective of the time step. This fact is hidden in the accuracy ratio
deﬁnition (4-15). In addition, the APDE term is redundant Since it is being used in the
norm computations in the numerator as well as the denominator of the error ratio
formula (4-15).
This issue become more critical in two dimensional analysis. The minimum
number of nodes to yield accurate results increases geometrically when going from one
to two dimensions. Based on preliminary calculations over a square domain, numerical
integration results are accurate only when the total number of nodes exceeded 25 for any
time step value. In addition, the AODE solution was computationally demanding both
in time as well as storage and was difﬁcult to handle on a micro computer.
For those reasons, the AODE was eliminated from the accuracy term and the
NODE and APDE were compared directly. This change in the error criteria eliminated
the coarse grid from consideration and made the computations more efﬁcient.
The comparisons between the numerical ODE and analytical PDE solutions were
made using an error ratio deﬁned by
Z": 2'": | NODE, -APDE, |
NODE, 15-6)
j-l i-l
e:
mn
where NODE is the numerical solution for the system of ODE’s, APDE is the analytical
solution for the PDE, m is the number of sampling points in the space domain, and n is
64
the number of sampling points in the time domain. Other deﬁnitions could be used such
as the L2 norm instead of the L1 norm, that is, the absolute value of the difference of
NODE and APDE can be replaced by sum square of the difference.
Problem Statement
The unit step change problem was used in this section for the determination of
time step estimate equation. The unit step change problem is numerically more difﬁcult
to solve and it has the full set of eigenvalues in the APDE solution. This study is based
on the rationale that if the time step estimate produced accurate results for a step change
problem, that time step will be sufﬁcient for other problems when the grid and material
properties are the same and the time to steady state increases.
The step change problem is deﬁned as:
Initial conditions : u(x,y,0) = 1 ; t = 0
Boundary conditions : u(0,y,t) = u(1,y,t) = u(x,0,t) (5'7)
=u(x,l,t)=0;t20,0$x,y$l
In other words, the initial values of u are one. When t > 0 the values of u on
the boundary are set to zero.
The analytical solution of the above problem is:
"(1050’- 1—622 (2n+1)1(2m 1)([sin(2n+1)1rx][sin(2m+1)-uy])e401M?«an-0111558)
lln'o III-0
This information was generated from information provided by Ozisik (1980). Due to
Sylnmetry about x = y = 0, the problem could be solved over one quadrant of the
65
domain. The actual problem that was solved is deﬁned
ﬂu
an
u(0,y,t) = u(x,0,t) = 0 t > 0
=0 atx=y=1t20
(5-9)
u(x,y,t) =1 t =0
The analytical solution of the above equation is:
,: , ,1), smA, x srnAm y
= on on -(
u(x,y,t) 4E£e A,.)».
two m-O
(5- 10)
where A = M
" 2
Determination of Eigenvalues
By deﬁnition the eigenvalues of the system are the solution of
det{1K1 - A 161} = 0 ‘5‘“)
The eigenvalues, A’s, can be determined analytically for few problems and numerically
for any problem using the Jacobi. The Jacobi method will determine the complete set
of eigenvalues. Other methods can be used to determine only the ﬁrst and the last
eigenvalues which correspond to the smallest and the largest values.
Numerical computation of eigenvalues requires meshing the space domain and
building the global system of equations comprising of [K] and [C]. The total number of
eigenvalues equals the total number of nodes in the mesh minus the nodes with known
Values of u(x,y,t), which are the boundary nodes in this study. The Jacobi’s method was
used once for each set of boundary conditions using a coarse grid. During the
Computation the eigenvalues are sorted lowest to highest. The effect of the grid on the
66
lowest eigenvalue was found to be negligible. In addition, the ﬁnite element and ﬁnite
difference methods gave A, within 2% from each other. In the current problem A, was
evaluated using a grid of 25 nodes. The lowest eigenvalue is reported since it is needed
in the computation of the time to steady state. The largest eigenvalue was used in the
stability criteria in Euler’s scheme.
The lowest eigenvalue A, is a parameter that captures the changes in boundary
conditions as well as material properties. It does not change with initial conditions. The
time to steady state is not affected by initial conditions as long as the boundary conditions
are the same. For this reason A, is used in the time step estimate equation that are
discussed in a later section to characterize the different problems. A, for the current step
change problem is computed to be 4.76.
Time to Steady State
The response of the parabolic equation to a certain initial and boundary condition
is considered to have reached steady state when the variable u(x,y,t) does not change
with time. The time to steady state (t,,) was deﬁned as
4
t =_ 5-12
,, A] r 1
where A, is the lowest eigenvalue for the system.
This deﬁnition comes from the analytical solution of the PDE. At steady state
only the ﬁrst term contributes to the solution, and for e", 98% of the exponential term
representing steady state solution is reached. The time to steady state was used to
67
determine for how long the computation nwded to be continued as well as sampling point
location in the time domain. The time to steady state for the current step change
problem was 0.84 sec.
Time Domain Sampling Points
The sampling points for the time domain for the step change problem were:
0.08, 0.16, 0.24, 0.32, 0.4, 0.48, 0.56, 0.64, 0.72, and 0.8 seconds corresponding to
9.5, 19, 28.6, 38.1, 47.6, 57.1, 66.7, 76.2, 85.7, and 95.2% of t,, , respectively.
An exception was made for the time step = 0.1 where the sampling points were at:
0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 seconds corresponding to 12, 24, 36, 48, 60, 71,
83, and 95% of t“ , respectively.
Also when the time step was of 0.16 sec the sampling points were at 0.16, 0.32, 0.48,
0.64, and 0.8 seconds corresponding to 19, 38, 57, 76, and 95% of t,, , respectively.
Space Discretization or Griding
The step change problem deﬁned earlier was solved for six different grids. Each
grid was identiﬁed by the total number of nodes (N) over the entire domain x,y (0,1).
The nodes were equally spaced in both directions to generate square elements. The
square element was selected because this shape allowed an easy and direct way of
comparing the ﬁnite difference and ﬁnite element methods, and also the square element
i 3 considered the most accurate four node quadrilateral element. The square domain was
divided into 1, 2, 4, 8, 16, and 32 divisions in each direction giving rise to grids with
68
a total of 4, 9, 25, 81, 289, and 1089 nodes, respectively.
Time Step Fstimation Process
The posed problem was solved for each of the six meshes and the six numerical
schemes. For each of the schemes and mesh sizes a series of time steps were used and
a corresponding error value was generated. This value represented a point on a time step
versus error graph. The graph for a certain mesh and numerical scheme shows a time
step value below which there is no signiﬁcant gain in accuracy. This point was chosen
to be at 5% accuracy level. For each of the numerical schemes and mesh combination,
a set of accurate time steps or points exist. Since A, is discussed earlier to be a problem
dependent parameter and captures the variations among different problems and is
independent of the grid size, including it in any general time step estimate equation is
necessary. Therefore, the product of time step and A, was regressed against N for the
six schemes. These equations are meant to be used as a time step criteria to satisfy
numerical accuracy, stability, and oscillations. These equations are presented as a
problem dependent estimates and were tested against a different set of initial and
boundary conditions in order to prove their validity for other than the step change
problem.
69
RESULTS
The time step estimate regression equations for the six schemes are presented in
this section. These estimates will be evaluated in the next section.
The data generated from the numerical experiments correspond to optimal time
steps for each of the six grids and the six numerical schemes. These points were
generated from the numerically evaluated response curves similar to the one in Figure
5-1. This ﬁgure shows the response line of the error as a function of time step for the
ﬁnite element central difference scheme. The curve crosses the 0.05 error at a time step
of 0.017 seconds. Response lines were generated for all other mesh size and scheme
combinations, Figures 5-2 through 5-7. These ﬁgures show the response lines of the
error against time step for various meshes. Figure 5-8 show the computed optimal time
step as a function of N for the six schemes. The optimal time steps shown in this ﬁgure
were multiplied by A, which is evaluated to be 4.75 for the step change problem. The
data were then best ﬁtted by the power model: AtA, =aN’° where a and b are parameters
to be determined for each of the six schemes. These were determined using the linear
programming tool "Solver Function" available in MicroSoft Excel for Windows version
4.0. The sum of the differences between actual and estimated AtA, was reduced to zero
by changing the value of the two parameters a and b. The "Solver" iterates by changing
the initial values of a and b until it converts.
Another regression approach was followed by linearizing the above model to
log(AtA,) = log(a) + b log(N). The linear model was then regressed using Lotus 123
70
GALARKIN CENTRAL DIFFERENCE N = 1089
0.2
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
ERROR
0 0.005 0.01 0.015 0.02 0.025 0.03
TIME STEP
Figure 5-1. Changes of the error as a function of time step for GCD scheme and a grid
size corresponding to 1089 total number of nodes.
71
«0.0
.oEoﬁm 000 .595 moﬁm 52: get? not. 02m 2:: we c2850 a ma Coho 2: 06 80:20 .~-m oSwE
000.. uz
512 toll murz III-ll
I191 oouuz o
05.0
r
d
.35 02:.
o 5.0 v 5.0 N 5.0 5.0 000.0 000.0 25.0 «00.0 0
a 4 w w w w w O
./o .- Ed
-. 85
.. no.0
80883 39V83AV
1 v0.0
1. 00.0
-r 00.0
002000005 Om<>>m00 2.5.0.20
SEQ—um 000 60:: 8% :88 32:; 5.. 08m 050 06 5083 a we. Stu 20 we 895:0 .m-w 2:30
000—"..2 IIOIII OQNHZ o Fonz 11.011 mNuZ III-III
Lmhm us:
p .0 00.0 00.0 no.0 00.0 00.0 #00 00.0 N00 .00 0
t u A u a A u u . . 1.... o
n 111111111.- .11 . n 1... 8.0
\ .. 3
a
o -.. . u
0— 0 O
a
1 N0
5. 0N0
[ﬁr ”.0
0025—0020 0<~=ZmU 25.00.30
73
.oEon3. 000 50:: 8% :88 got? :8 02m 2:: 06 506:8 : m: Ste 2: 00 30525 Jim 2:30
O
IT ommnz 512 IIUI mmlz it;
..me 0.2:
05.0 v5.0 N50 50 000.0 80.0 wood «00.0 0
_ _ _ p _ _ 0 _
. A _ . . _ _ . o
8088!
002020005 0055—04; 2.0.5.20
2:22 000 :25: 8.0» :85 32:9 :8 02m 2:: .8 5088 m m: 8.06 20 06 80:20 .00 2:00
owe—"Z I‘ll omwuz IIAYII 5u2 o mwuz IDII 002 tell
my...» 0.2:.
o p .0 v — .0 N P .0 — .0 00.0 00.0 .00 No.0 0
.. n n a .- n _ . - o
.l
U\ n V a. 00.0
75
mUmemuEE 0:2.”— ._<~_._.ZmU
.8828 0:0 2: 08 Z .00 000050 0 00 020 2:: .0808 0800500 .w-m oSwE
_ 00$ lall. 00.0 IITII our. IT 090 cl 000 IT 006 Ill _
00002 ".0 mun—2:2
89 82 80 80 00v 8m 0
0 .. ._ . _ . ._ o
a WWW”? 6.0
. . .. 8.0
-- 00.0
,- 00.0
-- 00.0
77
«1315 SW“.
A
I
50.0
00015.2 .55 302.0 x5 00“. #5 0.00 1:3 zO:.<_~_<> $5 05:... mh<~50<
78
spreadsheet version 3.0 where log(a) and b were evaluated. Both approaches gave results
with insigniﬁcant differences. The grid with four nodes was excluded from the
regression because numerical integration of the ODE was at best 25 % in error compared
to analytical PDE irrespective of the time step used. This particular information would
not have been evident if an error criteria involving comparison of NODE and AODE had
been used. The error in the grid with four nodes was not due to numerical integration
in the time domain, it was due to space discretization. The regression equation each of
which for the six different schemes shown in Figure 5-8 are
Galerkin Forward Dzﬁerence At = Vii/€707 N 2 25 (5-13)
Galerkin Central Dtﬁerence At = NIT?” 2 25 (5-14)
Galerkin Backward Dzﬁerence At = xiii—501 N 2 25 (5-15)
Forward Finite Dzﬁerence At = _l'_19_ N 2 9 (5-16)
xlNlm
79
Central Finite Dzﬁerence At = _19_ N 2 9 (5-17)
*1 No.55
Backward Finite Dzﬁerecne At = {0% N 2 9 (5.18)
1
DISCUSSION
The time step equations as related to various problems and grids as well as the
Froude or Courant number are discussed in this section.
Finite Element and Finite Difference Methods
The element matrices for the two methods were quite different. Although the two
methods could be forced to match by selecting a regular grid pattern (Zienkiewicz 1977)
the intent here is to compare the two methods on a basis that they normally exist.
Comparison of the two methods suggests that the two are identical in terms of accuracy
as shown by the identical time steps for both schemes. This statement is a generalized
fact that might not apply exactly to all cases. The ﬁnite difference method showed some
superiority over the ﬁnite element method for coarse grids as well as for small At. This
feature did not carry on for ﬁner grids and larger time steps. In fact the differences were
not enough to show up in the regression equations. Although the regression equations
80
for the two schemes are very similar, the ﬁnite difference method showed acceptable
accuracy for 9 node grids while ﬁnite element method rejected this grid as being
inaccurate. In general the two schemes were identical and this is conﬁrmed in a later
section during comparisons with the Froude number. Figures 5-9 through 5-11 show the
computed optimal time step variation as a function of the number of nodes N for the
forward, central and backward difference schemes, respectively. Figure 5-9 shows that
the two forward difference schemes, GFD and FFD, compare well for large N. Some
differences were seen for coarse grids. The central differences schemes compared in
Figure 5-10 shows that GCD and CFD schemes produce identical results. Figure 5-11
shows a similar trend for the backward differences. The slight differences between the
GBD and BFD diminishes at large N. The ﬁnite difference method shows superiority
at small N for the backward difference scheme.
Comparison of Numerical Schemes
This section will highlight some of the major features of each scheme and how
it compare with the others. The central difference schemes show superiority in terms of
allowing a larger time step for any particular solution while still maintaining a high
accuracy. This is true for both the ﬁnite element and the ﬁnite difference methods,
Figures 5-12 and 5-13.
Comparing the two forward and the two backward difference schemes shows that
the backward difference schemes are more accurate for both the ﬁnite element and the
ﬁnite difference schemes, Figures 5-14 and 5-15.
com—
81
350000 000 0:0 000 08 move: 00 50:5: 05 00 020008 0 ma Swag—0 020 2:: 020000 389000 .06 2:30
out IIDII 00.0 lull
$002 “.0 mun—2:2
89 80 80 8v Bu 0
. _ b p o
—.|l u 0 _ _
-.. «00.0
.1 «50.0
d.
- 5.0
t «5.0
t 30.0
.. A20.0
05.0
.0 so
4318 3W“
1
1
moo—tus— mmhm 302....“ 02:. «Om mam 0:5 1:3 zO_._.<.~_<> ..um ms...— mh<~50<
.ononom 000 :5 000 .80 move: 00 50:5: 2: ...o 026:5 0 00 30:20 not. 0E: REES 08:0:50
ONE ID]. 000 +
$002 ".0 awn—2:2
8m: 80— 80 80 8v 8N
[ p 0 _ 0 p
_ q _ . 0 0
83
00015.2 “—0.5 0002.0 03.— 20”. 000 0&0 1:3 zO_.—<_~_<> my; 0.2:. w._.<~_00<
1
Y
L
r
I
50.0
000.0
000.0
v8.0
000.0
000.0
“00.0
000.0
000.0
5.0
.:-n 200E
dais ilNll
85
8:02.00 8:08.06 2.5.. 025 2.0 :8 83: 0o 80:5: 2.0 5.3 0:00:09, 0000 0:5 02:95.0 .m.-w 050.0
000 Iol 000 '0'. 000 .II-II
00002 “.0 000202
89 80 80 8v SN 0
u n t. a .r T u o
. . i/. 5.0
. mod
1 8.0
.. cod
.. nod
. 8.0
.. 8.0
.. 8.0
. ood
. ..o
«Ills 3WD.
005.0 2. 002000005 00.2.0 02<
05:... z. 00050.2 0.0.5 0.62.0 000:... 00.... 0N5 0.00 1...? zO...<_~.<> 00.5 05.... 0h<¢=0<
86
This could be due to the stability problem that the forward difference scheme
exhibits especially for ﬁne grids (large N). It was impossible from an experimental point
of view to separate the effect of stability on the forward difference schemes. Note that
Figures 5-14 and 5-15 show comparable optimal time step for the two schemes at small
N values
Sensitivity of the Error to Changes in Time Step
The response of the optimal time estimates to the different grid sizes changes with
the particular scheme. Central difference schemes, both ﬁnite element and ﬁnite
difference, show a clear point at which the slope changes. This is where a larger time
step produces a larger error. The curve below this point is very well behaved, Figures
5-1, 5-3, and 5-6. The backward difference schemes, both the ﬁnite element and the
ﬁnite difference methods, do not show any point of change in slope. Figures 5-4 and 5-7
show that the smaller the At the smaller the error. This could be explained by the fact
that central difference schemes are more accurate and reach a high level of accuracy at
relatively large At compared to the backward difference schemes. The Galerkin forward
difference scheme, GFD, in Figure 5-2 did not show a change in slope point where the
forward ﬁnite difference scheme FFD in Figure 5-5 does show. This could also be
explained by the same accuracy argument especially for small time steps although the
overall conclusion is both schemes are on the same accuracy level.
87
Effect of Element Size on Accuracy
The general trend that was observed in one-dimensional problems repeats itself
in the two-dimensional schemes. As the grid is reﬁned (larger N) the time step required
to maintain accuracy within a ﬁxed limit (5%) must be smaller, Figures 5-2 through 5-
15. The degree of sensitivity changes with the scheme. The forward difference
schemes, both the Galerkin and the ﬁnite difference formulations, are more sensitive to
grid size than the backward difference schemes. A major observation found is the
insensitivity of the backward difference scheme to the grid size. Figure 5-11 shows no
change in optimal time step as N grows. Unlike other schemes backward difference time
step equations did not show sensitivity to changes in N. This feature could be either at
a advantage for large N, or a disadvantage for small N.
Comparisons with Froude or Courant Number
The Froude number is a common criteria used in the selection of time step in
transient problems. It is mostly used in transport problems particularly in heat transfer.
Chapter three discusses some of the literature that uses this number as a criteria for
selecting time steps.
The Froude or Courant Number is deﬁned as
kAt
pch2
Froude N0. = (5-19)
where k, c, and p are given parameters in the PDE. The three variables are often take
as unity for simplicity.
882.8 000 0:0 000 .50 8.00: 00 808:: 2.0 5.3 0:38.83 0000 08.0 80:08.00 .30 0.80.".
000 IIDII 000 III
00002 00 002202
89 80— 80 80 00: Ba 0
O
p _
l
I
000.0
.500
l
I
88
nT .. 000.0
5.0
N 5.0
t v 5.0
o 5.0
l
I
4313 3W“.
1
I
l
j
l
I
05.0
002000005 00<3¥0<0 02< 0055000 20030.00 200—003—200
.8828 000 88 0.00 80 move: 00 82.8.8 05 53 acoumtg 08m 080 380800 67m 8sz
000 IT 000 III
00002 00 100552
80
n
d
80 _. 80— 80 8v SN
1P
-_
up
89
002000005 923 x0\nb)
where
tanhnb = H
taana = H
>‘n
3..
The ﬁrst 25 eigenvalues were evaluated using a simple Newton Raphson based
computer program. The summation up to 25 terms was not needed fort > 0, since most
of the contribution to the solution is in the ﬁrst six terms. This is particularly true at
large values of time. Nevertheless, the 25 terms solution improved ﬁtness of the solution
to the initial conditions. The value of H was taken to be 10.
96
First Eigenvalue and Time to Steady State Computations
A square grid of 25 nodes was generated to numerically evaluate the ﬁrst
eigenvalue. Different grids give similar values A, that was determined to be 3.4. The
time to steady state was then computed using (5-12) and was 1.18 sec.
Sampling Points
The error determination procedure was the same as in the methodology section.
The error is determined at speciﬁc intervals of time for all nodes in space. For the
derivative boundary problem the error was evaluated at 0.08, 0.16, 0.24, 0.32, 0.4,
0.48, 0.56, 0.64, 0.72, 0.8, and 0.88 corresponding to 6.8, 14, 20, 27, 34, 40, 47, 54,
61, 68, and 75 % of t,,, respectively.
Space Discretization
The following table gives information about the relation of the grid to actual
eigenvalues computed and the degrees of freedom that each grid exhibit. The number
of degrees of freedom is deﬁned as the number of independent unknown variables.
Nodes of known values or those duplicated from symmetry are not included. Only
results for the 81 node grid are presented in this section.
97
Table 5-1. Number of eigenvalues and degrees of freedom for the various grid sizes used
in the derivative boundary condition problem
II N No. of Eigenvalues Degrees of Freedom
P 4 4 3
9 9 6
25 25 15
81 81 45
289 289 153
1089 1089 1 122
Evaluation of Results
The derivative boundary condition problem presented above was solved over a
grid with 64 elements and 81 nodes for each of the six schemes. The time step estimate
was evaluated for each problem. Several runs were made for each scheme using
multiples and fractions of the time step estimate. Results are shown in Figures 5-19
through 5-24 for all of the schemes. Results of ﬁnite element and ﬁnite difference
schemes are identical and are discussed together.
The graphs in these ﬁgures are plotted as number of calculated time steps in the
x-axis and the error in the y-axis. The x-axis is computed as the actual time step used
98
Estimated time step performance under a DBC problem: GFD
method
0.03 ~~
0.025 '—
0.02 -*
Error
0.015 -~
0.01 v
0.005 -- .f
0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9
’ Number of calculated tlme steps
Figure 5-19. Error propagation with the time step for GFD scheme. Time steps are
reported as ratios of actual time step to the estimated one for the speciﬁc scheme
and grid used. DBC is abbreviation of derivative boundary condition.
99
Estimated time step performance under a DBC problem: FFD
method
0.014 ~-
0.012 Jr
0.01 --
0.008 0
0.006 --
0.004 ‘-
0.002 ‘-
Error
0 0.2 0.4 0.6 0.8 l 1.2
Number of calculated tlme step
Figure 5-20. Error propagation with the time step for FFD scheme.
100
Estimated time step performance under a DBC problem: GCD
method
0.9 ~-
0.8 --
0.7 4»
0.6 ~~
0.5 *r
0.4 *-
0.3 --
0.2 -r
0.1 -~
0 -———+—l——t-—l**lr
Error
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Number of calculated time steps
Figure 5-21. Error pr0pagation with the time step for GCD scheme.
101
Estimated time step performance under a DBC problem: CFD
method
0.9 --
0.8 ‘-
0.7 ‘-
0.6 w-
0.5 ~-
0.4 --
0.3 4.
0.2-1
0.1 *-
0 : I s I . ‘ .L .
.
Error
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Number of calculated tlme step
Figure 5-22. Error propagation with the time step for CFD scheme.
102
Estimated time step performance under a DBC problem: GBD
method
0.18 ~-
0.16 ~-
0.14 --
0.12 0
0.1 ..
0.08 ..
0.06 0 /
0.04 ~-
0.02 .. ./"./
O 1 2 3 4 5
Number of calculated tlme steps
Figure 5-23. Error propagation with the time step for GBD scheme.
103
Estimated time step performance under a DBC problem: BFD
method
0.16 r
0.14 ‘-
0.12 ‘r
0.1 -r
0.08 v
0.06 --
0.04 --
0.02 --
O : s s
Error
0 1 2 3 4 5
Number of calculated time steps
Figure 5-24. Error propagation with the time step for BFD scheme.
104
in the computation divided by the estimated time step drawn from (5-14) through (5-19)
for each of the respective schemes. A value of one means that the same time step as the
estimated value was used. The x-axis is scaled to each scheme and differs from one
scheme to the other since each scheme has a different estimated time step for the same
grid and problem speciﬁcation.
Forward Difference Schemes
Figures 5-19 and 5-20 show the error propagation with the time step for the
forward difference schemes (GF D and FFD schemes, respectively). Both schemes are
shown to be unstable for At greater than the estimated time step. For time steps less than
or equal to the estimated time step both schemes were very accurate. This indicates that
they are always accurate when stability criteria is met. Therefore, the stability
requirement would be a good criteria to use for an optimum time step when using the
forward difference scheme. Galerkin forward difference shows a point of change in
slope which is not present in forward ﬁnite difference. This was discussed in an earlier
section.
Central Difference Schemes
Both central difference schemes, GCD and CFD have a clear point where the
slope changes, Figures 5-21 and 5-22. These ﬁgures show that for At less or equal to
the estimated time step the schemes are accurate to within 5 %. As At grows beyond the
estimated time step, accuracy starts to decrease and the error increases above the 5%
105
mark. Both central difference schemes are very similar and have almost identical
responses.
Backward Difference schemes
Neither the ﬁnite element or ﬁnite difference formulations show any sign of
change in slope point, that is, the smaller the time step, the smaller the error. This was
observed in the one~dimensional study and was attributed to the lack of accuracy
inherited in the scheme. This is shown in Figures 5-23 and 5-24. Time steps equal or
below the estimated resulted in numerical solution within 5% error. For larger time
steps, numerical solution diverge from the reference solution and the error grow with the
value of the time step.
In summary, the validation procedure provides an evidence that the time estimates
work for each of the six schemes. The time step estimate were not conservative in
several of the methods, therefore, it is recommended that the values be round down to
nice clean numbers.
CHAPTER SIX
THE EFFECT OF SHAPE AND SIZE ON STABILITY
FOR TWO-DIMENSIONAL QUADRILATERAL ELEMENTS
USING THE MAXIMUM EIGENVALUE
This part of the study compares lumped and consistent formulations of the
capacitance matrix relative to their maximum eigenvalue and analyzes the effect of
element size and shape on the maximum eigenvalue for the two-dimensional four node
quadrilateral element. The primary focus is on the four node quadrilateral element, but
a comparisons is made with the three node triangular element and a recommendation of
which element to use in transient two-dimensional ﬁeld calculations is made.
MATHEMATICAL BACKGROUND
The analytical and numerical solution of a system of ordinary differential equation
requires a knowledge of eigenvalues for the eigenproblem
([K] - A [C] ){U} = 0 (6-1)
where [C] is the global capacitance matrix, [K] is the global stiffness matrix, {U} is the
unknown vector and A is an eigenvalue.
106
107
The above general eigenvalue problem has received limited discussion in
literature. Most of the mathematics literature, Boyce and Diprima(1986), Johnson et.al.
(1989), Leon (1986) discuss the standard eigenproblem
( [C‘] [K] - X [11 ){U} = 0 (6-2)
The discussion of the general eigenproblem centers around a transformation that converts
(6-1) into (6-2). From there on the problem is solved using the many available tools for
the standard eigenproblem.
The general eigenproblem of (6-1) arises in several areas of engineering, most
notably, mechanical and structural vibrations, elastic stability of continuous bodies and
heat transfer. The best discussion of the direct solution of (6-1) occurs in an engineering
oriented book by Bathe and Wilson (1976). The authors present several methods for the
direct solution of the eigenvalues that eliminate the nwd to transform it to the standard
problem.
The generalized Jacobi transformation technique was used in this study because
it calculates all the eigenvalues in a straight forward fashion that is easily implemented
on the digital computer. Computer code provided by Bathe and Wilson (1976) was
modiﬁed for use with a standard ﬁnite element program.
The largest eigenvalue of the generalized eigenproblem of (6—1) is needed for the
determination of stability and whether ampliﬁcation factor oscillations will occur. An
excellent method for estimating hm exists when using the ﬁnite element method. The
method and its proof are discussed by Fried (1979). The estimate goes as follows: Let
A“) be the largest eigenvalue for the element eigenproblem
108
([k‘°’] - >\‘°’[<>‘°’]){ll"”} (6-3)
then
M 5 Ann“) (6-4)
The largest eigenvalue for the element is an upper bound for the largest
eigenvalue of the entire system. Equation (6-4) is an excellent estimate when the
problem is solved using a uniform grid. It becomes very conservative when the grid
contains a mixture of element lengths (Fried, 1979). Equation 6-4 reduces the study
from the entire system to one element.
Stability and Oscillations
Stability is the study of the numerical behavior if a small error is introduced to
the solution scheme. If the error grows then the system is unstable. Conversely, if the
error remain bounded then the system is stable.
When oscillations occur the numerical value appear on both sides of the correct
solution. The oscillations criteria for the single step methods is
Ats — (65)
where c = l, 2 and 3 for 0 = 0, 0.5, and 0.6667 respectively. The backward
difference scheme never oscillates.
The factors affecting Am are mesh geometry (size and shape), initial and
boundary conditions, and the material properties. For a given problem, the only
parameters that can be changed are those pertaining to the element geometry, changing
109
the other parameters changes the problem under consideration. The question is how to
keep km at a small in order to avoid numerical difﬁculties.
The element matrices [k] and [c] are easily evaluated for a triangle, rectangle and
square. The analytical evaluation for a four nodal quadrilateral element is not possible.
Gauss-Legendre numerical integration scheme is required, Segerlind (1984) and Bathe
and Wilson ( 1976). The following listing gives the [CM] and [km] matrices for the
various elements used in the current investigation:
For a triangular element [01"] is
2 1 1
Consistent: [em] = BIL/1 1 2 l
1 l 2 (6-6)
1 O O
Lumped: [c“’] = 93—14 0 l 0
0 0 1
For a rectangular element [CM] is
110
“4212'
36
L2124_
DtA
.p
EIS’
Consistent: [c “’] = '
Lumped: [c“’] = —
The stiffness matrix for a triangular element is
D b} b, bl. b, b,
[k“’] = ﬁ b, b} bf bj b,
b, b, b] b, bf
where:
m=&n-&u
aj = kai ’ xiYk
at "' Xin ' iji
bi = Yj ' Yk
bi = Yk ' Yi
bk = Yi ' Yj
Ci = 'Xj + Xk
CI = -Xk ‘1' xi
ct = 'Xi + Xj
and
- 1 X, Y,-
A = det 1 X} Y}
_ 1 X, Y”
Km and Ym are x and y coordinates of the nodes i, j, and k, respectively.
2421
1242
'1000'
0100
0010
_0001‘
(6-7)
(6-8)
111
The stiffness matrix for a rectangular element is
'2-2-11' '21-1-2‘
[km] = D: a ’2 2 1 "1 + Dyb 1 2 "2 "1 (6'9)
6b -112-2 60 -1-221
L1-1-224 .‘2‘112.
where a and b are the sides of the rectangle. A square element is derived by allowing
a=b.
METHODOLOGY
A numerical experiment was designed to study the effect of capacitance matrix,
and element geometry on hm. In addition, practical examples were solved to gain an
understanding of the impact of meshing on Am. The numerical experiment that lead to
deﬁning the research objectives is presented in this section.
A FORTRAN code computer program by Bathe and Wilson was modiﬁed for use
to determine the eigenvalues for the various scenarios. The respective eigenvalues were
used as the criteria for comparing the various scenarios. The goal was to reduce km, by
changing the element shape. Elements with lower eigenvalues are superior than those
with higher eigenvalues.
The parameters involved in the experiment were the shape and size of the four
node quadrilateral element, and type of capacitance matrix. Triangular elements were
modeled as three and four node elements. The fourth node in the later element was
ﬂoating along one of the sides. More detailed discussion on the triangular elements is
112
in a later section.
The effect of area on Am was studied by selecting two element shapes, a square
and a trapezoid. For each of the two element shapes, the two capacitance matrices, and
various element sizes, Am was computed.
The effect of the element shape on Am was also studied. For a ﬁxed element
area the shape was varied and the corresponding AM was computed. Two different
element areas were considered, 1 and 4 area units (L), and two [01"] matrices, lumped
and consistent.
Two of the practical examples commonly encountered in the mesh generation
process were considered. Those pertain to dividing a certain element into two or more
elements. Original element geometries of rectangular and trapezoidal shapes, and further
dividing these elements into triangles and squares is investigated. Finally, modeling a
three node triangle using the four node quadrilateral is analyzed.
Shape Parameter (Sp)
In order to study the effect of element size on the maximum eigenvalue the area
was taken as a strong measure of the size. There is a need to develop a measure for the
element shape on the basis of which elements are compared. The shape parameter, Sp,
Figure 6-1 was found to be a useful parameter and was used as a measure of the element
shape. Figure 6-2 displays four different shapes and their Sp values. A square has the
least Sp value among the four nodal quadrilaterals. As the square is deformed, its angles
are no longer 90" and its Sp value becomes larger.
113
Shape Parameter
m
Sp = (im/ik) + (mj/ki)
i .
Sp = (ratio of sides) + (ratio of diagonals)
All ratios are largest over smallest
Figure 6-1. Shape parameter deﬁnition for a four node quadrilateral element.
114
Shape Parameter Values for Unit Area Shapes
U)
'0
II
N
(D
'0
II
(D
(D
'5
ll
0
Figure 6—2. Shape parameters values for some selected quadrilaterals
l 15
RESULTS AND DISCUSSION
The effect of the capacitance matrix, size and shape of the element on Am is
discussed here.
Effect of Element Size
The effect of element size on the maximum eigenvalue is shown in Figures 6—3
and 6-4. Figure 6-3 shows a plot of the area against the maximum eigenvalue for a
square element and the two [01"] formulations. The ﬁgure shows that the maximum
eigenvalue varies as a linear reciprocal function of the element area. This is true for
both lumped and consistent formulations. As the area increases the maximum eigenvalue
decreases which will increase the time step requirement to satisfy the stability
requirement of (6-5).
Figure 6-4 shows the same trend as Figure 6-3 for a trapezoidal element. The
slope of the line is still -1 which yields the relation
kw = f(Area)" (6-10)
The regression equations for the lines in Figures 6-3 and 6-4 are shown in (6-11).
25 .70 A " : Square Consistent
03.98 A‘1 : Square Lumped
m“ 48.60 A " : Trapizoid Consistent
£08.24 A '1 : Trapizoid Lumped 1
(6—11)
where A is the element area.
116
Area Effect on the Largest Eigenvalue
Square Element
1000
Eigenvalue
r l
‘-
o'oo1r 1 11111111 1 11111111 1 llllllll I llllllll l lllllllI
0.1 1 10 100 1000 10000
Area
Figure 6-3. Effect of area on the largest eigenvalue for a square element
117
Area Effect on the Largest Eigenvalue
Trapezoidal Element
Eigenvalue
10000
0.1
0.01
0.001 L 1 1111111 1 11111111 1 11111111 1 1 Lllllll 1 11111111 J 11111111 1 1111111
0.001 0.01 0.1 1 10 100 1000 10000
Area
Figure 6-4. Effect of area on the largest eigenvalue for a trapezoidal element
118
Effect of Element Shape
The effect of element shape on the maximum eigenvalues is shown in Figures 6-5
through 6-7. Figure 6-5 plots the shape parameter Sp against the maximum eigenvalue
for a unit area and the two [01"] formulations. Figure 6-5 shows that maximum
eigenvalue increases with Sp. The maximum eigenvalue for the lumped formulation
shows less sensitivity to Sp than the consistent formulation. On the other hand, the
maximum eigenvalue for the consistent formulation varies linearly with Sp. The
difference between the two formulations increases as Sp gets larger. Figure 6-6 shows
the same effect for an area value of four. Figure 6-7 is a comparative chart for some
selected shapes all with a unit area. For all elements with four nodes, the square has the
lowest eigenvalue. As the shape becomes deformed and its Sp value gets larger, the
maximum eigenvalue increases. That increase is considerable and affects the behavior
of the entire system. The triangular element is a unique case and could plot any where
in the range of shapes presented.
The regression equations for the lines in Figures 6-5 and 6-6 are shown in 6-12).
6.02 + 9.24 Sp : Area=1 Consistent
x = 2.77 + 1.38 Sp : Area=1 Lumped (5.12)
m 1.51 + 2.31 Sp : Area=4 Consistent
0.69 + 0.34 Sp : Area =4 Lumped
The regression equations show the same trend in both element areas. The
correlation coefﬁcient for the linear regression of the lumped formulation curve for both
areas was relatively low (R2 = 0.53) compared to the consistent formulation (R2 =
0.99). The actual results indicated that the maximum eigenvalue for the lumped
119
"at
Shape Effect on the Largest Eigenvalue
Unit Area
0 Eigenvalue
Figure 6-5. Effect of shape on the largest eigenvalue for a unit area
120
Shape Effect on the Largest Eigenvalue
Area of 4
Eigenvalue
1o- ........................................................................................................................
Shape parameter
"N
Figure 6-6. Effect of shape on the largest eigenvalue for an area of four
121
Comparative Chart for Selected Shapes
Unit Area
Eigenvalue
Square Trapeziodz Rhombus Trapezoid 4 Trapezoid 3 Rectangle trapezoidt
Element Type
DUI
Figure 6—7. Comparative chart for some selected shapes as to their maximum eigenvalue
122
formulation did not grow for Sp values beyond 15 as was indicated earlier.
Effect of [cm] Formulations
For all geometries that were studied, the lumped formulation had a lower
maximum eigenvalue than the consistent formulation. The difference between the two
formulations is signiﬁcant in all cases especially for large values of Sp as shown in
Figures 6-5 and 6-6.
The lumped formulation is also shown to be less sensitive to changes in element
shape. As the element shape is deformed, the maximum eigenvalue changes more slowly
than the change with the consistent formulation. This is another advantage to the lumped
formulation. The maximum eigenvalue for the consistent formulation increases linearly
with Sp. The lumped formulations results were assumed linear although there was no
increase in km beyond Sp of 15. The linear ﬁt of the lumped was very poor compared
to the consistent formulation.
Practical Examples
Figures 6-8 and 6-9 present a common problem in the grid generation process.
Figure 6-8 shows a rectangle that is split into two squares each of which is further split
into two triangles. The consistent formulation results shows that splitting the rectangle
into the squares will increase the maximum eigenvalue by about 30%. Further splitting
of the square into two triangles will reduce the maximum eigenvalue by about 30%
123
Dividing a Rectangle into
Squares and Triangles
(1,1)
(0,1) - (2.1)
(0,0) (1.0) (2.0)
Max. Eigenvalue
Consistent Lumped
Rectangle 16.09 4
Square ' 25.75 4
Triangle 12 9
Figure 6-8. Practical example 1 : dividing a rectangle into squares and triangles
124
‘ Dividing a Trapezoid into Triangles
(1,2)
(0.1) (1.2.1)
(1.0)
Max. Eigenvalue
Consistent Lumped
Trapezoid 43.75 7.422
Triangle 1 1.7 6.407
Figure 6-9. Practical example 2 : dividing a trapezoid into triangles
125
below the original AM of the rectangle. That indicated that when using a consistent
formulation, cutting a rectangular element into squares is not favorable. On the other
hand cutting that square into two triangles is desirable. Therefore, when using the
consistent formulation, triangles close to an equilateral shape are more desirable than a
square element. The square element remains the best of the four node elements.
The lumped formulation has the opposite results. Cutting the rectangles into
squares did not change the maximum eigenvalue. This could be explained by the fact
that reducing the area was offset by the enhancement in the element shape. Cutting the
resulting square into triangles increases the maximum eigenvalue by more than double.
Therefore, when using the lumped formulation, square elements are the desirable
element shape. Changing these squares into triangles is not desirable. The only
exception is when equilateral triangles are used. The later are shown to be the best
elements for reducing the maximum eigenvalue. A three node equilateral triangle with
a unit area has a km of 3.36 which is around 15% less than a four node square element.
This will in fact be in favor of equilateral triangles as a superior meshing element.
Figure 6-9 presents a similar problem encountered on curved boundary. For both
consistent and lumped formulations, dividing the trapezoid into two triangles reduces the
maximum eigenvalue by about 75% and 15%, respectively. Therefore, it is
recommended to use triangles over a curved boundaries on the condition they be as close
to equilateral triangles as possible.
This is explained by the fact that no element having interior angles exceeding 90°
is desirable. This result is in agreement with Segerlind (1984).
126
What is Unique About a Triangle
A triangle can be approximated as a four node quadrilateral. In that case there
is a fourth node that is ﬂoating and needs to be deﬁned for the evaluation of element
matrices. Figure 6—10 describes a case where the triangle shown has been treated as a
three node triangular element and as a four node quadrilateral element. Equation 1 has
three and four solution points for three and four nodes element, respectively. This by
itself brings an advantage of triangular elements over a four node element. The above
conclusion is clear looking at the results shown in Figure 6-10. The maximum
eigenvalue for a three node triangle is less than the maximum eigenvalue for a
quadrilateral element. This is true for lumped as well as consistent formulations.
127
Comparison of a Quadrilateral and a Triangular Element
(0.4)
Max. Eigenvalue
Element Type Consistent Lumped
Oaudrilateral 73.04-33.32 9.6-4.55
Triangular 6.45 1.614
(0.0) (2.0)
In the quadrilateral case the hypotenous
has an extra ﬂoating node
Figure 6-10. Comparison between a quadrilateral and a triangular element
CHAPTER SEVEN
SUMMARY AND CONCLUSIONS
This study presented a functional estimate for the time step to be used in the
numerical solution of diffusion problems using the single step schemes. The estimate
equations were developed for the unit step change function through a numerical
experimental procedure by varying the problem mesh size and ﬁnding the optimal time
step that accurately integrated the problem. These points were ﬁtted by a power function
and a regression equation determined the unlarown coefﬁcients. The process was
repeated for three schemes; the forward difference, central difference, and backward
difference in time. A coarse grid was generated over the solution domain to determine
the numerical value of the lowest eigenvalue. This parameter was used to estimate the
time to steady state and time domain sampling points used to collaborate an error term.
The time step estimate equations were then checked against other problems not used in
its development.
The time step estimates developed in this work were correlated to a well known
criteria, the Froude number. A major observations made during this correlation process
was that the Froude number is problem independent. It only changes with the mesh size
and material property. On the other hand, the time step estimate of this study changes
for each problem boundary condition as well as mesh size and material properties. These
128
129
later characteristics were captured by the smallest eigenvalue. Therefore, any particular
correlation between the Froude number and the time step estimates changes with the
problem being solved. A regression ﬁt between the two time step estimates for various
problems ( M) was presented.
As a conclusion for this work, the most important observation was that the two-
dimensional results were an extension of the one dimensional analysis. In fact, most of
the ﬁndings and observations found to be true for the one—dimension problem were also
true for the two-dimensional problem.
The time step estimate equations for the one-dimensional problems are
Euler Scheme At = .27N-l.6
Al
Central Difference 1.13N-1Js
At = _
xl
Galerkin Scheme 70 N -3.79
At =
A1
Backward Difference At = 30.6N‘3‘9‘
Al
In each case N is the number of nodes and should be greater than seven
130
The time step estimates for the two-dimensional problems are
Galerkin Forward Dtﬁerence At = 1'8 N 2 25
x1 N104
Galerkin Central Dtﬁerence At = 1'6 N 2 25
x1 No.55
Galerkin Backward Dzﬂerence At = T1275“ 2 25
1
Forward Finite Difference At = L?— > 9
x1 N101
Central Finite Dtﬁ'erence At = _1_°6_ 9
x1 No.55
. . . _ 0.05
Backward Finite Dlﬂerecne At - A N“ N 2 9
l
The most important conclusions from the one-dimensional study are:
-Time step criteria equations are valid for problems other than those included in their
development. The estimated time step value ensures accuracy as well as stability
for all problems tested. The use of each equation should be limited to the mesh
sizes that fall in the region under which they were developed. For most problems
found in the literature this region seems to be inclusive.
- For a ﬁxed number of elements, the central difference schemes is less sensitive to
increasing time step than the other schemes.
131
For a ﬁxed time step the central difference scheme is less sensitive to increasing the
number of nodes.
- For a particular problem and a certain accuracy level the central difference scheme
allows a larger time step than the other schemes
The backward difference scheme is no better than Euler’s scheme. For larger time
steps where Euler’s scheme is unstable the backward difference scheme is
inaccurate.
- The Froude number that gives accurate results changes with the boundary conditions
as well as the single step scheme used.
The ﬁner meshes (N = 16), produced smaller error (compared to N =8), provided the
proper time step values were used.
The specific conclusions for the two-dimensional study are:
Analytical solution of the system of ordinary differential equations is not practical
calculation in two dimensions. This is due to the large size of the matrices and
the number of computations involved.
Applying the element matrices and stiffness procedure to the ﬁnite difference method
increases the efﬁciency of the method computations and allowed using the same
code that is used in the ﬁnite element analysis. This gives the ﬁnite difference
method the same advantages that have been attributed to the ﬁnite element method
such as: writing a general code (problem independent code), ease of handling
boundary conditions, etc.
132
- The error ratio that was used for the one—dimensional problems could not be used with
two-dimensional problems. Another error criteria was developed. The former error
ratio formulation assumed that the numerical solution of the system of ordinary
differential equations could be at best equivalent its counterpart the analytical solution of
the system of ODE. The danger of the approach is in using it for coarse grids where the
analytical ODE solution is not acceptable due to large space discretization error.
- The time step estimate equations prove to work satisfactorily for problems that were
not used in the determination of the estimate.
- The two numerical methods, ﬁnite elements and ﬁnite differences gave time estimates
that are very close to each other suggesting a similar accuracy. The ﬁnite
difference method was more accurate for a coarse grid. This slight superiority
did not carry for ﬁner meshes and was not signiﬁcant enough to be detected in
the regression equations.
- The central difference scheme was the most accurate single step scheme.
The backward difference scheme were more accurate than forward difference scheme.
This superiority was not signiﬁcant and was attributed to the instability of the
forward difference scheme for many of the time step values.
The backward difference scheme time step estimate shows less sensitivity to the grid
size than the other schemes.
Central difference scheme shows a clear point where the slope changes rapidly in the
optimal time step vs. grid size curve. The other schemes, namely the backward
and forward differences, did not. This was attributed to the accuracy of the
133
scheme.
- For all schemes the optimal time step changes with the mesh size. The rate at which
it changes varies with the scheme.
- Comparison of the time step estimate equations with the Froude number shows
that for each problem there exists a Froude number equivalence that will give a
time step equivalent to the presented time step estimates. This equivalence
changes as A, of the problem changes.
- Central difference schemes exhibit large variations in Froude number equivalence
corresponding to different h, compared to the other schemes.
The element size and shape are two important mesh parameters that affect the
dynamic behavior of numerical solutions of transient systems. Some of conclusions of
this part of the study are:
- In all cases studied, the lumped formulation had a lower ism than the consistent
formulation.
- The shape parameter, Sp, was a good measure of how the shape of the element
behaves because hm grows with Sp, km for consistent formulation grows
linearly with Sp. The km for the lumped formulation is less sensitive to the
element shape. As the shape is deformed, hm does not increase as much as the
km for consistent formulation does.
- hm = f(A"), as the area increases Am decreases and the minimum acceptable time
step gets larger. This is true for both lumped and consistent formulations.
- For a four-node element with a ﬁxed area, the square has the lowest hm. Depending
134
on the element shape a three-node triangle may have a lower hm than a square.
- Dividing an element into two or more elements can increase or decrease hm,
depending on type of formulation used for the capacitance matrix and the
original shape.
RECOMMENDATIONS
Based on the results presented above the following major recommendations are made:
Use the time step estimates presented for the particular scheme when solving the
parabolic diffusion problem.
Use either ﬁnite element or ﬁnite difference methods for space discretization.
Use central difference scheme for time discritization.
FUTURE RESEARCH
The following research needs to be addressed in the future:
- Study the effect of element shape on the time step.
Use the methodology to ﬁnd regression equations for three-dimensional and
axisymmetric problems.
REFERENCES
Abdalla, H. and R. Paul Singh. 1985. Simulation of Thawing of Foods Using Finite
Element Method. Juumﬂ gf Fug! Prmess Enginﬁring. 7: 273-286.
Alagusundaram, K., D.S. Jayas, and WE. Muir. 1991. A Finite Element Model of
Three-Dimensional CO2 Diffusion in Grain Bins. Paper presented at the 1991
International Winter Meeting of the Amsrjgu Smisty uf Agricultural Engingrs,
St. Joseph, MI.
Allaire, P. E. 1985. i ini El men
Tgausfer, auu Fluiu Mghauius. Wm.C. Brown Publishers.
Baker, J. (1993). Personal communication.
Bathe, Klaus-Jurgen and Edward L. Wilson. 1976. m ri M h in Fini El m n
Analysis. Printice-Hall, Englewood Cliffs, New Jersey.
Blacker, T.D., M.A. Stephenson, J.L. Mitchiner, L. R. Phillips and Y.T. Lin. 1988.
Automated quadrilateral mesh generation: A Knowledge System Approach.
Amsrigu Sguigy gf Mghauiual Enginﬁrs paper 88-WA/CIE-4.
Boyce E.W. and RC. Diprima. 1986. Elementary Diffsrenuﬂ muau'gns. Fourth
Edition. John Wiley & Sons, New York.
Bykat, A. 1976. Automatic Generation of Triangular Grids: I-Subdivision of a General
Polygon into Convex Subregions. II-Triangulation of Convex Polygons.
Intematignal Juurnal fur Numgriﬂ Methods in Engineering Vol. 10, 1329-1342.
Cavendish, J .C. 1974. Automatic Triangulation of Arbitrary Planner Domains for the
Finite Element Method. Intemau'gnal Jgumal fur Numeriual Methﬂs in
Engingring Vol. 8, 679-696.
Celia, M. A., Bouloutas, E. T., and Zarba, R. L. 1990. A general mass-conservative
numerical solution for the unsaturated ﬂow equation. W r R r R h,
25(7), 1483-1496.
Cheney, W. and D. Kincaid. (1985). Numsriﬂ Mamsmau'us aud Campuu'ng. Second
Edition. Brooks/Cole Pub. Co. Paciﬁc Grove, Ca.
135
Churchill, Ruel V. and Brown, James Ward. 1987. W
Emblems, fourth edition. McGraw-Hill Book Company, New York.
Chu, S. T., and Fustrulid, A. 1968. Numerical Solution of Diffusion Equations.
Tmsaetious at the ASAE, 705-715.
Cleland, A.C. and R.L. Earle. 1984. Assessment of Freezing Time Prediction
Methods. W. Vol. 49.
Cooley, R. L. 1983. Some new procedures for numerical solution of variably saturated
ﬂow Problems. MW, 12(5). 1271-1285.
Dhatt, Gouri and Gilbert Touzot. 1984. The Finite Element Methetl Displayﬁ. John
Wiley and Sons.
De Baerdemaeker, J ., Singh, R. P., and Segerlind, L. J. 1977. Modelling heat transfer
in foods using the ﬁnite-element method. leumal ef Feed ﬂuxess Engingring,
37-50.
Elnawawy, O. A., Valocchi, A. J., and Ougouag, A. M. 1990. The cell analytical-
numerical method for solution of the advection-dispersion equation: Two-
dimensional problems. W r R r R h, _2_6(11), 2705-2716.
Fong, F. K., and Mulkey, L. A. 1990. Comparison of numerical schemes for solving
a spherical particle diffusion equation. W r R r R h, 26(5), 843-
853.
Franca, A.S., K. Haghighi, L. Oliveira, and E. Kang. 1992. Adaptive Finite Element
Analysis of Transient Heat Conduction Problems. ASAE Paper 92-3525.
Presented at the 1992 m ri i ri 1 En in r Winter Meeting
Nashville, Tennessee.
Fried, Isaac. 1979. ri l ' f iff r n ' ' n . Academic Press, New
York.
Gear, G. Williams. 1971. Numerieal lnt'u'al Vﬂue Preblems in minaty Diffetenu'al
Eguau'gus. Printice-Hall, Englewood Cliffs, New Jersey.
Giannakopoulos, A.E. and A.J. Engel. 1988. Directional Control in Grid Generation.
1011005239an0901.0010: 74. 422-439-
Goudreau, G.L. 1970. Evaluau'eu ef Diserete Methtuls fer the Liugt Qynamie Resugng
f El ' Vi l ' ' . UC SESM Report 69-15. University of
California, Berkley.
136
137
Gresho, P.M. and R.L. Lee. 1979. ”Don’t Suppress the Wiggles - They’re telling you
Something", in Finite Element Methods for Convection - Dominated Flows,
AMD-Vol.34, Edited by T.J.R. Hughes. Am, See, ef Mggh, Enggs., New York.
Haghighi, Kamyar and Larry Segerlind. 1988. Modeling Simultaneous Heat and Mass
Transfer in an Isotropic Sphere-A Finite Element Approach. 11201150911905.9110:
Am ri ' f ' l E in r , Vol. 31(2): 629-637.
Henrici, Peter. 1977. Errer Meagatlun fut Differeng Methgus. Krieger Pub. Co.
Huntington, New York.
Huebnes, Kenneth H. 1975. jljhe Finite Element Memm fer Engt'ueers. Wiley
Interscience, Nwe York.
Hughes, Thomas J.R. 1987. Th Fini El m n M Lin
WM. Prentice-Hall, Englewood Cliffs, New Jersey.
Irudayaraj, J. 1991. Moving Evaporative Front Prediction Using the Finite Element
Method Paper presented at the 1991 International Winter Meeting of the Amerieau
Seeiety ef Agrieuluttd Engineers, St. Joseph, MI.
Irudayaraj, J ., Haghighi, K., and Stroshine, R. L. 1990. Nonlinear ﬁnite element
analysis of coupled heat and mass transfer problems with an application to timber
drying. Dtying Tghnelugy, 8(4), 731-749.
Jaluria, Y. and K. Torrance. 1986. W. Hemisphere Pub. Co.
New York.
Johnson, L.W., R.D. Riess, and LT. Arnold. 1989. Innueu'en tg Linea; Algebta.
Second Edition. Addison-Wesley Pub. Co.
Kam Liu, W. Belytschko, T., and Fei Zhang, Y. 1984. Partitioned rational runge kutta
for parabolic systems. lnternau'enal leumal fer Numerical Methtﬂs in
Engingring, 20, 1581-1597.
Kang, E. and K. Haghighi. 1992. A Knowledge-Based A-Priori Approach to Mesh
Generation in Thermal Problems. MW
in Engingring Vol. 35, 915-937.
Kreyszic, E. (1962). Atlvauegl Engingring Mathemau'es, John Wiley & Sons, New
York.
Lapidus, Leon and George F. Pinder. 1982. Numerigl Seluu'en ef Pasu'al Differential
i n in i n En i rin . John Wiley and Sons, New York.
138
Leij, F. J ., and Dane, J. H. 1990. Analytical solutions of the one-dimensional advection
equation and two- or three-dimensional dispersion equation. WateLRestuuces
Research. 26(7). 1475-1482.
Leon, SJ. 1986. W. Second Edition. Macmillan Pub.
Co. New York.
Lewis, R.L., A.S. Usmani, and LT. Cross. 1991. Adaptive Finite Element Analysis
of Heat Transfer and Flow Problems. Neullu, Camp, Meeh,-State ef the Ag.
Prof. Stein Memorial Volume, Springer-Vertag.
Lewis, Roland W. and John C. Brunch, Jr. 1974. An Application of least Squares to
one-Dimensional Transient Problems. lnt, leur, fer Numeriﬂ Methguls in
Engg., Vol. 8, 633-647.
Lynch, D. R. 1984. Mass conservation in ﬁnite element groundwater models. Aduaums
iu Water Reseurgs, 7, 67-75.
Maadooliat, Reza. 1983. Element autl Time Step Criteria fer Selving Time Demntlent
Field Preblems Using The Fiuite Element Methed. Unpublished Ph.D.
Dissertation, Michigan State University.
Mavriplis, D. J ., Jameson, A., and Martinelli, L. 1989. Mulu'gu'gl selutien ef the
MW Paper Presented at the 27th
Aerospace Sciences Meeting, Reno, Nevada.
Misra, R. N ., and Young, J. H. 1980. Numerical solution of simultaneous moisture
diffusion and shrinkage during soybean drying. Ttausaeu'ens 9f the ASAE, 1277-
1282.
Misra, R. N., Young, J. H., and Hamann, D. D. 1981. Finite element procedures for
estimating shrinkage stresses during soybean drying. Tt'ausaetiuns ef the
WWW 751-755.
Mohtar, R. H. and L. J. Segerlind. 1991. Effect of Element Geometry on the
Eigenvalues. Paper No. 91-7551. Paper presented at the 1991 International
Winter Meeting of the American Society of Agricultural Engineers, St. Joseph,
MI.
Mohtar, R. H. and L. J. Segerlind. 1992. Time Step Criteria for Solving Unsteady
Engineering Field Problems. Paper No. 92-3523. Paper presented at the 1992
International Winter Meeting of the American Society of Agricultural Engineers,
St. Joseph, MI.
139
Myers, GE. 1977. The Critical Time Step for Finite Element Solutions to two
Dimensional Heat-Conduction Transients. ASME paper 77-WA/HT-14. Am,
W New York.
Myers, Glen E. 1971. WWW. McGraw-Hill,
Inc, New York.
Narasimhan, T. N. 1978. A perspective on numerical analysis of the diffusion equation.
Atlvaugs in Water Reseurees, 1(3), 147-155.
Ortega. J. M. 1990. WWW. Reprinted by the Society
for Industrial and Applied Mathematics part of the Classics in Applied
Mathematics Series (No. 3)
Ortiz, M., and Nour-Omid, B. 1986. Unconditionally stable concurrent procedures for
transient ﬁnite element analysis. WWW
Engineering, 58, 151-174.
Ozisik, Necati M. 1980. Hgt Cendueu'en. John Wiley & Sons, New York.
Patankar, Suhas V. 1980. Numerieal Heat Ttausfer mu Fluitl Elew. Hemisphere Pub.
Co. New York.
Patankar, Suhas V. 1991. m i n f n ' n Fl w T f r.
A textbook featuring the computer program CONDUCT. Innovative Research,
Inc. Maple Grove, Minnesota.
Peraire, J ., Peiro, J ., Formaggia, L., Morgan, K., and Zienkiewicz, O. C. 1988. Finite
element Euler computations in three dimensions. MW
Numerigl Methms in Engingring, 26, 2135-2159.
Potter M. and J. Goldberg. (1987). Mathemau'eal Methﬂs. Second Edition. Prentice
Hall Inc. Englewood Cliffs, New Jersey.
Powers, David. 1987. MW, Third Edition. Harcourt Brace
Jovanovich, Pub. Academic Press.
Reddy,J.N. 1984. An n ' h Fini Elm n M . McGraw-HillBook
Co. New York.
Rigal, A. 1990. Numerical analysis of three-time-level ﬁnite difference schemes for
unsteady diffusion-convection problems. WWW
Methegs in Engineeg'ug, 39, 307-330.
140
Roberts, David L. and M. Sami Selim. 1984. Comparative Study of Six Explicit and
Two Implicit Finite Difference Schemes for Solving One-Dimensional Parabolic
Partial Differential Equations. Intematienal leumal fer Numerieal Methms in
Engingring, 2Q, 817-844.
Rushton, K. R. , and Tomlinson, L. M. 1971. Digital computer solutions of groundwater
ﬂow. W. 12. 339-362.
Sarker, Nripendra N. and Otto Kunze. 1991. Finite Element Prediction of Grain
Temperature Changes in Storage Bins. Paper presented at the 1991 International
Winter Meeting of the Am ri i ii 1 n 'n r , St. Joseph,
MI.
Schreyer, H. L. 1981. Nonlinear ﬁnite-element heat conduction analysis with direct
implicit time integration. NW, 4, 377-391.
Scott. E. P. 1987. W
W Doctoral dissertation. Michigan
State University.
Segal, A., and Praagman, N. 1986. A fast implementation of explicit time-stepping
algorithms with the ﬁnite element method for a class of nonlinear evolution
problems. lnternatienal leumal fer Numerigl Methguls in Engingring, 23, 155-
168.
Segerlind, L. J. 1993. Personal Communications.
Segerlind. L. 1.. and Scott. E. P. 1988. W
W. Paper presented at the 1988 International
Summer Meeting of the American Society of Agricultural Engineers, St. Joseph,
MI.
Segerlind, L. J. 1987. Observations of and Recommendations for the Finite Element
Solution of an One-Dimensional Parabolic Differential Equation. Unpublished
report for Mathematics 890 Michigan State University.
Segerlind, Larry J. 1976. Anglia Fiuite Element Analysis. John Wiley & Sons, New
York.
Segerlind, Larry J. 1984. Appliﬂ Finite Element Analysis, second edition. John Wiley
& Sons, New York.
Shih, Tien-Mo. 1984. Numerigl HeatTtausfer. Hemisphere Publishing Corporation,
New York.
141
Shayya, W. H., Segerlind, L. J ., and Bralts, V. F. .1991. Optimization analysis of the
four-level-time schemes. ' '
We 11, 1113-1119°
Singh. R. P-. and Segerlind. L. J. 1974. WW2
Paper presented at the 1974 Annual Meeting of the American Society of
Agricultural Engineers, Stillwater, OK.
Smith. G. D. 1985. WWW
_DﬁfemMethms, third edition. Oxford Applied Mathematics and Computing
Science Series. Clarendon Press, Oxford.
Stoer, J. and Bulivsch. 1980. lntrﬂuetien te Numerieal Analysis. Springer-Verlag,
New York.
Strang, Gilbert. 1976. WWW. Academic Press, New
York.
Sun, N. Z. 1989. Applications of numerical methods to simulate the movement of
contaminants in groundwater. BMW, 83, 97-115.
Todini, E. , and Bossi, A. 1986. PAB (Parabolic and Backwater) a unconditionally stable
ﬂood routing scheme particularly suited for real time forecasting and control.
leumal ef Hydtaulie Reﬂmh, _2_4, 405-425.
Tsuboi, H., Tanaka, M., and Misaki, T. 1988. Computation accuracies of boundary
element method and ﬁnite element method in transient eddy current analysis.
IEEE Tmsaetiens ef Magneties, 2_4(6), 3174-3176.
Utnes, T. 1990. A ﬁnite element solution of the shallow-water wave equations. Appliﬂ
W02. 14. 20-29.
Visser, W. 1965. A Finite Element Methud fer the Qeterminatien ef Nen Statienaiy
Temmtature Distribution autl Thermal Defetmau'ens. Proceedings of Conference
on Matrix Methods in Structural Analysis, Air Force Institute of Technology,
Wright Patterson Air Force Base, Dayton Ohio.
Williams, W.E. 1980. P ' Diff r n ' ' n . 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. Nueleat Enginﬁring ahtl Design, Vol. 4,
276-286.
142
Wood, W. L., and Calver, A. 1990. Lumped versus distributed mass matrices in the
ﬁnite element solution of subsurface ﬂow. Water Reseurees Remeh, 26(5),
819-825.
Wood, W. L., and Lewis, R. W. 1975. A comparison of time marching schemes for
the transient heat conduction equation. In rn ' urn f r N m ri
Methegs in Engingring, 2, 679-689.
Wood, W. L. 1990. W. Oxford Applied Mathematics
and Computing Science Series. Clarendon Press, Oxford.
Yadava, R. R., Vinda, R. R., and Kumar, N. 1989. One-dimensional dispersion in
unsteady ﬂow in an adsorbing porous medium: An analytical solution.
Hyurelogig Prg‘esses, a, 189-196.
Yen, D. (1993). Personal communication. Department of Mathematics. Michigan State
University.
Yu, C. C., and Heinrich, J. C. 1987. Petrov-Galerkin method for multidimensional,
time-dependent, convective-diffusion equations. Intematienal leumal fer
Numeriw methﬂs in Engiuﬁring, 2_4, 2201-2215.
Zienkiewicz, O.C. 1971. The Finite Element Methed in Engineering Seieng. McGraw
Hill Book Co. London.
Zienkiewicz, O.C. 1977. The Finite Element Methm. McGraw-Hill Book Co.
London.
for: UNIV. LIBRARIES
111111111111|11111111111111111
301 219040
nrcuran 3
1111111111111
3129 0