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 fulfillment 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 field 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 field 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 finite element, and three finite difference schemes for two- dimensional problems. Comparisons between finite element and finite 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 finite element and finite 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 flexibility 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 defined 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 azalfituuan. éED‘DJHJaémJS‘WD.‘ Jagchlluflldamlul éflkéé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 financial support during the first 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 financial 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 fills 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) Justification 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 finite element schemes. page 84 Computed time step variations with the number of nodes for the three finite 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 fitting the GCD time step estimate equation and the Fr. time step estimate for a fixed boundary condition and material properties i.e. fixed 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 specific 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 definition. 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 fields are governed by a class of equations called Parabolic Equations. In the engineering literature they are often referred to as Difl‘icsion Equations. These equations have the general form: 6U _ = V. (1-1) cm k (VU) where c is the capacitance coefficient, k is the conductivity coefficient and U is the unknown variable: temperature, moisture content, pressure head, etc. Equation 1-1 applies to transient heat conduction in solids, gas diffusion, flow of fluids 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, finite elements (FE) or finite 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 finite 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 briefly introduce some of the processes that are governed by the diffusion convection field (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 flush 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 flow 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 specification. 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 flow 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 (justification) 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 finite element and finite difference literature. The solution of the time dependent field problems using the finite element method was discussed only briefly in early finite element books. Heubner (1975), discusses the derivation of the capacitance matrix [C], for transient heat transfer but never discusses the solution of the resulting system of 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 significant. Until recently, most application oriented books in heat transfer and groundwater flow centered their discussion of numerical solutions on finite difference methods. Base and length of presentations are in favor of the later method. Jaluria and Torrence (1986) discuss the three node triangular finite element for solving heat transfer problems but do not make any comparisons with a two dimensional finite difference solution. These authors do not discuss any of the other types of two-dimensional elements. Their discussion could lead one into thinking that the three node 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 finite 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 finite 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 finite difference and finite 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 define 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 finite element software when solving nonlinear problems. This method suffers from the need to repeat the calculations if the time step is too large and also does not give any information on how to select a starting value for At. Reddy ( 1984) has a section on time dependent 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 defined 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 defined by smith does not include material properties since cp/k term was defined as one. Allaire (1985) called Smith’s r term the Courant number. Allaire’s variable included the material properties. In addition to illustrating stable and non-stable schemes, Allaire defined an oscillatory stable scheme as having spatial oscillation that eventually dies out with the solution converging to the correct steady state values. Allaire showed the following criteria to be true for the single step methods: 0 < r S 0.25 No oscillation 0.25 < r 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) defined 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 defined 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 satisfied by any scheme in order to be mathematically sound. Briefly, 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) defined 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 define the numerical value for the time step that will minimize these errors. Rushton and Tomlinson (1971) used the alternating direction approach as a numerical scheme. They studied stability and found that for different boundary conditions the Courant number, C, that 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 satisfied stability requirements. Williams used a term equivalent to the Courant number and stated that it should be less than 0.5. Fried used the stability criteria (At 5 2/maximum eigenvalue). Haghighi and Segerlind (1988) solved the coupled heat and mass transfer equations using the finite 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 finite element solution for temperature distribution in a storage bin. They used the Crank-Nickolson scheme for the time domain. They presented comparisons between numerical and exact solutions. There was no mention of time step in their paper Irudayaraj (1991) and Irudayaraj et. a1. (1990) applied the finite element method to the solution of a coupled heat and mass transfer problem. Both papers used the stability criteria for selecting the time step. There was no check whether this time step ensured accurate results. The 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 modified Runga Kutta method to solve the parabolic system. Their work did not discuss solution accuracy. Peraire et. a1. (1988) studied finite element solution of fluid flow. 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 fluid, and K is a constant. Alagusundaram et. a1. (1991) applied the finite 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 finite 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 finite difference methods. They ensured accuracy by reducing the time step until the numerical results converge to a consistent value. They encountered a stability problem and a physical reality violation that they 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 finite 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 first to deal with the time step estimates from the accuracy perspective. They presented a time step estimate for one- and two-dimensional problems that 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 flow and advection-dispersion equation. He applied a variation of the linear finite element method and compared his solution with analytical ones. No indication was given as to what time step was used in his analysis. Scientists reporting new time stepping schemes seem to discuss stability, and oscillations only. 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 define 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 sufficient to ensure stability and accuracy. Although this estimate might be a good starting point for some problems, there is no justification for its use nor any experimental work that back this claim was given. Maadoliat (1983) studied stability and physical reality oscillations of the finite 15 element numerical solution. He concluded with a set of conditions that must be satisfied in order to avoid both numerical problems and recommends a time step estimate accordingly. He did not consider the accuracy criteria. 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 quantifies C and its variation with time step this relation is useless. A complete analytical approach to finding 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 sufficient 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: finite element or finite 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 finite 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 findings are the basis for redefining 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; amplification 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. - finite difference time discretization is superior to the Galean finite element time discretization. The following restrictions were incorporated in this study: - use finite 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 satisfies 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 first 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 finite element and finite 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 field problems. CHAPTER FOUR ONE-DIMENSIONAL FIELD PROBLEMS Parabolic differential equations governs a variety of time-dependent problems in science and engineering. One-dimensional field 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, flow of fluids and transport of solutes in porous media, to name a few. U is the temperature in a heat conduction problem, pressure in a flow 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 finite element (FE) or the finite 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 definite symmetric matrices. [K] is singular before boundary conditions are imposed. The matrices [C] and [K] are built from element contributions whose coefficients 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 coefficient 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 Stifi‘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 definite all the eigenvalues are real, positive and distinct and the associated eigenvectors are also real and linearly independent. If we define [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. Define 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 final solution is a linear combination of the Z,(t) solutions 24 (1,0) = Ele(t) + 15,2220) +... + Efinzno) (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, finite element or finite difference methods can be used to solve (4-2) in the time domain. Although the finite element method shows clear advantages over the finite 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 confined to the single step schemes defined 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 finite difference method, while 0 = 0.6667 is associated with the finite 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 first 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 satisfies an accuracy criteria. 2. Develop recommendations on the solution of the parabolic partial differential equation based on the first 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, defined 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 definitions 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 defining a time step estimate. The infinity 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 defined 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 definition in (4-16) comes from e’“, the first 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 finite element program was modified 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 specification 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 specification 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 defined earlier in (4-15) and used in this investigation is the ratio of these Ll norms. Each case, defined earlier in this section, produced one error ratio value. Varying N and time step for a certain case creates a set of points that fits 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 figure 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 fitting 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 find 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 figure 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 coefficients. Similarly, the backward difference and Galerkin schemes run parallel but with different coefficients. This behavior may or may not have any significance. 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 .. fi -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 defined in (4-15). These estimate define 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 coefficients of the analytical solution except the first 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(sinfl)(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 coefficient 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 finite 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 ’“fl-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 first 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 defined 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 five problems defined 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 define 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 defined 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 figures 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 figure 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 confirm 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 fixed and the time step was optimized. Another approach to accuracy is to select the proper mesh that produce accurate results for a fixed 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 defined 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 defined 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 figure 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 refined. 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. /fi _.-;———— 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 defined 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 fixed 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 fit was performed for each scheme and the Optimal Froude number was determined for a fixed 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 defined 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 confirm 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 infinite series of such curves the figure Shows important features and characteristic of the comparison between the two time step estimate approaches . The figure shows that the central difference scheme is superior to the other schemes. The difference between central difference scheme and the other schemes is significant. The other important feature is that the Froude number equivalence is not a fixed parameter and is not sufficient 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 significant. 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 significantly 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 define 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 defined 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 figure 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 finer 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 field equation have the following general form: 2 2 Dfl+Dfl+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, flow of fluids and transport of solutes in porous media, to name a few. U is the temperature in a heat conduction problem, pressure in a flow 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 finite element (FE) or finite 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 definite 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 coefficients 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, finite element or finite difference methods can be used to solve (5-2) in the time domain. Although the finite element method shows clear advantages over the finite 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 finite element and finite 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 finite element in space, forward finite difference in time. 2. Galerkin Central Difference GCD: Galerkin finite element in space, central finite difference in time. 3. Galerkin Backward Difference GBD: Galerkin finite element in space, backward finite difference in time. 4. Forward Finite Difference FFD: Explicit finite difference in space, forward finite difference in time. 59 5. Central Finite Difference CFD: Explicit finite difference in space, central finite difference in time. 6. Backward Finite Difference BFD: Explicit finite difference in space, backward finite difference in time. The Objectives of the Two-Dimensional Study are: 1. Develop a time-step estimate for six numerical schemes given above that satisfies an accuracy and stability criteria. 2. Develop recommendations to the solution scheme for the parabolic partial differential equation based on the first 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, finite element and finite difference. METHODOLOGY The methodology is briefly 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 defined accuracy limits but too small for efficient 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 verified 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 finite element and the finite difference methods. Each method uses a different approximation or integration scheme to develop the element matrices. The finite 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 finite difference on the other hand approximates the slope of the function between two nodes using the Taylor’s expansion series. In finite 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 finite element and finite 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 finite element books, Segerlind (1984), Reddy (1984). The derivation of the finite difference element matrix is not a common presentation in the literature. The element capacitance matrix for the finite 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 finite 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 finite 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 finite 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 defined 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 definition (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 difficult 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 efficient. The comparisons between the numerical ODE and analytical PDE solutions were made using an error ratio defined 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 definitions 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 difficult 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 sufficient for other problems when the grid and material properties are the same and the time to steady state increases. The step change problem is defined 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 defined flu 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 definition 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 first 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 finite element and finite 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 defined as 4 t =_ 5-12 ,, A] r 1 where A, is the lowest eigenvalue for the system. This definition comes from the analytical solution of the PDE. At steady state only the first 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 defined earlier was solved for six different grids. Each grid was identified 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 finite difference and finite 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 significant 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 figure shows the response line of the error as a function of time step for the finite 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 figures 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 figure were multiplied by A, which is evaluated to be 4.75 for the step change problem. The data were then best fitted 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 .oEofim 000 .595 mofim 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 [fir ”.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 insignificant 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 Dzfierence At = Vii/€707 N 2 25 (5-13) Galerkin Central Dtfierence At = NIT?” 2 25 (5-14) Galerkin Backward Dzfierence At = xiii—501 N 2 25 (5-15) Forward Finite Dzfierence At = _l'_19_ N 2 9 (5-16) xlNlm 79 Central Finite Dzfierence At = _19_ N 2 9 (5-17) *1 No.55 Backward Finite Dzfierecne 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 finite difference method showed some superiority over the finite element method for coarse grids as well as for small At. This feature did not carry on for finer 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 finite difference method showed acceptable accuracy for 9 node grids while finite element method rejected this grid as being inaccurate. In general the two schemes were identical and this is confirmed 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 finite 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 finite element and the finite 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 finite element and the finite 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 fine 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 finite element and finite 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 finite element and the finite 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 finite 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 refined (larger N) the time step required to maintain accuracy within a fixed 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 finite 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 defined 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 first 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 first six terms. This is particularly true at large values of time. Nevertheless, the 25 terms solution improved fitness 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 first 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 specific 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 defined 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 finite element and finite difference schemes are identical and are discussed together. The graphs in these figures 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 specific 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 specification. 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 finite 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 figures 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 finite element or finite 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 field 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 modified for use with a standard finite element program. The largest eigenvalue of the generalized eigenproblem of (6—1) is needed for the determination of stability and whether amplification factor oscillations will occur. An excellent method for estimating hm exists when using the finite 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 difficulties. 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“’] = fi 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 defining the research objectives is presented in this section. A FORTRAN code computer program by Bathe and Wilson was modified 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 floating 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 fixed 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 definition 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 figure 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 coefficient 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 significant 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 fit 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 floating and needs to be defined 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 floating 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 finding the optimal time step that accurately integrated the problem. These points were fitted by a power function and a regression equation determined the unlarown coefficients. 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 fit 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 findings 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 Dtfierence At = 1'8 N 2 25 x1 N104 Galerkin Central Dtfierence At = 1'6 N 2 25 x1 No.55 Galerkin Backward Dzflerence At = T1275“ 2 25 1 Forward Finite Difference At = L?— > 9 x1 N101 Central Finite Dtfi'erence At = _1_°6_ 9 x1 No.55 . . . _ 0.05 Backward Finite Dlflerecne 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 fixed number of elements, the central difference schemes is less sensitive to increasing time step than the other schemes. 131 For a fixed 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 finer 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 finite difference method increases the efficiency of the method computations and allowed using the same code that is used in the finite element analysis. This gives the finite difference method the same advantages that have been attributed to the finite 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, finite elements and finite differences gave time estimates that are very close to each other suggesting a similar accuracy. The finite difference method was more accurate for a coarse grid. This slight superiority did not carry for finer meshes and was not significant 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 significant 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 fixed 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 finite element or finite 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 find 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. Juumfl gf Fug! Prmess Enginfiring. 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 Enginfirs paper 88-WA/CIE-4. Boyce E.W. and RC. Diprima. 1986. Elementary Diffsrenufl 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 Numgrifl 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 Methfls 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 flow equation. W r R r R h, 25(7), 1483-1496. Cheney, W. and D. Kincaid. (1985). Numsrifl Mamsmau'us aud Campuu'ng. Second Edition. Brooks/Cole Pub. Co. Pacific 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 flow Problems. MW, 12(5). 1271-1285. Dhatt, Gouri and Gilbert Touzot. 1984. The Finite Element Methetl Displayfi. John Wiley and Sons. De Baerdemaeker, J ., Singh, R. P., and Segerlind, L. J. 1977. Modelling heat transfer in foods using the finite-element method. leumal ef Feed fluxess 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 Vflue 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 finite 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 Methtfls 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 Numerifl Methguls in Engg., Vol. 8, 633-647. Lynch, D. R. 1984. Mass conservation in finite 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 finite 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 Methfls. 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 finite 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 flow. 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 finite-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 finite 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. Applifl 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 _DfifemMethms, third edition. Oxford Applied Mathematics and Computing Science Series. Clarendon Press, Oxford. Stoer, J. and Bulivsch. 1980. lntrfluetien 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 flood routing scheme particularly suited for real time forecasting and control. leumal ef Hydtaulie Reflmh, _2_4, 405-425. Tsuboi, H., Tanaka, M., and Misaki, T. 1988. Computation accuracies of boundary element method and finite element method in transient eddy current analysis. IEEE Tmsaetiens ef Magneties, 2_4(6), 3174-3176. Utnes, T. 1990. A finite element solution of the shallow-water wave equations. Applifl 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 Enginfiring ahtl Design, Vol. 4, 276-286. 142 Wood, W. L., and Calver, A. 1990. Lumped versus distributed mass matrices in the finite element solution of subsurface flow. 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 flow 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 methfls in Engiufiring, 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