“'3' ' MllGSATCHANT ll!!! “MN” ”Will/INNIllllllllllllml 01573 0181 L'BFLAngs “2:" l V‘chw STATE UNIVERSITY This is to certify that the dissertation entitled Element- 4“! +152 54:? W -FOr Soil/1A7 TM — Defiwdad fab! Frau»... anti—D71: M am New presented by RESQ. ”M0100” has been accepted towards fulfillment of the requirements for FIND degreein AE' r ' profe or Date WM MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771 MSU RETURNING MATERIALS: Place in book drop to LJBRARJES remove this checkout from .—_—. your record. FINES will be charged if book is returned after the date stamped below. v‘ "i“ "g‘e’iJ’éi‘; t“‘} 5 7? #77 *1 L 3*? L..- l z,~.. :3 watchman. ELEMENT AND TIME STEP CRITERIA FOR SOLVING TIME-DEPENDENT FIELD PROBLEMS USING THE FINITE ELEMENT METHOD by Reza Maadooliat A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Agricultural Engineering 1983 ABSTRACT ELEMENT AND TIME STEP CRITERIA FOR SOLVING TIME-DEPENDENT FIELD PROBLEMS USING THE FINITE ELEMENT METHOD BY Reza Maadooliat The method of estimating an optimum integration time step for irregular finite element grids was established by use of subregion analysis. The lumped and consistent for- mulations, a study of element matrices relative to the con- cepts of Positive Coefficients Rule and the determination of which elements are the most suitable for use in time dependent field type problems were also investigated. The matrix stability analysis and Dusinberre's Criteria were used to study solutions to [C]{(f)] + [K]{¢} - {F} = {0}. The matrices [A] and [P] in the finite differ- ence solution in time, {CD}1 = [Al-1[P1{¢}O + [Al-1{FL must satisfy the following criteria in order to avoid numerical oscillations and ensure physical reality. 1. [A] must have positive diagonal coefficients and negative off-diagonal coefficients. 2. [P] must be positive definite and have positive coefficients. Other conclusions of the study include: a) Quadratic elements can not satisfy the requirements placed on [A] and their use will always result in physical unrealistic values Reza Maadooliat for some of the nodes; b) the interior angles of the linear triangular element must be less than or equal to 90°; c) the aspect ratio for the bilinear rectangular element must be less than /_; d) dividing a rectangle into two tri- angular elements significantly reduces the allowable time step; e) the convection matrix must be lumped if the cri- teria on [A] is to be satisfied for all values of the con- vection coefficient. The maximum allowable integration time step for the consistent finite element formulations is small compared to the lumped formulations. The satisfaction of Positive Coefficients Rule is also very difficult if not impossible for consistent formulations. The maximum eigenvalue of a subregion consisting of the five smallest elements in a grid yields a good estimate for the integration time step. Approved Major Professor Approved Department Chairman Dedicated to: Susan, Sohila and Mehdi ii ACKNOWLEDGMENTS I would like to thank Dr. Larry J. Segerlind, my academic advisor, for his advice and encouragement and editorial help throughout the long months of research and writing. Thanks are also due to members of the Guidance Committee and invited Examiner: Dr. J.F. Steffe, Dr. D.C. Wiggert, Dr. N.L. Hills, and Mr. Marvin Church. I would also like to thank Lou Anne Alderman for her dedi- cated typing of the final manuscript. iii TABLE OF CONTENTS LIST OF TABLES . LIST OF FIGURES CHAPTER I. INTRODUCTION 1.1 NUMERICAL OSCILLATIONS . . 1.2 CONSISTENCY AND CONVERGENCE 1.3 OBJECTIVES . . . II. LITERATURE REVIEW . III. IV. .1 NNNNN N O‘U’l-l-‘UJN 2.7 FUNDAMENTAL CONSIDERATIONS . . . SOME PROPERTIES OF [k(e?|MATRICEs MODAL ANALYSIS . . . . . DIFFERENCE RECURRENCE FORMULA PROPERTY OF [A] [P] MATRIX . . . STABILITY AND OSCILLATION ANALYSIS . 2.6.1 General Approachs . . 2. 6. 2 Heuristic Stability . 2.6.3 Von Neuman's Method . 2. 6. 4 Dusinberre Concept 2.6.5 Matrix Stability STABILITY CRITERIA FOR GENERAL SCHEME. EVALUATION CRITERIA . 3.1 3.2 3.3 MATRIX STABILITY . . DUSINBERRE CONCEPT . EVALUATION CRITERIA CONSISTENT AND LUMPED FORMULATIONS 4.1 4.2 HEAT FLOW IN A ROD . 4.1.1 Consistent Formulation 4.1.2 Lumped Formulation . OPERATING REGION FOR A ONE- DIMENSIONAL UNIFORM GRIDS 4.2.1 Consistent Formulations 4.2.2 Lumped Formulations . iv Page vi vii (”\lN [-4 CHAPTER VI. VII. 4.3 OPERATING REGIONS FOR TWO- DIMENSIONAL UNIFORM GRIDS . 4.3.1 Square Grids . 4. 3. 2 Equilateral Triangular Grids : ACCEPTABLE ELEMENTS Ul U1UI UIU'IU1U1 O‘U‘I #WNH LINEAR ONE-DIMENSIONAL ELEMENT . . QUADRATIC ONE-DIMENSIONAL ELEMENT . THE LINEAR TRIANGULAR ELEMENT . TWO- DIMENSIONAL PARALLELOGRAM ELEMENT . . LINEAR QUADRILATERAL ELEMENT . . . QUADRATIC AND CUBIC QUADRILATERAL ELEMENTS . . . CONVECTION BOUNDARY CONDITIONS ESTIMATION OF A TIME STEP VALUE O‘O‘O‘O‘O‘O‘ O‘U‘l-l-‘UDNH METHOD OF ESTIMATION LINEAR UNIFORM GRIDS . . . UNIFORM TWO- DIMENSIONAL GRIDS . . NON- UNIFORM ONE- DIMENSIONAL GRID THE RADIAL ELEMENT . . . . TWO- DIMENSIONAL GRIDS . CONCLUSIONS LIST OF REFERENCES Page 55 55 58 64 64 67 75 82 87 87 90 91 91 93 95 96 101 106 108 TABLE 2.1 6.1 6.2 LIST OF TABLES [A] and [P] for the popular choices of O . . . . . One-dimensional grids used to calculate ASIAN . . . Two-dimensional grids used to calculate xS/AN . . . vi Page 23 97 107 FIGURE HI—Il—l PLAN J-‘J-‘b-L‘ .10 .11 LIST OF FIGURES Stable Convergent solution . Stable and oscillating solution Unstable solution Heat transfer analysis in an insulated rod General two-dimensional heat conduction problem Variation of A9 with respect todifor different values of O Insulated rod with heat at node 1 Nodal temperatures for the consistent formulations (Analytical Solutions) Nodal temperatures for the lumped formulations (Analytical Solutions) Stability and non-oscillation criteria (consistent) Non-oscillation region (consistent) Stable Operating region for consistent formulation . Operating regions for lumped and consistent formulations A grid Of four square elements . Operating regions for square grids . Equilateral triangular grids . Operating regions for equilateral triangular grids . . . . vii Page 13 . 34 . 40 . 43 . 45 . 49 . 51 . 52 . 54 . 56 . 59 . 60 . 63 FIGURE UUU'IU'I UlU'IU'IU'IUIU \lO‘U‘I-L‘UONl—i O‘O‘O‘O‘O‘ Ln-l-‘WNH CD .10 .11 One-dimensional quadratic element . Two-dimensional triangular element Triangular element Triangular element Variation of )‘n with respect to a Parallelogram element . Operating region for parallelogram element . . . . . . . . . . Variation of An with respect to 2'. Linear quadrilateral element Quadrilateral elements Quadratic elements with corresponding stiffness matrices Right triangular element One-dimensional non-uniform grids . Linear radial element . Unequal triangular grid . Two-dimensional non—uniform girds . viii Page 66 68 70 73 74 76 81 83 84 86 88 94 98 99 . 103 105 CHAPTER I INTRODUCTION An analytical solution of the heat conduction for mixed boundary conduction defined on irregular geometries is in general very difficult if not impossible to Obtain. ‘Many of the practical heat conduction problems have no known analytical solution. The finite difference method and the finite element method have been used by many people to solve the heat conduc- tion equation. The finite element scheme has received much recent attention because of its ease in handling irregular nodal point locations, and mixed boundary conditions. The unique aspects of finite element solution is the integral formulation and the associated approximation equation. Applying a finite element scheme to the heat conduc- tion equation or a similar equation reduces the partial dif- ferential equation (PDE) to a system of algebraic or differ- ential equations. The solution of the system of differen- tial equations is not simple and straightforward. Many solution schemes exist and there are difficulties associated with each scheme. Some of these problems are discussed in the following sections. 1.1 NUMERICAL OSCILLATIONS The first problems are those of stability and numer- ical oscillations. This aspect of the numerical study has nothing to do with the original partial differential equa- tion (PDE) but rather concerns the investigation of errors in the arithmetic Operations needed to numerically solve a system of ordinary differential equations (ODE). Although, there is an analytical solution for the system of ODE's, the system is often solved using a numerical method. In this process, the ODE's are transformed into difference equations which have an error associated with their solution. The total error is the numerical error plus the discretization error. The numerical error comes from rounding. According to Welty (1974), this error is negligible in modern com- puting machinary. The discretization error originates from.the replace- 'ment of an ordinary differential equation by a finite dif- ference equation. This error may be reduced by reducing the size Of integration time step (ITS). The following three figures illustrate some of the situations associated with stability and numerical oscil- lations when calculating the temperature distribution in an insulated rod of unit length whose ends are kept at a con- stant temperature of zero. The initial temperature distri- bution is: T=2X 0sxs1/2 and (1.1) T=2(1-X) 1/zsxs1 Figure 1.1 illustrates a typical stable and convergent solution. Figure 1.2 is stable, but oscillates about the exact solution. Figure 1.3 shows an unstable finite differ- ence scheme. The temperature values oscillate about the exact solution and become worse with increased time. The size of the integration time step is responsible for the numerical oscillation and the unstable solution. The question at hand is, "How large can the time step be and still avoid the numerical oscillation?" The answer for some finite difference methods are available, Myers (1971), Yalamanchili (1973). The difficulty arises when the region of the body does not have a nice shape and is composed of several materials or has mixed boundary con- ditions. Finite element grids are nearly always irregular and little information for estimating the allowable time step for irregular grids is available. For example, the ANSYS finite element program.manual Swanson (1981) suggest that the finite difference estimate for a uniform.grid be used to estimate the integration time step. It is also pointed out that the presence Of heat transfer mechanisms other than conduction, such as radiation, etc. may require additional guidelines for selecting ITS. For example, abrupt changes in the ITS (i.e., changes by more than a .cOHuDHOw unmwum>coo manmum mozauoo>aoo 14 2 2 1T = [1/2[Kx(g—;I;) + Ky(%—'}r; - 2(Q -pC—g%—)T] dv (2.7) + Iqus + Ila-(T-TdeS 51 32 where Tf is the ambient temperature, V is the volume of the body, and S is the surface area. Defining two matrices T __ 3T 3T {g} - [a_x a—y-] (2.8) and K .0 [D1 = x (2.9) 0 Ky (2.7) can be written as 1! = f1/2 [{g}T[D] {8} - (Q " OC%%)T] dv v (2.10) + Iqus + I%(T—Tf)2ds s1 82 The function for T is continuouslnnzdefined by several equations over the region. The individual equations, T(e), are defined over subregions called elements and the above equation should be separated into integrals over the individual elements, therefore, 15 E N = f 1/2 [{g(e)}T[D]{8(e)}]dv e=1 v‘e) (e) _ T(e)[Q(e) _ p(e)C(e)—3—E- ]dV (2.11) v(e) (e) + J[ T(e)q(e)ds + Jf %. (T(e)-Tf)2ds (e) 51 82 or E 21w e=1 where E is the total number of elements. The function u is minimized with respect to a set of nodal values contained in {T}. Minimization produces an a. E (e) E 3W(e) m =mzlfl = 1 my 5' 0 (2.13) e= e= The temperature in each element is given by (e) - (e) (e) (e) T -Ni Ti+Nj Tj+... +Nm Tm (2.14) or T(e)= [N(e)]{T(e)} (2.15) 16 where [N(e)] is a row vector of shape functions (inter- polation functions) and {T(e)} are element nodal tempera- (e) T _ tures {T } - [Ti Tj ...Tm]. Some properties of the interpolating functions include: 1. The interpolating functions must sum to one at every point in the element. Ni+Nj+ooo+Nm=1 2. The derivative of the shape functions with respect to a variable sums to zero, i.e., 3N. 3N. aN __£ _.1 _II;=_§_ , ____a_ = 3X + X +°°'+ 3x ax(N1+Nj+"°+Nm) 3x(1) 0 (2.16) The temperature gradient can be written in terms of the interpolation functions as {g(6)} = [3(3)] {T(E)} (2.17) where [3(6)] is a gradient matrix *BN- 3 . . . BNm __1 _N.'1 _ [B] = 8x 8x 3x 3y 3y ... 78; Introducing (2.14) and (2.17) into the functional and performing the minimization, gives the matrix differential equations, Segerlind (1976). 17 [c1§5{%}+ mm = {F} (2.19) where the element contributions to [C], [K] and {F} are the element capacitance matrix [0(3)] = ‘fpc[N]T[N]dv (2.20a) V The element stiffness matrix [k(e)] = fiBlTlD][B]dv + fh[N]T[N]ds (2.20b) V 82 and the element force vector {£(e)} = - JfQ[N]Tdv + Jrq[N]Tds - jrhrftules (2.20c) V 81 32 The three integrals are evaluated over each element and the element contributions are added using the direct stiffness procedure. A superscript (e) on the left of the equal sign implies that all the quantities on the right hand side are to be interpreted on an element basis. 2.2 SOME PROPERTIES OF [k(e)] AND [c(e)] MATRICES The stiffness matrix (conductivity matrix), [K], and the heat capacitance matrix, [C], have some interesting properties. l8 1. Writing (2.20b) using indical notation and deleting the surface integral, which means either no losses through boundaries or a specified boundary temperature, a coefficient in the element matrix is - BM 3.); f 3N1 flu ki,j - [KX— 3x dX+ y-fi 3y dy (2.21) V V th Taking the sum of the i row of [k(e)] 2: ki j sum of individual coefficients in ith row j'1 u x 02' NLJ 1":\ M22 M'u LO) :2 \f_,/ o. < + ‘< 0’0) “4;; z’_‘\ :z '6 0’0) ~2. 34 o .0 mo mmDHm> 6:36.33 .46: 8 6883 :33 4 mo 833.35 N.N 85w: Dine m N a a p p b p . h p p p p b L h _ _ .. one .. 1 N\HH® .l m\~u® 1 1- ,lr'llill 'I'I 'Ill'l-|l, r) 1! Huo CHAPTER III EVALUATION CRITERIA As mentioned previously, the analysis of the time dependent heat conduction equation may be regarded as con- sisting of three parts: (i) Discretization process, (ii) numerical stability, and (iii) Dusinberre's concepts. The various kinds of stability analysis and the discretization concepts were discussed in the previous chapter. In this chapter the matrix stability along with the requirements of Dusinberre's criteria are explained. 3.1 MATRIX STABILITY The non-oscillation criterion of equation (2.51) requires that the eigenvalue of [A]-1[P] matrix, A9,1x>be bounded between zero and one. The eigenvalues, A0, are all positive if the matrix [P] remains positive definite. From.Eq. (2.64) the eigenvalue, A is always less 9’ than one provided [P] remains positive definite. Therefore, the non-oscillation criterion reduces to finding the minimum value of At that makes the matrix [P] = [Cl-(1-0)At[K] singular. As a result, the analysis boils down to finding 35 36 the smallest eigenvalue of the equation [C]{X} = a[K]{X} (3.1) where a = (l -O)At. Consequently, the allowable time step should satisfy the following inequality Atsa/(l-O) ' (3.2) which is exactly the same as Eq. (2.66), (AN=l/a). 3.2 DUSINBERRE CONCEPT The positive coefficients rule, discussed earlier, considers the stability criteria from the physical point of view. The coefficients a pn’ in (2.50), which come from [A]'1[P] have to be non-negative as well as the elements in [Al-1. One way to gurantee that the positive coefficients rule is satisfied is to require negative off-diagonal terms in matrix [A] and positive coefficients in matrix [P]. This statement is based on the following analysis. Knowing that every symmetric positive definite matrix has a unique Choleski Decomposition, the matrix [A] can be decomposed into the product of a lower and an upper tri- angular matrix. [A] = [L][L]T (3.3) 37 or P ‘ n q I .1 311 a12 aln 111 0 111 121 lnl 321 6122 a2n = ¥21 122 122 162 an1 an2°°°ann Llnl 1n2"°1nn[ L0 lnn L ' - Using the general formula for Choleski factorization, i-l 2 g 111 = (an ’ Z lik k=1 (3.4) 11j= aij - X likljk Jenning, (1977). It can be shown by induction that every diagonal of [L] is positive and all off-diagonal coeffi- cients are negative. Also, the inverse matrix of [L] is a lower triangular matrix with all positive coefficients. AS a result , the [A] matrix which is equal to [vL"'1]T[L]'l has posi- tive coefficients. 3.3 EVALUATION CRITERIA The evaluation criteria used in this study are: l. The coefficients of [A-1][P] and [Al'-1 'must be positive thus [A] must have positive diagonals and negative off-diagonal coefficients and [P] must have positive coefficients. 2. The matrix [P] must be singular or positive definite. 38 CHAPTER IV CONSISTENT AND LUMPED FORMULATIONS This chapter concerned with the consistent and lumped heat capacitance matrices. The usual finite element method produces a non-diagonal heat capacity matrix often called the consistent formulation. The diagonalization of the heat capacity matrix, however, has been considered by several researchers. The first to lump [C] were Nickel and Wilson (1966). The diagonalization of [C] can be done by simply add- ing the coefficients in each row and placing the sum on the diagonal, Myers (1977). Some special numerical inte- grations also, have been suggested for making [C] diagonal, Zienkiewicz (1977). The lumped formulations reduces the computer storage requirements and the number of computations in the step-by- step solution. Users experience has shown that a larger allowable time step may be used before introducing the oscillations in the solutions. The objective of this chapter is to compare the ana- lytical and numerical solutions of the two methods and determine the operating region for the numerical solution of each procedure. 39 4.1 HEAT FLOW IN A ROD This section compares the analytical solutions for the consistent and lumped formulations using the temperature distribution in an insulated rod, Figure 4.1, as an example. Heat input occurs at the left end and the right end is insulated. The initial temperature is zero. The element matrices are given by 2 - 8 4 [k(e)] = ; and [C(eh = -2 2 4 8 The global matrices for two uniform elements are ' 2 -2 o” '8 4 0' [K] = -2 4 -2 and [C] = 4 16 4 _ o -2 2‘ b0 4 84 and the forcing term is {F}T = [5 o 01 since the heat input is a point source. 4.1.1 Consistent Formulation Substitution of [K] and [C] into (2.34) yields r- . 1p ‘ 8 4 0 2 -2 0 5 4 16 4 {i} + -2 4 -2 {T} = 0 Lo 4 8J _o -2 2] o (4.1) 40 Heat 4 Input 1t— 2cm i 2cm fl (1) 2 (2) _ 4 watts - cm Kx — CU _ watts Pc _ 12sec - cm — C° 2 A = 1cm Figure 4.1 Insulated rod with heat at node 1. 3 41 The Laplace transform method is used to obtain the exact solution to the system of ordinary differential equations in (4.1). This method gives 'CST(s) + KT(S) = F(s) or (4.2) T(S) = (65 + K)-1F(S) From which '3s+2 45-2 0 (C3 + K) = 43-2 16s+4 43-2 Lo 43-2 88+2‘ (cs+K)'1 = 4(2852+zos+1) -4(2s-1)(4s+1) (43-2)2 1 2 4(2s-1)2 -(ss+2)(4s—2) 4(2832+zos+1)_} Also, {F(s)}T = (5/5 o 01 Multiplying (CS+K)-1{F(s)} gives 42 2882+208+1 (4s+1)(s+1)s2 ' 1-2S T(s) = -——————2 (4.3) (S+1)S (25-1)2 (4s+1(s+1)s2 Taking the inverse Laplace transform from (4.3), gives the nodal temperature distributions. T1(t) = 5/48(15 + t - 12e' 25‘ - 3e't) T2(t) = 5/48(-3 + c + 3e‘t) (4.4) T3(t) = 5/48(-9 + t + 12e'°25t - 3e't) The nodal temperature distributions obtained using (4.4) are shown in Figure 4.2. The temperature history for node two and three have negative values. Since the equa- tions in (4.4) are the analytical solutions to the system of ODE's it must be concluded that the consistent formula- tion produces unrealistic results at least for the first few time intervals. 4.1.2 Lumped Formulation The use of lumped heat capacitance matrix saves con- siderable computational efforts in determining the nodal temperature {T}. In this case, [C] is 43 Amcowusaom Hmofluhamc 1/2, the system is unconditionally stable regardless of the integration time step and physical properties. The above inequalities are plotted in Figure 4.4. It can be seen from this plot the Euler's method (9 = 0) requires the smallest value for At in order the system be stable. There— fore, this scheme is the most restrictive one for solving transient field problems. Dusinberre's criteria requires negative off-diagonal coefficients in [A] and positive elements in [P]. Lemmon, (1969) states that the coefficients in (4.7) have to be positive which is not the correct interpretation of Dusinberre's criteria. Dusinberre's analysis requires 1/6 - re s 0 or r 2 l/6O (4.12) and 2/6 - r(1-O) z 0, or “3714—6)— (4.13) 49 .ADCoumflmcouv wfluoufluo cowumaawowOuCoc paw Nuaawnmum q.q ouswwh '4 (I) a I IrII 4 l I ii: A ’II’I T —-—-—-~-—-——1 50 Equations (4.11) along with (4.12) and (4.13) are plotted in Figure 4.5 to obtain the operating region for a non- oscillatory system and similarly inequalities (4.10), (4.12), and (4.13) are plotted in Figure 4.6 to show the operation region for a stable system. It is clear that Euler's method will not yield good results because it is not in the operating region for the values of 0 less than 1/6. 4.2.2 Lumped Formulation The element heat capacitance matrix, [C(e)], for linear temperature distribution and a lumped formulation 1 0 (e) Ax [C l = and the global matrix [C] for three identical element is is '1 o o o” [C] = 4255 o 2 o o o o 2 o 9 o o 1. The global conductivity matrix [K], remains unchanged and is the same as that for consistent formulations. Following the same procedure illustrated for consistent formulations, the finite difference recurrence formula for the lumped case is 51 .Aucoumwmcoov coawou CoaumaaflomOIcoz m.q muswwm a ma 6 NA R. 3% m cowumaawomouaoc oomaoa \\ GowumHHHomOIGOG mucmflowmmooocw>wwMMMa t\ uaoumwmcoo u u . o IIII i I II L-” IIII A A A A 52 .cowumHDEHom acoumwmcoo How aofiwmu wcHUMHomo manmum o.¢ ouswfih a n: a mu 7w. 9 0. >4 9 muamwowmmooo m>fiuflmom mucofloflmmmoo o>flufimom .l; 41!!!! uamumfimcoo unoumme£D\\1 ~-——-”-—-w-‘Ffi_' 53 r 1 l + 2rO -2rO -2rO 2 + 4rO -2rO {T}n+l -2rO 2 + 4rO -2rO . -2rO l + ZrO‘ F (4.14) 1-2r(1-O) 2r(1-O) ' 2r(l-O) 2-4r(l-O) 2r(l-O) mu 2r(1-O) 2-4r(l-O) 2r(1-O) _ 2r(l-O) 1 - 2r(1-O) The largest eigenvalue of the system (4.8) is A ___. 40: (4 15) N (Ax)2 ' Therefore, the non-oscillation, stability and Dusinberre's requirements are l non-oscillation requirement r‘<4TI:§) (4.16) stability requirement r<<7?I%757 (4.17) Dusiberre's requirement . rw<§TI%574 (4.18) The operating region for non-oscillatory and Dusinberre criteria are determined by the single inequality of (4.16). This operating region is shown in Figure 4.7 along with the same conditions (Dusinberre and non- oscillation requirements) for consistent formulations. It shows that consistent formulation is more restrictive rela- tive to selecting a time step than the lumped formulation. 54 .4 ME \ N31§|S\OIA 3%,)? nowumaawomoucOc pmmfiba coaumaafiomoncoc ucoumamcoo .mcowuwansnom uaoumflmcoo paw pmaESH How mGOHon wcwumuomo m.q ouswwm mquflonmooo m>HuHmoa ucmumwmaoo 55 4.3 OPERATING REGION FOR TWO DIMENSIONAL UNIFORM GRIDS A pair of uniform two-dimensional grids were analyzed using the criteria stated in Chapter III. One grid con- sisted of square elements while the other consisted of equilateral triangular elements. 4.3.1 Square Grids A typical set of elements for a uniform square grid are shown in Figure 4.8. The element matrices for each element are (e)_D [K 1’; and 4 2 1 2 [Cm] “43564 2 4 2 1 1 2 4 2 L2 1 2 4. where A is the area of the square and D is thermal conduc- tivity. The largest eigenvalue of the system (4.8), An, is A.., = “D n ECTA- (4.19) 56 v.— 31. ‘81 (1) (4) 1 3 f— #Tsi W (2) (3) L6 2] ‘73 Figure 4.8 A grid of four square elements. 57 Substituting the computed value into (2.66) gives the non- oscillation criteria 1 r (m (4.20) where r = (DAt)/pcA. Dusiberre's criteria is related to the coefficients in global matrices [A] and [P]. Using the direct stiffness procedure, the equation for the center node is cA . . . . . . . . . %F(‘T1 + 4T2 + 413 + 4T4 + 16T5 + T6 + T7 + T8 + T9) (4.21) D _. + §(-Tl-T2-T3-T4+8T5- T6-T7-T8-T9) - O The evaluation criteria derived in Chapter III, requires the positive diagonals and negative off-diagonal coeffi- cients in [A] and positive coefficients in matrix [P]. These requirements are satisfied when r and O satisfy the inequalities r > 31—9 3 for [A] (4.22) r < 1' - for [P] (4 23) The element heat capacitance matrix for the lumped formu- lations is 58 '1 o o 0" (e) pcA 0 l 0 0 [C ']= —z— 0 0 1 O L_o o o 1.. which from the eigenvalue analysis yields the inequality 1 I3é: ; for consistent [A] (4.29) Figure 4.11 shows the operating regions for equi- lateral triangle grids for consistent and lumped formula- tions. It shows that the consistent formulation is highly restrictive. The inequalities are satisfied only for values of O greater than 0.75. 63 .mpfluw Hmaswcmwuu Hummumawswo How szHwou wcfluMHomo HH.q munwfim a R. a b b p P A... \W\ .. N r. S III m o7 . 11. .- PN. 5W 4. \\ . .. \ 11... g L. CHAPTER V ACCEPTABLE ELEMENTS One aspect of this study was to investigate which elements are appropriate for a finite element discretiza- tion. Elements which fail Dusinberre's Criteria produce mixed coefficients in the general difference equation (2.41) and give results which violate physical reality. A group of the commonly used one- and two—dimensional elements were analyzed relative to Dusinberre's Criteria. The analysis was restricted to the lumped formulation. 5.1 LINEAR ONE-DIMENSIONAL ELEMENT The element matrices for the one-dimensional linear element are 1 - 1 0 (e) (e) AX [k ]=-q- ; and [C ]=-2— (5.1) The investigation of positive coefficients in Eq. (2.41) is straight forward. The off-diagonal terms in [k(exl are negative, consequently the coefficients in [A]'li:1(2.4l) are all positive. The positive coefficient rule is satis- fied provided all elements of [P] in (2.41) are positive. 64 65 Since [C(eh is diagonal and all of the off-diagonal coefficients in [p ] are positive, [P] has positive coefficients. 5.2 QUADRATIC ONE-DIMENSIONAL ELEMENT The element matrices for the one-dimensional quad- ratic element are (5.2) ' 14 -16 2] '1 o o] [k(e)] = f—x -16 32 -16 ; and [cm] = 935 o 4 o L 2 -16 14_ [p o 1: The stiffness matrix [k(e)] contain positive off-diagonal terms thus [A] contains positive off-diagonal values and Dusinberre's Criteria is not satisfied. It is impossible to change the positive coefficients k(e)]. in [ The shape functions for the configuration in Figure 5.1 are N. = (1 - §)(1 - §), 1 E = x(x - L) “1 g—(g—rm ‘5-3’ Nk _ X(§ - X) - L(€ - L) 66 .uaoEoHo ofiumu upwsv Hmcowmcmawvuoao H.m muswum The coefficient k which is zero and 1,3 = 1,3 dN. de [szrij—d" V _____ E Kx AoL m.m onswwm 75 of An for the equilateral triangle is not appreciably dif- ferent from.the values for a = 0 or a = l which corresponds to right triangles. 5.4 TWO-DIMENSIONAL PARALLELOGRAM ELEMENT The study of negative off-diagonal terms in the stiff- ness matrix of a parallelogram element, Figure 5.6, is interesting and is a starting point for a discussion of the quadrilateral and rectangular elements. The temperature filed within the element is T(x,y) = a + azx + a3y + a4xy (5.15) 1 The four constants ai’ (i = 1,2,3 and 4), are determined by nodal values T1, ..., T4. T a -+a x + a l l 21 3Y1 + “4x1y1 T2 = al-Iazxz + a3y2 + 0(4ny2 (5.16) T3 = a1-+a2x3 + a3y3 + a4x3y3 T4 = ol-l-ozx4 + a3y4 + 04x4y4 which yields a1 = T1 (5.17) , T2 ' T1 a — 76 9K1 Figure 5.6 Parallelogram element. 0‘ W 77 - 8(T4 - 121) + (:(T4 - T 3 ab 3) (5.17) T1 - T2 + T3 - T4 “4 = ab substituting the above values for a through a4 into (5.15) 1 and rearrangement produces the shape functions. 17:: 2: JLIL (5'18) from.which 2 ll 1 1 - x/a - y/b + xy/ab 2 ll 2 x/a — xy/ab Z [I 3 -cy/ab + xy/ab 2 ll 4 cy/ab - xy/ab + y/b The gradient matrix {g} is , 3T 3N1 1 .53; 4 “-5-; (e) {g} = = Z {T 1 _ BN 3: 1—1 __1. 3y _ 3y i (5.20) {g} = [131(1‘9’} 78 where (5.21) ' 1 1 “I [B] ___ (-a _+ 5%) (3 - g) (335) 45155) 1 x x x c c x 1 ['64‘35’ ('25) (3'36) (36‘36+6)] An individual coefficient of [k(eh is obtained by evaluating 8N. ON. ON. EN. = __1-. _l_l ki,j foax‘z‘fiildA+ [K M A A For example, 3N 3N. k i 2 1 2 1,1 1185—3") dA + [Ky—3y) dA and - a k1,1‘ Kx(‘33) + Ky('3'é'5 " 215 + ’33) 79 Similar computations yields element stiffness matrix as ' 2 -2 -1 1] ' 2 1 -1 -2' K (e) 3ch -2 2 1 -1 a 1 2 -2 -1 [k l 7? a + 7g‘b -1 1 1 -1 -1 -2 2 1 L_1. -1 -2 2_ L42 -1 1 2‘ (5.22) ' 2 -2 -1 1‘ :1 o 1 01 K 2 c K —- -2 2 1 -1 c o 1 o -1 +‘6zab +7115 -1 1 2 -2 1 o -1 0 L1 -1 -2 2‘ Lo -1 o 1. Eq. (5.22) shows that the elements of the stiffness matrix are functions of the material properties and the element geometric parameters. The diagonal terms of (5.22) are always positive and k2,4 is always negative. The following inequalities are required for satisfaction of Dusinberre's Criteria. k3’4 = k1,2 < 0 (5.23a) k1,3 < 0 (5.23b) k2,3 = k1,4 < 0 (5.23c) Taking 3% = 3, (5.23a) becomes .1. m + ooh-a I A O or 80 assuming Kx = Ky. The inequality holds for values of greater than /272 and B/I72—-32 (5.24a) Using a similar procedure (5.23b, c) simplify to §<1/2(3 - W) , if e < /§/2 (5.24b) and -§ /2 and any value of a and The above inequalities are plotted in Figure 5.7. he operating region is shaded. The ratio c/a for a rect- angle is zero and the aspect ratio may vary from /2/2 to f2. The determination of An for a rectangular element indicates the square element has the smallest value and, therefore, the largest allowable time step. The analysis for the rectangular elements went as follows. The element stiffness matrix and the lumped heat capacity matrix for a rectangle with sides a and b are 82 r- '5 2(6+%) «23%) 43%) ( é) 1 1 2 1 [km] 2% (-23+§) 2(3)?) (Bug) -(B+~g) 1 2 1 1 -(B+—B-) (e-g) 2(s+§) (43+?) 2 1 1 1 _(s--§) 48+?) (-ZB+-B-) 2%?) J and 1 0 0 0 l 0 0 (e) ocA [c ]= T o o 1 o 0 O 0 l L A where K is thermal conductivity and 8 is the aspect ratio. The maximum eigenvalue of [K]{X} = A[C]{X} for above matrices can be obtained for different values of 8 using the computer to evaluate).n given 8. The results have been plotted on Figure 5.8. It shows the square element has the minimum value for An and therefore, the maximum allowable time step. 5.5 LINEAR QUADRILATERAL ELEMENT A general linear quadrilateral element with a natural coordinate system is shown in Figure 5.9. The use of a natural coordinate system allows the side of the quadri- lateral to rotate relative to xy coordinate system and still satisfy the continuity requirements. The shape func- tion for this element are given by (Segerlind, 1976). .m\n ou Dooammu £DH3 GK mo cofiumwpm> w.m ouswam 83 F: a: and M <0 m v 84 .ucmEmHm Honoumawnpmnv Hmocad m.m muswflh ‘. ~~ J ~ — I"|‘ 1‘ ~ ~ ~~ 85 N1 = %(1-g)(1-n) N2 = 21;(1+£)(1-n) (5.25) .. 1 _ 1 N3 - I(1+£)(l+n) N4 - [(1-E)(l+n) The stiffness matrix must be evaluated using numerical inte- gration techniques. As a result, the investigation of negative off-diagonal in [k(e4]is not simple and straight— forward. To obtain some information about the coefficients in [k‘eil,this matrix was evaluated for some simple vari- ations from a square. The base was fixed and the upper corners were moved together or separately to see how much the element could be deformed before positive coefficients appeared in [k(eh . Figure 5.10 shows five different shapes where the quadrilateral deformed to a point where positive coef- ficients appear in the off-diagonal location of [k(eU. The coordinates of each node are shown close to the corners. Figure 5.10a shows, that the upper corner can be moved until the side is about 47 percent of its original length. Plot 5.10b shows the movement of point three in a horizontal direction to the right. The side can be increased in length and about 56 percent. Figure 5.10c indicates that node three can be moved to the left 50 percent of its original length. Figure 5.10d shows that node three can be moved along a 45 degree line (the line y=x) toward the origin until the square became a triangle. The last figure, 5.10e indicates a significant deformation can occur 86 .muCoEoHo Hmuoumawupmsc oa.m muswwm Amy Ao.0HV Ao.ov oa.ov ANm.va on Ans Ao.oav Ao.ov Ao.CHVAo.oV Aoa.mv AoH.ov AoH.o.mHv Aoa.ov Ago Ao.ofiv Aciov Am.mv AcH.ov A68 Ao.oav Ao.ov AoH.mm.NVAoa.mo.Nv 87 along the line y=x and away form the origin before positive coefficients occur in the off-diagonal values of [k(eh. 5.6 QUADRATIC AND CUBIC QUADRILATERAL ELEMENTS Numerical calculation shows that the quadratic and cubic elements, both squares and equilateral triangles, have positive off-diagonal coefficients in [k(84] and can not satisfy Dusinberre's Criteria. The need to numerically integrate these elements prevents an analytical investi- gation. Sample matrices for each type of quadratic element are given in Figure 5.11. Since the square and equilateral triangle were the optimum linear two-dimensional elements, deformed quadratic elements were not investigated. 5.7 CONVECTIVE BOUNDARY CONDITIONS The contribution of convective heat loss to the stiff- ness matrix is given by £3h[N]T[N] ds. The element stiff- ness matrix for a one-dimensional fin problem is AK 1 —1 2 l (e) _ x h L [k 1 - T [_1 1] + --§- [1 2] (5.26) where p is the perimeter and h is the convection coefficient. It is clear from.(5.26) that positive off-diagonal coefficients can occur in [K(eh for some values of h andpn This possibility should be eliminated by using the concepts of lumping. The convection contribution should be lumped to give .mmoauumfi mmmcfiflum wcfivaommouuoo Sufi» mucofioao 08283.85 :6 ouawwm 88 s x. 2: - M u x . . 36.68 8.8. 8.8. 66.6 8.8. 88. a 8.3 8 3 8.8. 8.8 8.8. 88 66.6 «6.8 8.8. 8.8. 66.68 8.8. 8.8. 8.6 66.6 «6.8 66.6m. 66.6w 68.6n. ~6.m 8.8. 66.6 8.8. 8.8. 66.8" 8.8. r88. 88 8.6 «6.8 8.8. 8.8- Aocm..m.v fi00.~nfl 0N.NQI 00.0 0H.nnl 0m.nn 0.1th 00.0 0a.”! a". I co.“ . - . C U C . - C 8.3 8.3 a a 6:6 668 S a S a 3 a 668 00.0 0N.N0I 00.nnN 0N.Nfll 00.0 0H..nh[ 0n.nn 0.1th 6n.~n. 66.6“ 6«.~6- 6n.n- 6N.«.- 66.6» 6H.~n- 6_.~n £ 06.nn 0~.~nl 00.0 0N.N0l 00..~MN 0N.N0l 00.0 0H.—.nl 6~.—6- 6~.~n 6~.~n. 66.6“ 6~.«6- 66.nnfi 6«.«6- 66.6“ 66.6 6~.~m- 66.8” 6~.~n. 66.6 6«.«.. 66.nnu 6u.~6- AHgHv AHuov g 0N.N0| 00.0% 0H.Hnl 0n.Hn 0H.~ni 00.0“ 0N.N¢O 00.an 1 -1 3 o (e) _ AKx h L [k 1 “—L‘ [a 1]+‘%'[0 3] (5°27) The off-diagonal terms remain negative and Dusinberre's Criteria is satisfied regardless of physical properties and geometrical size. The convection related matrices for all elements (in all dimensions) should be lumped to avoid positive off- diagonal values in [K(eh. CHAPTER VI ESTIMATION OF A TIME STEP VALUE An estimation of an optimum time step value is neces- sary for efficient numerical computation. The value of the integration time step for some uniform grids is known. The difficulty arises for irregular meshes with different types of materials and boundary conditions. These grids repre- sent the usual situation in finite element solutions. The maximum allowable time step can be obtained by knowing the maximum eigenvalue of the system [K]{X} = AICHX} (6.1) This computation is not a quick one and the algebraic com- putations are large for a large system of equations. Con- sequently, the idea of an estimated value for the time step has been investigated to eliminate the eigenvalue computa- tions. Myers(1977), based his estimation of time step on the Gerschgorin's theorem and derived the safe time step for a two dimensional region subdivided into the triangular elements. The idea of element computation of a maximum time step value has also been suggested by Welty (1974). 90 91 The methods of Myers and Welty are conservative. The idea of a subregion analysis for estimating the time step value is given in this chapter. 6.1 METHOD OF ESTIMATION A method for estimating the maximum eigenvalue is based on the eigenvalue bound theorem for (6.1) developed by Fried (1979). In this regard, the lowest (lst) and the highest (Nth) eigenvalue of the global system (6.1), are denoted by A1 and AN, respectively. The bound theorem states e 0min (over all elements, e) soo .onx um woumasmcH .u x x x .w x x x .w x x x .o Aavx um woumasmaa .onx um mouaom ummm .m x x x .m x x x .q x x x .m x x x .N x x x .a paw sumo um poamwoomm unnumwogaow .< ouwmogaoo mamawm Hmfipwm amwwmuumo anomaabuaoz Buomwas Emummm mcowuwpaou hnmvcdom mawauoumz oumafipuooo pwuo .z«\m< mamasoamo ou pom: mkuw Hmcoamcoawpnoco .H.6 magma 98 .mpauw EHoMHGSIGoc HmSOfimcmEHp-mco N.o mudwflh mBZMZmAm ho mmmZDz S m .6 N a — A 6 — 6 fl - I I 6 6 6 6 6 m a m .1- I Old (U OIIVH HOTVANHDIH 99 Figure 6.3 Linear radial element. 100 The element matrices for a radial element are, Segerlind (1976). 2'IIKr 1 -1 (e) _ and HI f’+ R. [C(e)] = IQCAR 1 (6.9) — + R ”I H for the consistent formulation and 2? + R1 0 (6.10) [C(eh =.12%A§ o 2—+R. r J for the lumped formulation where f = (R1 + Rj)/2' The thermal properties Kr’ p and c are assumed constant and AR is the length of the element (Rj - R1). The estimation of allowable time step requires the computation of the largest eigenvalue of (6.1). Using the element matrices for the radial element gives E + R. ? ZUEK l - r[‘ ]{X} = AM 1 (6.11) - 1 AR 1 _ r r + Rj from which the eigenvalues are and 101 36Krf 2 2 = pc(2r2 + Ri Rj)(AR)2 For large values of R1 and small element length, R, the approximate value of A2 is 12K ~ r ~ 12a A2 = m2 ‘—(Ax)2 ‘6'”) which is the same as obtained for the linear case in car- tesian coordinates. The analysis for the lumped cases is similar and the largest eigenvalue is (6.13) which is the value for lumped formulation Dione-dimensional cartesian linear element. For small values of Ri’ the elements have eigenvalues different from those in (6.12) and (6.13) and the problem should be treated as a non-uniform grid. 6.6 TWO-DIMENSIONAL GRIDS The analysis of some two-dimensional grids using the subregion idea are presented in this section. Eight dif- ferent triangular grids with different element size and shapes and different boundary conditions were used for this study. One of the non-uniform grids with insulated bound- aries is shown in Figure 6.4. Since the smallest elements 102 have maximum eigenvalues, the subregion in the lower left corner were combined to form subregions and AS was cal- culated for each subregion. Table 6.2 contains a partial description of all eight grids and the corresponding eigenvalue ratio with respect to element number has been plotted in Figure 6.5. The results are similar to those for the one-dimensional grids. A subregion consisting of the five smallest elements gives a satisfactory estimate of the largest eigenvalue. For five elements, AS/AN z 1.05. 103 .pnflw nmaowcmfluu Hmsvocz 6.0 muswwm A 1).? 104 mnmaaum N ><><><><><>< ><><><><>< .w .h mcowuwvaoo humvaaom o>wuoo>¢ou .m opoc umuam um monaom umom paw mowumpqaon woumaomaH .¢ wauoumz oncflm Bhomaa51aoz enemas: kuw maofiuwpcoo NHMpaaom ZK\m& UUNHfioHMO OH U099 mUHHw HMCOHmCQBHUIObK—n. "~.6 66666 105 .mpHHw Enomwcducoc Hmaowmcoawpuo3e m.o ondwwm mHZMSMAm mo mmmZDz am mu 3 m s O VII I OIIVH HDTVANHDIH CHAPTER VII CONCLUSIONS The objective of this study was to investigate the numerical difficulties associated with the solution of transient field problems. This study consisted of four parts; establishment of the necessary properties for [A] and [P] in [A] {T}n+1 = [P]{T}n, investigation of the con— sistent and lumped formulations, a study of individual elements and an estimation of the maximum time step which can be used and still avoid numerical oscillations. Specific conclusions from this study include: 1. The following properties are required for matrices [A] and [P] in order to meet Positive Coefficient Rules: (1) [A] must have positive diagonals and negative off- diagonal entries, (ii) [P] must be positive definite with positive entries. 2. The maximum allowable integration time step for the consistent finite element formulation is small compared to the lumped formulations. For example, in one dimensional (linear) uniform grids the time step for the lumped formula- tion is three times that for the consistent formulation. 3. It is very difficult to meet Positive Coefficient Criteria using the consistent formulation. 106 107 4. The positive Coefficient Criteria requires the convection part of the stiffness matrix to be lumped. 5. The analytical solution of the temperature dis- tribution in an insulated rod showed that the lumped form— ulations gave results which are physically realistic while the consistent formulation gave unrealistic values for small values of time. 6. The interior angles of triangular elements should be less than or equal to 90°. 7. The aspect ratio for rectangular elements should not exceed by f2. 8. The one- and two-dimensional quadratic elements fail to meet the Positive Coefficient Rule. 9. Division of a square into two right triangles decreases the allowable time step value by a factor of 1.5. 10. The maximum eigenvalue of a subregion consisting of the five smallest elements in a grid yields a good esti- mate for the time step. 11. The uniform radial elements (equal AR) can be treated as a uniform grid. LI ST OF REFERENCES LIST OF REFERENCES Dusinberre, G.M., ”Heat Transfer Calculations by Finite Differences," Second Edition, International Textbook Co., Scranton, Pa. 1961. Fried, 1., "Numerical Solution of the Differential Equa- tions," Academic Press, New York, 1979. Jennings, A., "Matrix Computation for Engineers & Scien- tist," John Wiley and Sons, Inc. 1977. Lapidus, L. and Pinder, G.F. "Numerical Solution of the Partial Differential Equations in Science and Engin- eering," John Wiley & Sons, Inc. 1982. Leech, J.W., "Stability of a Finite Difference Method for Solving Matrix Equations," AIAA Journal,;mn 2172-2173, 1965. Lemmon, E.C., and Heeaton, H.S., "Accuracy, Stability, and Oscillation Characteristics of Finite Element Method for Solving Heat Conduction Equation," ASME Paper No. 69-WA/HT-35, 1969. Myers, G.E., "Analytical Methods in Conduction Heat Trans- fer," McGraw-Hill Book Company, New York, 1971. Myers, G.E., "Numerically-Induced Oscillation and Stability Characteristics of Finite Element Solutions to Two- Dimensional Heat Conduction Transients," Report 43, Engineering Experiment Station, University of Wisconsin- Madison, 1977. Myers, G.E., "The Critical Time Step for Finite Element Solutions to Two-Dimensional Heat-Conduction Transients," ASME Paper No. 77-WA/HT 19. O'Brien, G.G‘. Hyman, M.A., and Kaplan, S. , "A Study of the Numerical Solution of Partial Differential Equations," Journal of Mathematical Pcs, 191, pp. 223-251. Patankar, S.V., "Numerical Heat Transfer and Fluid Flow," Hemisphere Publishing Corporation, McGraw-Hill Book Co., New York, 1980. 108 109 Segerlind, L.J., "Applied Finite Element Analysis," John Wiley & Sons, Inc., New York, 1976. Trent, D.S. and Welty, J.R., "A Summary of Numerical Methods for Solving Transient Heat Conduction Problems,” Bulletin No. 49, Engineering Experiment Station. Oregon State University, Corvallis, Oregon, 1974. Trujillo, D.M., ”An Unconditionally Stable, Explicit Algorithm for Finite Element Heat Conduction Analysis," Nucl. Eng. Des. 41(1977)l75-180. Wilson, E.L., Nickell, R.E., "Application of the Finite Element Method to Heat Conduction Analysis,” Nuclear Engineering and Design, Vol. 4, 1966, pp. 275-286. Yalamanchili, R.V.S. and Chu, S.C., "Stability and Oscil- lation Characteristics of Finite-Element, Finite- Difference, and Weighted-Residual Methods for Tran- sient Two-Dimensional Heat Conduction in Solids," Journal of Heat Transfer, Trans. ASME, Vol. 95, May 1973, pp. 235-239. Yalamanchili, R., "Accuracy, Stability, and Oscillation Characteristics of Transient Two-Dimensional Heat Conduction," ASME paper No. 75-WA/HT-85. Zienkiewicz, O.C., "The Finite Element Method," McGraw- Hill, London, 1977. ICHIGAN STATE UNIV. LIBRARIES I“WWIIHIIWWNWW“WWI“ 31293015730181