,, ..,I.1I....I., 1.51.11... , ..‘.,.I... .. , ,. -. .'. .1171 1 1' I..1>1'.',.,, 1 .1 . 1 1'..; ’. 1'.':I,'. ," "1618.13.13; 1"“ "L" . .,. 1, . .'. .’."..‘I.. My , .- 713:: .1I11,',:'1.'-I::;" ' I5. II II” (In! l‘4<,u..:.u '- ,117 ! H 3:! {1,111,310 1: (:1, '1 . ':..,. ”n.1,... v 1, up, n i 1 l ".1... ..... .- 1 r", .4. u. rIun- :- . 1.“ 'r i W" 141.. ..,i..,.1 . . {I} :lw 1-4 . 1‘;n-h'A-‘n “1.11.11.27.18 '."';‘ ..1..1.I1.... I . xl'v , hm"... \- ""n" m, 1 1,117.1»: "‘1 1. . div-'7 In .. .".1',1,.,1'1": I1 1.:tz;',‘.2,l1.1 1 fit“ ".3. .- 1.4. ....,I:I1.u;¢y-.. ..ym.1,11...\1m .1 "1'... 11.15.53 IIIII 0 '- 1'99". ‘11:..11‘g;.1.' 131317113 . 1wa. w ' 14 fr 1 m I {"1}: 'r a1§§§§gyisf 21% . . 1M" ‘V‘MIV‘M almmfirh éd” . u 1'71 n1» -- .mI 1‘ ‘1" EEK/:1, .I' 117;"? ’ I - .1 . \nm mar 1 ‘Nfrégz'u Ti: L"" " ‘u, ‘L 12.... . g h. :Ir-Jw “gs-({- , 1.. ”‘7 , 1m. 1 1 ”1" "I! , ,...‘11:<.-’ , A}. le-rtflv'J .11.... ".52.-.. 1.91:4. . «I1: rII 1-4 \ :Iau‘iq‘f.‘ .1 511‘ ‘1.“ ill-".13" -1I'v'1~v'n‘- n~ 'l {1- 1 1 1-1:: «I .' - 1 11.1. 11.1.. {1;51'2‘r': 1 "'1" . "1. .n ins? I" ' ’:.,” ,1. ”NEW?"""";‘,.,1..1'L-1 .’"' . 5". 1. 1 u ,1:- 1 ':‘JLiII'aC‘w ’ '1 n.- 1 1.1-1. ., 14.1 -.-1 “at“... 1: . ”III. :;'I 1 "1:11: .1": 1... z -...,. 1.1. . .1 :u‘ ' .- 1 n 1 Hr 1" n 5;“ I". 1 n . 1‘ '4 ;‘ 1'. .I.';2111.aI.<,‘ 11,.l'.'.‘1:.1r.'..."1..,1:_. ,. 1- , .1 . t 1.» $1,. 5 1 -- - ... cignm .~-n- "1:111:11. ”n ".3. m 21.“, 1 . .1Mm.1.."”'1L1I..,.I.L .1 an > -- -, '1'1 6'1: ~171ng v :5! '1' . 1, 1 1 1 -n. 3 v , 1‘15. 1 mun.“ . 1-1' '1--~4~ Mr- “.1 1 5'" .:.1.':!‘,.1.§1.t.1.... .... .., ..;,,.1.1..I. vulv- v} ”I :’i" u. 1 F" «11.1 4 Irv. 1151. aw raw: ,1. w- “van-On mu "1" "2!."1‘" . .1'.: w "' ".1‘,'..:."1'T. .1“ " ,« ..I.I... 3..".' x '1'." 111"?" I;.I:1 "1'21” *1 .- ;:1. ‘- 1;£, W‘ WITH; .1. .1“ :5. 31,1: 1 341,. «1""4131'95- H. (“313, “11.53 v. '1’”: I... :{1' "Item's-s CINHGA IIIIIIIII: 300897 5041 This is to certify that the thesis entitled -ON ACCURACY ESTIMATES FOR FINITE ELEMENT EIGENVALUE COMPUTATIONS presented by Yun—Jae Kim has been accepted towards fulfillment of the requirements for M.S. Mechanical Engineering degree in IIIIII/ MIME Major professor Date 5.35040 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution C LIBRARY Michlgan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATEDUE DATE DUE DATE DUE MSU Is An Affirmative Action/Equal Opportunity Institution c:\cIrc\dateduo.pm3-p.t 0N ACCURACY ESTIMATES FOR FINITE ELEMENT EIGENVALUE COHPUTATIONS BY Yun-Jae Kim A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 1990 (04(0" 3506' ABSTRACT 0N ACCURACY ESTIMATES FOR FINITE ELEMENT EIGENVALUE COMPUTATIONS By Yun-Jae Kim Accuracy estimates for finite element computations of elliptic eigenvalue problems are presented. From results of error estimates for elliptic boundary value problems, error estimates for elliptic eigenvalue problems are established for each mode based on interpolation error theory. An error indicator, which is defined as an element-wise approximation to the true error, is computationally implemented using two different techniques, a smoothing technique and a direct substitution technique using a differential operator. The accuracy of error indicators is checked using the eigenvalue problems of a uniform, Bernoulli—Euler beam. Two measures of accuracy are discussed. The first measures how accurately the estimator can capture the maximum true element-wise error over all elements. The second measures how accurately the estimator estimates the distribution of true error over the domain. Using the simple node-moving algorithm, improved grids are constructed keeping the total number of degrees of freedom constant. ACKNOWLEDGEMENTS The author would like to express his deepest gratitude to his advisor, Prof. A. Diaz, for his valuable guidence, encouragement during his years at MSU. The author also feels highly indebted and would like to express his most sincere gratitude to his family, especially to his father, Tae-sup Kim, and to his mother, Kun-ja Hong Kim. The author would like to thank other Korean students in ME Dept. at MSU, Keun-Chul Lee, Moon-suk Seo, Seoung—chool Yoo, Seung-bok Choi, Jung-k1 Lee and Sang-woo Kim, for their help and friendship. Finally, the author wish to give special thanks to Ji-young Sung and.his wife, Hae-jung Park Sung, for their support and friendship during these years at MSU. iii TABLE OF CONTENTS 1 Introduction 2 Review of Adaptive Finite Element Methods 2.1 Adaptive Methodology 2.2 Error Measurement 2.2.1 Residual Error Estimates 2.2.2 Error Estimates based on Interpolation Error 2.3 Adaptation Strategies 2.3.1 Elliptic Problems 2.3.2 Eigenvalue Problems 2.3.3 Time Dependent (Parabolic and Hyperbolic) Problems 3. Eigenvalue Problems and Finite Elements 3.1 Solutions in Infinite Dimensional Space 3.1.1 Differential Form 3.1.2 Weak (Variational) Form 3.2 Approximations in Finite Dimensional, Finite Element Subspaces 3.3 Summary A Error Indicators for Elliptic Eigenvalue Problems 4.1 Review of the Error Estimates for Elliptic Boundary Value Problems 4.2 Error Estimates for Elliptic Eigenvalue Problems iv 10 ll l6 17 20 20 20 23 27 31 32 32 3S 4.3 Error Estimators 4.4 Computational Implementation of Error Indicators 4.4.1 Smoothing Techniques 4.4.2 Direct Substitution Technique using a Differential Operator 4.5 Accuracy of Error Indicators 4.5.1 Accuracy on Maximum Errors 4.5.2 Accuracy on Error Distributions 4.5.3 Node Relocation : An Improved Grid 5 Conclusions and Recommendations 5.1 Conclusions 5.2 Discussions and Future Research List of References Appedices A Mathematical Background A.l Real Linear Space A.2 The Continuity Class Cm(0) A.3 The L2(0) class A.4 The Sobolev Class Hm(0) A.5 Equivalence of Two Norms B Interpolation Functions for Quintic Beam Elements B.1 Quintic Element Type 1 (3 nodes. 2 dof/node) B.2 Quintic Element Type 2 (2 nodes. 3 dof/node) C Exact Eigenpairs of a Uniform, Bernoulli-Euler Beam 37 39 39 42 44 46 61 74 88 88 89 92 96 96 97 98 98 98 100 100 101 103 C.1 Pinned-Pinned Beam (SS-SS) C.2 Clamped-Free Beam (CL) C.3 Clamped-Pinned Beam (CL-SS) C.4 Clamped-Clamped Beam (CL-CL) C.5 Normalization vi 103 104 104 104 106 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure LIST OF FIGURES Approximation error in elliptic boundary value problems (measured in the energy norm) Approximation error in elliptic eigenvalue problems (measured in the energy norm) ...... Relative percentage error in eigenvalue vs. Number of cubic beam elements Error in eigenfunction vs. Number of cubic beam elements (pinned-pinned boundary condition) ...... Error in eigenfunction vs. Number of cubic beam elements (clamped-free boundary condition) ...... Error in eigenfunction vs. Number of cubic beam elements (clamped-pinned boundary condition) ..... Error in eigenfunction vs. Number of cubic beam elements (clamped-clamped boundary condition) .... Relative percentage error in eignevalue vs. Number of quintic beam elements ........... Error in eigenfunction vs. Number of quintic beam elements (pinned-pinned boundary condition) u ..... Error in eigenfunction vs. Number of quintic beam elements (clamped-free boundary condition) ...... Error in eigenfunction vs. Number of quintic beam elements (clamped-pinned boundary condition) ..... vii 33 36 A 5’ ‘0'- o .53 .56 .59 Figure Figure Figure Figure Figure Figure 10. Error in eigenfunction vs. Number of quintic beam elements (clamped-clamped boundary condition) ... Percentage error vs. Location along the beam I (cubic beam element) ................. Percentage error vs. Location along the beam (quintic beam element) ............... Error in eigenvalue and eigenfunction vs. Nodal position .......................... Error distribution of uniform and improved grid . Nodal location of the improved grid ............. viii ..60 78-79 .81-85 CHAPTER I. INTRODUCTION The finite element method has become an effective and powerful tool for obtaining approximate solutions to engineering problems. Much of the power of the method is due to the freedom that it allows in the construction of the discretized model. However, the quality of the finite element approximation greatly depends on how the discretization is performed. For this reason, in the last ten years considerable effort has been devoted to adaptive finite element methods which are designed to automatically improve the quality of the finite element approximation. The goal of adaptive finite element methods can be achieved only after computationally useful measures of the ’quality’ of the approximation are available. This thesis discusses the construction of such measures for the simple class of 1- D eigenvalue problems. In finite element methods, two basic techniques for error estimation have emerged : one based on residual error estimates, introduced by Babuska and Rheinbolt [1-4] and another based on interpolation error, introduced by Diaz et.al. [14-16]. The first . . h class of error estimates is constructed from estimates of £(u-u ), whered‘e denotes a differential operator and u and uh are the solution sought and its approximation. The second technique is based on interpolation error theory and is constructed from an estimation of higher order derivatives of u. Having resolved the issue of error estimation, one can improve the quality of finite element approximations efficiently using local refinement or relocation. In local refinement, more degrees of freedom are added to elements where the approximation is<fiflower quality by either increasing the order of polynomial approximation inside elements (p-method) [7] or by subdividing elements (h-method) [1-5]. In node relocation, the quality of the approximation is improved by optimizing the disposition of the nodes while keeping the number of degrees of freedom constant [14-16]. This thesis presents accuracy estimates for finite element computations of elliptic eigenvalue problems. From results of error estimates for elliptic boundary value problems, error estimates for elliptic eigenvalue problems are established for each mode based on interpolation error theory. An error indicator, which is defined as an.element-wise approximation to the true error, is computationally implemented using two different techniques. a smoothing technique and a direct substitution technique using a differential operator. The accuracy of error indicators is checked using the eigenvalue problem of ainfiibrm, Bernoulli-Euler beam. Two measures of accuracy are discussed. The first measures how accurately the estimator can capture the maximum true element-wise error over all elements. The second measures how accurately the estimator estimates the distribution of true error over the domain. Using the simple node- moving algorithm proposed by Diaz et.al'. [14-16], improved grids are constructed keeping the total number of degrees of freedom constant. The thesis is divided into five chapters. Chapter two is a review of existing adaptive finite element methods. Error measurement techniques and adaptation strategies are reviewed with an emphasis on elliptic problems. Chapter 3 presents an overview of the solution to eigenvalue problems in infinite dimensional space and their approximation in finite element subspaces. In chapter 4, error estimates for elliptic boundary value problems are reviewed and error estimates for elliptic eigenvalue problems are established. As a computable form, two types of error indicators are proposed and their accuracy is tested. Node relocation is also performed to achieve an improved grid. Chapter 5 discusses the results of the research, presents conclusions and proposes further research. CHAPTER II. REVIEW OF ADAPTIVE FINITE ELfliENT METHODS 2.1 Adaptive Methodology The objective of adaptive finite element methods is to adaptively change the finite element model to improve the quality of the finite element approximation. A possible approach to improve the finite element solution is to increase the number of degrees of freedom by subdividing or increasing the degree of polynomial approximation. For example, new degrees of freedom can be added selectively to elements where the finite element approximation is poorer. In these elements, new degrees of freedom are added by subdividing the element or by increasing the degree of polynomial approximation within the element. This process can be repeated until a prescribed accuracy is achieved. We will refer to this process of selective addition of new degrees of freedom as optima]. refinement. Another approach to improve the quality of the finite element solution is to relocate the mesh by reducing or increasing the element size (length in 1-D or area in 2-D) to achieve the best possible grid for the given number of degrees of freedom. This method can be referred to as ogml relocation problem. Both the refinement and relocation approaches need the element information from previous solutions to decide where the approximation is poorer so as to add new degrees of freedom selectively. To implement an adaptive strategy, one must address the following issues 1. How to measure the quality of the approximation. The quality of the approximation is measured by the difference between the exact solution and time finite element approximation. There are two basic techniques to estimate errors in finite element analysis: the residual approach and the WW- 2. How to adapt the solution procedure to improve the quality of the approximation once an element—wise estimate of the error is known. Three methods are available involving refinement (h-method, p-method) and relocation (r-method). Combinations (such as h-p method) are also available. We outline in the following sectfinrexisting approaches to adaptive finite element methods and present a review of the adaptive finite element literature. 2.2 Error Measurement 2.2.1 Residual Error Estimates One popular way to estimate the finite element error is the so- called residual error estimate technique. Considertflm linear. elliptic differential equation. :3 u + q - o in n (2.1) with boundary conditions on 60 Replacing u by its finite element solution uh, one has h . ii u + q - r in 0 (2.2) where r is the residual. It is assumed that the approximationii . . . . . . h . satisfies the boundary conditions. If the approx1mation u is exact, the residual r is zero. Otherwise, denoting the finite element error as e-u-uh, subtracting (2.2) from (2.1) leads to i: e -j2(u-uh) - -r To measure an element-wise error, a local auxiliarv problem must be solved, stated as follows 33 w + q - o in element K (2.3) w - uh on element boundary 8K In element K, the solution w to (2.3) is treated as the true solution. The element-wise residual, rK, and the element-wise error in the * h . energy norm, ||e where e -w-u , can be defined as * IIE,K ' r -m(w-uh)] and K K * 2 h h h ||e I'E,K - fK(w-u ) rK dK - Ix (w-u )jE(w-u ) dK (2.4) Integration by parts using the boundary condition of the local auxiliary problem leads to the final form of the element-wise error. For example, in a 1 dimensional 2nd order differential equation, d2u d2 _dx__z_.+q-O Ofi-‘de 1110 (25) the local auxiliary problem associated with (2.5) is 2 '2';’ + q - 0 in element K (2.6) dx w - uh on element boundary 6K The element-wise error in the energy norm is * * d2e ||e*||2 - f e dK - f (w-uh) d2(w'uh) dK E,K K de K ""7"— dx Note that we have the boundary condition, w-uh = ex= O on 8K Integration by parts leads to the final form * 2 d h 2 He IIE,K=- IKI-a; 1 cm (2.7) In general, the solution w to the local auxiliary problem (2.3) * cannot be.computed exactly and, instead, only an approximation w to w can be computed. There are several ways to compute the approximation * w . For example, one may increase the polynomialOdegree of approximation inside K or, instead, one may refine the element K by ' * subdividing it into smaller elements. Once the approximation w to v: is available, the error in element K can be computed from eq (2.4). The residual approach is complicated by difficulties iJi solving the local auxiliary problem. In addition, the approach is difficult to extend to nonlinear problems and to choose the appropriate norm. For a discussion of these difficulties, see Oden et.al. [27]. 2.2 2 Error Estimates based on Interpolation Error A different technique for error estimation, interpolation error theory, can be used to construct error estimates for adaptjAna finite elements. This approach was proposed by Diaz et.al. [14,15,16]. The procedure is outlined below. Assume that the true solution u to an ordinary differential equation is a smooth function, and let u be an interpolator of u in a I finite dimensional space. By construction. u is equal to u at the I finite element nodes, i.e., u b. - u b. 2.8 I< J> < J> < > where bj are finite element nodes (j-l,2,...,NO) and NO is a number of nodes in the finite element model. When u is the SOIUCUNICD an elliptic boundary value problem, the error of approximating u by its finite element solution uh is bounded by the error associated with the interpolation u by uI. This follows from Cea's Lemma [45]. where c is a positive constant (2.9) h iIu-u II S clllu-ulll 1 The following result is available from interpolation error theory k+l-m llu-ulllm,K 5 c2 hK |u|k+l K (2.10) where m, k and hK denote the order of the variational form, the degree of the finite element solution and the diameter (length) of element K. respectively, and c2 is a positive constant. From (2.9) and (2.10), h k+1-m Ilu'u lIm,K S C hK lulk+l,K (2.11) I Inequality (2.11) provides the basis for element error estimates. For example, consider the solution u to a 2nd order differential equation approximated by the piecewise linear function uh and” m-l). The inequality (2.11) becomes, ||u-uh||1,K s c hK |u|2,K This technique is simpler conceptually as well as computationally than the residual technique. However, when we use this technique, we face one problem: since the true solution u is, in general, not available, we have to use available information, the finite element solution, to estimate it. To estimate the error, we need to calculate higher order derivatives of u since the term |u|k+l K includes (k+l)th O O O o h order derivatives. However, the finite element solution u on element K is a k-th order polynomial. The procedure to extract higher order derivatives needs very accurate post-processing for error estimations. A simple technique was proposed by Diaz et.al. [14-16], and is also 10 used in the present work. A different technique based on rigorous estimates was proposed by Babuska and Miller {6}, and extended by Demkowicz et.al.[5]. Applications of error estimation techniques based on both residual and interpolation approaches to linear problems in elasticity, fluid, and heat transfer as well as nonlinear problems can be found in the paper by Oden et a1. [27]. 2.3 Adaptation Strategies Suppose that we have already estimated the element-wise error of the approximation. We are now in position to change the finite element model in order to improve the accuracy of the approximations The methodology for adaptation can be roughly classified1hnx>two classes : local refinement and relocation. re 'nement Local refinement increases the total degrees of freedom in the finite element model by 1. increasing the number of elements while keeping the polynomial degree of local basis functions fixed (h-method),or 2. increasing the polynomial degree of approximation while keeping the number of elements fixed (p-method) node relocation 11 In the node relocation method (r-method), the total number of degrees of freedom remains constant. Nodes are relocated within a fixed number of elements and with a fixed polynomial degree of approximation. In time-dependent problems, some moving mesh methods (e.g. moving finite element method) have been proposed. These methods have similar basic features to the r-method. combination Some combinations of the above methods are possitflxa (e.g. h-p, h-r, and p-r method). The most popular method is the combination of h and p method. (h-p method) In the following sections, the different adaptation strategies are discussed based on the type of differential equation : time- independent (elliptic) and time-dependent (parabolic, or hyperbolic) problems. We also include discussion of eigenvalue problems, which are the main objective of this work. 2.3.1 Elliptic Problems Since the late 70's, adaptive finite element methods have been applied to linear, elliptic problems, using heuristic as well as rigorous mathematical justifications. Originally developed for linear elliptic problems, adaptive finite element methods have now been extended to some classes of nonlinear problems. We review adaptive methods in elliptic problems based on the h-, p-, h-p, and r-method. 12 2.3.1.1 h-version In the h-method, elements are subdivided into smaller ones while the order of polynomial approximation remains unchanged. Optimallv refined meshes are achieved by selectively subdividing elements where the error is large until a specified accuracy is achieved. This method is based on the fact that, as the mesh size h goes to O (as the number of elements increase), the error in the energy converges to 0 [44]. 'fiuzcomputational detail of this approach can be found in series of papers by Babuska and Rheinbolt [1-4], where the residual technique was used for error estimation. There, error indicators constructed from the finite element solution were used to identify elements where the approximation is less accurate. An impressive work using the h-version of mesh refinement was done by Demkowicz et.a1.[5], where interpolation error estimates are used.tna estimate the error and the extraction formula was modified to calculate highly accurate higher order derivatives. 2.3.1.2 p-version In the p-method, the number of finite elements remains unchanged, ‘while the order of polynomial approximation is selectively increased. The p-version of the finite element method is based on the notion that higher order polynomials can approximate a smooth function better than lower order ones. Distributing different higher order polynomials in regions where the approximation is poorer can produce better overall approximations. The higher order polynomial is used iritflua elements 13 where the error is large. Note that the original grid (element size) is not modified. The p-version was studied by Babuska et al. [7]. The authors obtained error estimates in terms of the polynomial of degree p, and showed that if the convergence rate of the h-method and p- method is expressed in terms of the number of degrees of freedom, the p-method cannot have a slower rate of convergence than the h-methodu In particular, when corner singularities are present, the convergence rate of the p-method is exactly twice that of the h-version. Numerous works have been published on p-methods, especially applications to fracture mechanics problems (see references in Szabo [10]) The p-method has been shown to produce better results than the h-method in elasticity problems with singularities. [7,11,12] .A popular way to formulate the p-version is the so called hierarchical finite element approach. This is a computationally efficient procedure in which the stiffness matrix corresponding to the hierarchically enriched mesh includes the stiffness matrix corresponding to the previous mesh as a submatrix. The hierarchical finite element method is discussed by Zienkiewicz et.al. [8,9]. The paper by Zienkiewicz and Craig [11] includes recent advances in the p- method and hieriarchical finite elements. 2.3.1.3 Combined h-p method Babuska and Dorr [13] studied the combination of h and p-methods, *where error estimates in terms of both the mesh size h and the polynomial degree p were explicitly obtained. An important result 14 there is that the optimal h-type refinement together with properly distributing p's can produce exponential convergence. In summary, it is generally believed that a faster increase in local accuracy can be achieved using the p-method, specially in problems with singularities. The best approach in terms of accuracy is the combined h and p method, which leads to exponential convergence. These results haveibeen restricted to one-dimensional problems and to linear elliptic problems in two dimensions. 2.3.1.4 r-version The h- and p-method improve the finite element approximation by increasing the number of degrees of freedom. In some problems, however, one may want to keep the number of degrees of freedom constant, i.e. improve an existing grid. This is the case, for instance, when finite elements are used as part of an iterative design optimization process. In.this context, the r-method was proposedlnr Diaz et.al. [14-16]. The paper [15] includes some theoretical aspects of the r-method. The work is summarized as an optimal relocation problem where an objective function B is derived based on interpolation error estimates and design variables are element lengths h in 1-D problems or areas A in 2-D problems. The optimal relocation problem is stated as follows (1) l-dimensional case Find the vector of element lengths h-(hl,h2,...,hN) that 15 . . . 2 2k 2 minimizes B (u,h) - N21 hK |u|k+l K (2.12.a) subject to KglhK-i , th o (2) 2-dimensional case Find the vector of element areas A-(A1,A2,...,AN) that . . 2 k 2 minimizes B (u,A) - Kg, AK |u|k+1 K (2 12.b) subject to KglAK-l , AK 2 0 where N is the total number of elements. The authors also derived optimality conditions for (2.12), fK- 1112 a“ < X) XX’ + 2 u x’ + ZX’ =- A u(x.y) (my) 6 0 (5) 6x 6x 8y 6y + boundary conditions on 60 22 In this case, we can get the general terms of the series solution only when the problem is defined over a simple domain. e.g. circular or rectangular with simple boundary conditions. e.g. clamped or simply supported along the boundary. The general eigenvalue problem (including beam. plate etc) can be written in a compact form as follows Find the pair (A,u), where uEECZT such that L u - A u in O M u - O on 60 (6) where L and M are linear differential operators and 0 is a smooth, bounded region with its boundary 80. The differential operator L is a self-adjoint, elliptic operator of order 2m and M is a compatible boundary operator of order m. We will call equathn1(6) the differential form. Equation (6) has an infinite number of solution eigenpairs (Lu) and since the equation is homogeneous in u, amplitudes of eigenfunctions are arbitary and only their shape can be determined uniquely. Admissible boundary conditions are either essential, resulting from the geometric compatibilities, or natural, resulting from moment or shear force equilibrium. For example, in the case of the Bernoulli- Euler beam of length 2 simply supported at both ends, boundary conditions are essential : u(O) - 0 u(t) - O 23 d2u(2) = 2 d “(0) - o E 1(2) 2 dx dx 0 . (7) natural : E I(O) In the case of the thin elastic plate clamped along its boundary an, essential : u(s) - 0 _g%_(_§l - 0 s E 80 (8) natural : none 3.1.2 Weak (Variational) Form The differential form (6) requires that u and all of its derivatives of order less than or equal to 2m be continuous at every . . . . m ’ . . . . . . pOint in 0 (i.e, u is in C2 ). This is a difficult condition to satisfy. A weak or variational formulation relaxes this requirement and facilitates computations. This weak form will be used to obtain a finite element approximation to the eigenfunction. Consider again the differential form (4) and boundary conditions (7), representing the behavior of the Bernoulli-Euler beam simply supported at both ends. We multiply a weight (test) function v defined on (0,2) to both sides of equation (4) and integrate tuner the domain such that the differential equation with boundary conditions is satisfied in the sense of weighted average, i.e., fév[%2-2—(E1(x)iz%) - Au 1 dx= o in (0,2) x dx + boundary conditions on 60 . (9) 24 The equation (9) includes 4th order derivatives of u, whereas no derivatives of v appear. Integration by parts of (9) leads to 2 d2u d2v 2 f0 E I(x)—7—7dx - A [O uv dx = o in (0,2), (10) dx dx provided that the solution u and the test function v belong to the class of admissible functions, denoted by V, defined as v - ( VEHZ I v(0)-O, v(2)-0 ) In words, \I is the set of HZ-functions (bounded energy) that satisfy the essential boundary conditions v(0)-0 and v(2)-0. A function is said to be in H2 if the function and all of its partial derivatives up to 2nd order are defined in a square-integral sense (see Appendix A for more detail), i.e., H2 - { v | |v|0,|v|l,|v|2 S M < w I 2 where |v|0- févzdx, |v| 1- f20(%2dx , |v| 2= [20(3-1 2dx and M is a dx constant. Note that essential boundary conditions are included in the definition of V and that the smoothness requirement on u is weakened. Now we can write the weak form of this problem as follows Find uEiV such that 2 d u d2v 2 fOEI(x)-——2- zdx-Afouvdx=0 forallvinV (11) dx dx where V is defined as before. 25 In the case of the thin elastic plate clamped along its boundary, the weak form can be derived follwing similar steps. After introducing V and integrating by parts, the corresponding weak form is find uEV such that 2 f0 V2u V v dxdy - A In uv dxdy - 0 for all v in V (12) 2 2 a2 62 where V is Laplacian operator defined by V =-——§ + -—7-. 8x 6y The weak form can be written in a compact form by introducing inner products. The general form of variational eigenvalue problems is find uEV such that a(u,v) - A(u,v) - 0 in 0 for all v in V (13) where a(.,,.) is a symmetric, bilinear energy inner product and ( ,.) is a symmetric, linear inner product. In the case of the Bernoulli-Euler beam of length 2sfihmly supported at both ends, 2 dzu dzv 2 a(u,v)-f E I(x)———-———-dx , (u.vr-f uv dx in 0(0,2) 0 dx2 dx2 0 u and v should satisfy u(O)-v(O)-O u(2)=v(2)=0 In the case of the thin elastic plate. a(u,v)-f0 Vzu V2v dxdy , (u,v)- [nuv dxdy in 0 2 2 2 6 a where V - + .;;§ 3;?- 26 8u(s)_dv(s):0 u and v should satisfy u(s)-V(s)-0, an an sGdO The differential form (6) and weak form (12) lead to an infinite number of eigenpairs (A,u). These pairs can be«ordered according to the magnitude of eigenvalues under the.assumption of no repeated eigenvalues, i.e., A1< A2< A3< 121 qi ¢1 (15) ‘where ¢§K»(i-l,2,3,4) are local shape functions (cubic polynomial) in element K. Patching together shape functions ¢1. ¢3 and d2, oé‘over the domain lead to the global basis function u(x) and m(x) (K) {K) and q3 represent the respectively. The nodal values q (K) 4 represent the and q displacements at nodes i and i+1, and qéK) slopes at nodes i and 1+1. Substituting the element-wise approximation (15) and the test function into the weak form (14) leads to an element matrix, 30 [K] q - Ah [M] q (16) where [K] and [M] are the 4x4 element stiffness and mass matrix. Summing up the contribution of each element, one can get the global formulation, [KIGQ - Ah [MJGQ (17) where [K]Gand [M]G are the N by N global stiffness and mass matrix. Equation (17) is a matrix eigenvalue problem. The approximated eigenfunction uh can be computed from the approximated eigenvector Q and global basis functions d, i.e., h ” '1g1 Q1¢1 The matrix eigenvalue problem (17) leads to N approximated eigenpairs (Ah,uh). As before, ordering them according to the magnitude of approximated eigenvalues, i.e., h h h h h A1 < A2< A3< ... < AN_1< AN the approximated eigenvalue A? and corresponding eigenvector Q3 are the 2-th mode eigenvalue and eigenvector The 2-th mode approximated eigenfunction u? can be computed from h “2 '1-1 Q2121 Equation (17) can be also written in terms of the 2-th mode as follows. 31 [K1G Q, - A? [MIG Q, (l7.a) 4.3 Summary and Discussion In the previous sections, we have gone from the differential form (6,6.a) to the matrix eigenvalue problem (17,17 a). Closed form solutions to the differential form (6) or to the weak form (13) are not generally available, which leads to approximations in finite dimensional subspaces Vh of the original space V. We selected Vh as a finite element space. In general, Vh differs from V and the finite element approximation is different from the exact solution. Our concern is to reduce the error between the finite element . . h . approx1mation u and the exact solution u. To reduce the approximation error, an error indicator should be constructed first so that we can estimate and reduce the error based on the error indicator. This question is discussed in the next chapter. CHAPTER IV. ERROR INDICATORS FOR ELLIPTIC EIGENVALUE PROBLEMS 4.1 Review of the Error Estimates for Elliptic Boundary Value Problems Let u be the weak solution to the variational form of the elliptic boundary value problem of order m, a(u,v)-(f,v) for a11.\r in V, where V denotes an admissible space. Then u minimizes the functional (potential energy) I(v)-a(v,v)-2(f,v) over V. Let uh be a finite element approximation to this problem. Then uh is the solution to a(uh,vh)u- (f,vh) for all vh in Vh, where Vh denotes a finite element subspace. It is also the minimizing function of the functional I(v) over Vh Theorem [44] (a) a(u-uh,u-uh) - hMin h a(u-vh,u-vh) for all vh in Vh v in V . . h . . TUnis means that, measured in the energy norm, u 15 the best pOSSlble member of all the members in the subspace Vh. h h (b) a(u-u ,vh) - O for all vh in V 32 33 The error u-uh is orthogonal to all the members in Vh, i.e., uh is the . . h . prOJection.of u onto the subspace V with respect to the energy inner product. See Figure 1.1. The distance between u and uh in the energy norm is bounded by c2h2(k+1‘m) 2 for uEHk+l(0) (l8) a(u-uh,u-uh) 5 |u]k+l where k is a polynomial degree of approximations in subspace Vh, c is a constant independent of h, m is an order of the variational form, denotes the semi-norm in Hk+1(0) and h is the size of the l'Ik+1 largest e1ement.[44] In l-D, h is the length of the largest element. In 2-D, h is the diameter of the largest circle inscribed in the largest element. v; -----n Figure 1.1 Approximation error in elliptic boundary value problems (measured in energy norm) 34 The exponent OflllIl(18) indicates the rate of convergence as meshes are refined. The error can be reduced if h becomes smaller via . h h mesh refinements : as h 4 0, a(u-u ,u-u ) -* 0. As more elements are used, the finite element solution uh converges to the solution u in the energy norm. This is the key to the h-method. The term [u]k+l reflects the smoothness of soluthnL Suppose that linear interpolation functions (k-l) are used to approximate the solution to the 2nd order elliptic problem (m=l). Frmn(18), the approximation error is bounded by c2h2|u|§ and the convergence rate in the energy norm is 2. If u is linear, the error is 0 since the |u|2 term vanishes. This is the case in truss problems with concentrated loads. ‘If u is quadratic, the approximation error is not 0 since the [u]2 term does not vanish. However, if we used quadratic interpolation functions (kr2), the approximation error would be bounded by czhalu]§ and, thus, the error would again be 0 since the |u|3 term would vanish. Higher order basis functions can approximate u better. This is the key to the p-method. Equation.(18) is the key for the construction of error estimates for elliptic boundary value problems. The same equation.will be applied to elliptic eigenvalue problems. 35 4.2 Error Estimates for Elliptic Eigenvalue Problems In elliptic boundary value problems, the finite element . . h . . . h . approx1matuniti is the prOJection of u onto V and is the closest member to u in Vh. However, in.elliptic eigenvalue problems, due in Vh is the closest approximation to the 2-th mode eigenfunction uz Rayleigh projection Pug. The Rayleigh projection Puz is defined as follows If u is the 2-th mode eigenfunction (solution) to the 2 variational eigenvalue problem in V, then the Rayleigh projection Pu£ is its orthogonal projection in the subspace Vh , i.e., a(ui-Pu£,vh) - 0 for all vh in Vh. By definititni, the Rayleigh projection Pug is the closest approximation in Vh to ufl measured in the energy norm. See Figure 1.2. Since the projection of ug onto Vh is Pug, the error bound (18) for elliptic boundary value problems suggests the following error bound for the 2-th mode eigenfunction uZ for elliptic eigenvalue problems c,h2(k+1-m) 2 (19) a(uz-Pu£,u2-Pu2) s Iu2lk+l 36 “I. II \/ I , I _.~' ] a(uj'Puth‘Pul) I I I I Pu] h “t Figure 1.2 Approximation error in elliptic eigenvalue problems (measured in the energy norm) USing the equivalence of the energy and Hm norm, equation (19) is equivalent to k+l-m [lug-Pu!”m 5 ch |u£|k+1 (20) where ||.||m denotes the Hm norm (Sobolev norm). The distance between u! and u? in the energy norm is expressed as h h h h a(ul-u£,u£-u£) - a(uz-Pu2,u2-Pu2) + a(Pug-u2,Pu2-u£) . (21) If Vh is a finite element space, the Rayleigh projection Pug is not computable in general, but is 'close' to u?, i.e.. u? z Pug. Thus, if 37 the term a(Puz-u$,Pu2-qr) in (21) can be neglected, from (19) and (21), a(uI-uh ug-uz) S c'h2(k+l m) £. (22) In 12 2 k+l This error bound (22) allows us to set up the explicit form of error estimators for eigenvalue problems. 4.3 Error Estimators To improve the quality of approximations efficiently, new degrees of freedom should be added in a selective manner to elements where the approximation is poorer. This requires error information at the element level. Define an element-wise estimator of the true error as . 2 . an error estimator and let 6K denote the error estimator of the 2-th mode in the element K. Then, the error bound (22) suggests the following form of the error estimator. 2 2a 2 . 6K - c hK lu2lk+1,K in l-D (23.a) 2 a 2 . 6K - c AKlu2|k+l,K in 2-D (23 b) where hK is a length of element K, AK is an area of element K, 2 is a mode number, k is the order of local basis functions, 38 a is the order of convergence (=k+l-m). c is a constant independent of hK' A and 2, K k+l . denotes the H semi-norm on element K. and |'|k+l,K Note that the constant c depends only on the element type. In general, this constant is not available. For example, suppose that Hermite cubic interpolation functions (standard beam element, 2 nodes and 2 dof per node) are used to approximate solutions (eigenfunctions) to the vibration of a uniform, Bernoulli-Euler beam. Then the 2-th mode error estimator on element K is 4 2 2a 2 4 d u 2 ex - c hK |u|k+1’K- c hKfK( t/dxa) dx a II rx.) The error estimator of the form (23) is not still usable for computation since it includes the exact solution ug. To make the form (23) computable, it is necessary to replace ufi by a known function A . . . . h that can be computed from the finite element approximation ui. Let u£ be such function. Then, the computable form of the error estimator, denoted by 0:, is 2 2a A 2 "K - hK lu2lk+1,K a-k+m-l . (24) We call this an error indicator. The error indicator should be A available from finite element approximations. Different forms of u will be discussed in the following sections. This is the main result of this work. 39 4.4 Computational Implementation of Error Indicators Recall that the error indicator has the form of 2 2a A 2 0K - hK Iu2lk+l,K a=k+m~l . (24) A The function u should be, first, computed easily from the finite 2 A and, second, such that the ratio lu2'u2lk+l,K / element solution u: |u2|k+l K is small and asymptotically correct. i.e.. converges to zero as h a 0. In this section, two techniques to construct A U2 are discussed. 4.4.1 Smoothing Technique This technique was originally proposed and successfully tested by Diaz et.al. in [14,15,16]. Using this technique, the function.uy is constructed directly from the finite element approximation. For simplicity, let us illustrate only the l-dimensional case here. Irrea l-D problem, an approximation u? is a k-th order polynomial over each element where k is the polynomial degree of local basis functions. The k-th derivative of uh is a constant that may vary from element to 2 element, and the (k+l)-th derivative vanishes inside the element and A is undefined across elements. The new function ufi is defined such 4O dkui dkuh that the difference between k and vanishes in the weighted d x d x average sense, i.e., dku dkuh 2 2 f0v( k— k)c1x==o d x d x where v is a weight function and k is the order of the basis functions of us. Let us illustrate the procedure using an example. Consider cubic beam elements over the domain 0-(O,l). The approximation u: is a 3rd order polynomial on the element. Over the domain, the 3rd order . . h . . . . derivative of u! is a pieceWise constant function. The function u satisfies 2 l d3u£ d3u3 f0v(——3-— 3 )dx-O. (25) d x d x To simplify the notation, let 22 and p2 be the 3rd order derivative of 62 and “3' Then eq (25) becomes fév< «a; 2‘“ 1 s 1 l H a- 1... .. H b- ] b—i ‘- H e- 1 1 r V—. r r r ' I I r —‘ J 5 7 9 11 13 number of elements CL quintic 1000.01 umnf 1 M 100.01 \‘ >4 3 1 \ '\ Q 5001 \\ :i . ‘1 E? 1 \ \ \\N\ #4 1 a) 3 \ u .5 10.0? . M .5 V 5.0-;1 V X 1 x «ed a” 1 \ e- 1.01 X 0.5.; H 5", \ 1 I-I b-l \\. H a-s ‘ H 41-4 \., \ 01 '3‘“ fifi number of elements 555'5f1‘1' CL quintic 1000.0 1 500.03 \ \ A\\\ q\ g\ k \ 100.01 x ' 3 W \ \ \ 50.0: \ \ . \ . 1 x \D 10.01 - \ \ \ 5.0.3 \ \ .. 1 \ \ \ ‘ \ \ 1.013 \ . ‘1 0.53 H 0"? \ \ 1 H #0 \ 1 0—0 0-6 X 1 H 0" ‘1 ‘ 01 ’3‘? 1 ,. 9 , r I! 5 7 9 11 13 number of elements CL quintic 1000.0 500.0 1000 \\ \\ R\\ 500 X \ ‘\n \ \ \‘ . x 100 \ \\ \ , 5.0-1 \,‘ \ \ 10 \\. 3 read 1 (L5 I—O b-O “\ H 4-5 '\ hi‘“ ‘ 0.1 ‘1 db: f j' T'— I Y r f T 3 7 9 11 13 number of elements Figure 5.2 Error in eigenfunction vs. Number of quintic beam elements (clamped-free boundary condition) 59 CL-SS quintic M H C-C H 6-! H H ‘ e-vo 4-3 _ 41 0 11 1’0 13 '74 number of elements 28‘ "IRX (K 8 A A -1-“ a. O A 4 AlAAAA. *- CL’—SS quintic CL-SS quintic “ \\\\ E JALJ §'§‘ mflx (interp 1)K A “5.5.3 A n max (interp 2) -35.} 101 \ 10? 1'»i . sf 1 H h. + H H * rasbfl H H 1 ‘ H H H H ‘ ’:. “3* fi' , I v 1 1 a ‘1’, ' t r f l fl 41 5 5 1'0 12 14 4 a a 10 12 14 number of elements number of elements Figure 5.3 Error in eigenfunction vs. Number of quintic bean elements (clamped-pinned boundary condition) max (interp 1)K 60 CL—CL mnnfio : \\\\\_ a: 50 W X Q4 E 4 101 I 5.1 1 ‘HH 4H” HH ,HJ-a ' I rj—V I V fl :1 7 o 11 13 number of elements CL_CL quintic CL" CL quintic é u 0 .5 o 3 me (interp 2)K (I a H 5-0 H H H H ' H 6-8 ' ' v f ' I ' I v I w ' 5 1'1 13 11 i 0 11 13 number of elements number of elements .1111 0.3111 1\\\\ ‘0 \\X Figure 5.4 Error in eigenfunction vs. Number of quintic beam elements (clamped-clamped boundary condition) 61 General conclusions concerning estimates of maximum errors using indicators can be made as follows (1) The maximum error of the l-th mode eigenfunction is the same for all modes 2 when the number of elements N = (22 + b), (32 + b), etc. is used. The constant b depends on the boundary conditions. It equals 0 in 55-38 case, -l in CL case, and l in CL-SS and CL-CL cases. This is true for both cubic and quintic elements and is true for the maximum error estimator and maximum error indicators. (2) Both error indicators can estimate accurately the maximum error in the i-th mode eigenfunction when more than N= (32 + b) elements are used. (3) Since, in general, fewer elements are preferred for practical reasons, the numbers of elements between (23 f b) and (32 + b) elements may be used when a certain amount of error can be accepted. 4.5.2 Accuracy on Error Distributions In this section, the accuracy of error indicators is tested based on how accurately they can estimate the distribution of true errors over the domain. As will be discussed in the next section, the node relocation to achieve an improved grid within a fixed total degrees of freedom requires accurate estimates of the distribution of true errors over the domain. 4.5.2.1 Cubic Beam Elements 62 Figure 6 includes results concerning the percentage error distributjxni1ising cubic beam elements. Let EK(%) and nK(%) denote the percentage error in the energy norm and percentage error indicator corresponding to the K-th element. The percentage error distribution is computed by EK(%) - BK / Emax X 100 , nK(%) = UK / ”max x 100 In Figure 6, three kinds of bars are shown . The height of the bar at a given location corresponds to the percentage measure ( EK(%) or nK(%)) at that point. In symmetric cases (SS-SS and CL-CL cases), only half of the domain is shown. The legend of three bars is descibed below. EK(%) = BK / Emax x 100 K 2;??? (interp 1)K/ max (interp 1)K x100 l (interp 2) / max(interp 2) x100 K K K _ .—.__. ._-—_——— Legend for Fig. 6 and 7 63 The interpolation type 2 indicator can eStimate the distribution of the true error more accurately than the interpolation type 1. In the 88-83 case, the interpolation type 2 indicator can estimate the distribution of the true error accurately regardless of N. However, in the CL case, N - (3.2 + b) elements are needed to accurately estimate the distribution of the true error. This can be observed also in the cases, CL-SS and CL-CL. We can conclude that the interpolation type 2 indicator can estimate accurately the distribution of the true error when, in symmetric cases, more than N = (22 + b) elements is used and, in unsymmetric cases. more than N - (3£ + b) is used. 6M 3rd mod. SS—SS 2nd mode SS—SS ll/l/l/l/l/I/lJ/I/I/I/I/ 7.0.0... co. «0... as. unoca. lupus-o. vooooo’oooooooooooo coon-...... on. capo-0.0.0.0.... .oonou-u ...-o. H‘filfl 1.43.1 41”” 2 ,yvvvuV/w/vuw/VJJIVJIV/é ...................r.............................u.n.u.u.u.. w w w w m o Mammmmmwzzzazorz Aaazz ..........-a..yo.on...uononouonouououooouououo.ouoaouonououf mwmwwmzmwMme on uouum Location along the beam Location along the beam ne1x=10 nelx=10 SS-SS III/IIIIIIIIII/I’II’IIll/’III/a w w w 1. 2 any nouum w o ..o. «h mode MmmmwmwzzZZZZZoZazav .. o o o o . o o o o o u.ao.o.a.........ououoo...u...oo..¢ool o. ,v77%.?vayvvvvvvuzvvvvvH ..292wmmmm A: nonum e b e h t 80 “1| 0: 1X 51 e n on .1 t a w. e b e h t 0 gl n: 10.x am nn 0 .1 c a w. 61!: mod. 83-88 £7.[u/vvvvvvu/uvvvvvvvw/vm 7thxnodc 38-88 .u/VJ/V/V/l//////////.////.I. on a. .nffo.co..ono..uononanon-nououououonononouovouc- A.v “out” .zmeMJMJmmmmmm ........... I.J/V/////////////////////////A o. o o coo...500......coo..........-oucuouoooouoo-oooooooco. V/III/lll/I/l/ll/II. hon."on.”.uouhouououououonouou.nfiononououh on wouum Location along the beam Location along the bean nelx=10 nelx:10 m a e b e h t 00 n O 1 a) 5 nt 0“ .18 cm 88 C1 we 59 Vb «LC 0.1 Yb ru CC I.\ e g a t n e C r e P 6 e m. i F number of elements nelx 65 «a and. fiZQZQZéZfiZéZéZZZfiZZZ 3...}... .I...f.... .ZQZQZQZQZZZZZ . u u u a . ...oual. ...................ulo...o... ...... A.c tout” Location along the beam nelx=15 Location along the beam ne1x=13 owl/V5417 I/l/t/f/Jl/x/I I47 747/ llllll x} ...... m 2 .. . _ e n; .h 2J llllllllllllll t In Coho n 1 avg/ova. m = O a m w n m o m .1 t m a AOV nouum Adv uOuum Location along the beam ne1x=10 Figure 6. (cont'd) 66 5th :11on CL . a a a o a u o . covOOOIODOOCOOtiu 7///////////IZ “hon.” aonouououonououonon.uonoconouonononouououououooo-ouou OOCOOQIOQQ r3333??!:? Vl//////’///’///////A ononononoufiouououo”anon.”ounuouounnouououou I‘OOOIOOOOUOCO.‘ ...»...u....o.o..v....o.o.’o.o. o.oa.~.c.cooaon. .00....novooooououosovovo-o-o-o. I v/é/J/vvvvvvvvvua .....cccoo. ..OODOIOIOOOOQOODOUO~ .oanccacaqaocca-oaoac-a. O.I....0000...IOOCDOIOOOOOOOOOD.IOOOOODOOOCODO- OHSMGOE hOHMO Location along the beam nelx:10 6th mode CL anon-o- ...-ocooo-oqou ...... ....oooo ...-acqacoooau .oopsuo- ....... .... ...-noo-coouoooo Ohflnflva MORHO . ......cuaancn. «35.45................. le/I/ll/IIII/ll. a .o-o...o.-oocaoo-o-unoooooooooo coccaooqoo: ......ou....-..o...o ~ .aa ccccccc nonsense-cocoon... 0.0.9.0.... coo-on.o.a..0uooonocope-05.05.990.338 ......cooc.a:.aoo. O.COO....FODODOOOOOOOOODOI.DOIO'OUOIQ...’DOO.DOOO’ODO~ QOOIIIOIICIO ......0IODOUOOODOOO.ODO- - o c ..... o a o a o o a o a o o o o o u o o o 0.0 o o o o o.o.to....ooovouooo-o-cnoeouo. 2,2,xzvvvvvvzvvvvvvvvvvvva o o o o . o..........-oc...c.o..o¢oa.ooao. 1O ... ..... O‘OCIIOQlCCOOO'OICOCQO oOto...QO.IObODOooOObOOOIOO.IOuODOOO>OIOIODODODCDODOCODOOOo 7 o. v. so... Location along the beam ne1x=10 CL uth mode 7/4 .................. .... CDC00......-...IDOOIODODOI...IOO§OOO.IO&OIO- .................... vvvvvvvvvvvvvvvvvvvvvvm .IV/l/l/l/l/l/l/l/l/l/I/A cocooaoo.ococoaaocncuo. v.0O.OtCDOIODODOOOOODOVODOOOUOubb.....DQDODODOIO. W/vvvVVVVVVVVVV/VVVVVVVVVM VII/ll/lll/l/Z/l/Vl/V[VIII/Ill. ...... ...................... 0.0.0.0.....0’050-ooc00-O-0nOooso.on...-ooowoottUDODOOOnoton .I III 1 Iii ...... o a o a . 55.33.5555... a a a a a . w w w w 030606 HOP—0 Location along the beam nelx=1u (cont’d) Figure 6. 67 (BK—I) 5th mod. OI ..... n Ill ...WWWNHQAWE 3 . ... m yvvvuvvvvvvvvvééé .5.» mum”We.......u.u.wm.u.u..m.w.m.wu.u 1 III—1 .oooo W «IvVVVVVVVVV/III/I/l/I/II/é on... o a o o ooo o o on.oouonononoucnououououououonon.“ .......... .u..n.u.n.u.x...u.u I a f I 5 .rVVVVVVVVV/t/If/I/I/l/I/IA ....................3".”.N.H.H.n.u.n.u.n.u.u.u.s.” yaw/a 32595fi avg/u,vvvvvvvvvvvvvvvvvm avvvvvvvvvvvvvvm va uouum ‘ Location along the beam nelx=1u 8th mode | Z////////l. ........... ......aoocauoacscz DOUG-...... OIOIOIOI'OOIODIQOIODOooworOOOIOIOAOOOOO-C-OIO- vvvvvvvvvvvvu/yvvvvyvvvvvvm 7/////.//////.////Jlm ............. {COO-.0...DOOODOIOIOOQIOOODOOOIO- aI,xxx/zyvvvyvvvvvvvvu/vm I VIA w m 0 mu» m Aav uouum ‘ 1| . . . . ......r a i Location along the beam nelx:1u (cont’d) Figure 6. 68 3rd nod. CL— 33 3rd mod. CL-SS UMSMGOE hOhhv Location along the beam nelx 7 :6 Location along the beam nelx 3rd nod. Ohflfldoa hOLHO Location along the beam nelx:8 (31+!) 81-.“ Gr- Okflsoa kg. Location along the beam :10 nelx (cont’d) Figure 6. 69 (21+!) Cb-CL 4th mod. w w w m Ohdmflvg hOth OHDnGOE hOth Location along the beam Location along the beam nelx=9 (31+!) =5 nelx uh mod. Cir-CL 4th mod. Cb—fl. Ohfinfloa hOhho Location along the beam Location along the beam nelx=13 =7 nelx (cont'd) Figure 6. 70 4.5.2.2 Quintic Beam Elements Figure 7 includes results of error distributions using quintic beam elements. Using these elements, the interpolation type 2 indicator estimates very accurately the distribution of the true erroru Even in the CL case, N - (2 + 1) elements are enough to estimate accurately the error distribution for the 2-th mode. Conclusions concerning the accuracy of the distribution of error indicators over the beam are as follows 1. The interp 2 type indicator can estimate more accurately the distribution of the true error than the interpolation type 1. When higher order elements are used, the interpolation type 2 estimator can estimate error distributions very accurately. 2. Symmetric boundary conditions have an important effect on the number of elements needed to estimate accurately the error distribution. Using cubic beam elements, the interpolation 2 type indicator can estimate the distribution of the true error when N-(22 t 1) elements are used in symmetric cases and N=(3£ : 1) elements are used in unsymmetric cases. 71 6th node enauuda 88-33 3rd node 100.. go. ‘0 ‘0 20- 0 .66.....o...o.o.o.o.o.o.o.o.o.o..... ..Kflaovvvvm.u.m.u.u.va.u&.~.4 o...o.o.....o.o.o.o.o...o.o.o.o.o.o.......o.o.o.o.o.o....... o4Sow.wowowowowowoaowomomowfl.45odomowomomomowomomoaouowon (41414411410 w w w m o iOO‘ on uouum Momma ‘ 100‘ on nouum Location along the beam Location along the beam Location along the beam :6 nelx =5 nelx =u nelx 7u::nodo onitnoda on nouum id» .... .......‘O'.............C.............C.... ., .94.”.a.”5&3.movmonowomowowomomflhono?» .. q i a 4 a l a 4 1 4| 14'" m w w w 0 0 2 ‘ Avv nouum a i a 4 . ¢ a 4 q 41. m w w m N _o 1 Location along the beam Location along the beam 8 nelx =7 nelx e b C h t to n O) 13 Ct nu he on 13 c1 ‘9 w.” e b C 1 t n 1 U q I\ Figure 7. Percentage error vs. 72 CL 3rd mode CL 311! mode (21 ’I) 100: 100‘ ‘ J ‘0‘ w. E . s . '0‘ V .0 3 ‘ .2 1.0. l H z" 4 a 401 .3 33, 40- 7 4 a, 4 20- .5 20.. a ':I < .I . .I o .I 7 7 f o« . a 1 :2 .l 4 o I :1 . a as Location 81008 the beam Location along the beam nelX=u nelx=5 c1. 4th mode or. 4th mode {21-1) 1M‘ r ‘w‘ I . 5 4 004 2 so« A , .. . 5 , ' a I V 00-4 I; {5 00-1 *- 1 .:: :5 , 2 ‘ 5: 5;; f 35 “H :3 2'; ‘°‘ o: o:’ l 2 25 l 204 f: 5:? 204 :-: _o:I < z. 3; ‘ 0' 0" 7 .. , o‘ 7 "o a a 4 4 : 1 4 n u .v Location along the beam Location along the'beam “813:5 ne1x=7 Figure 7. (cont'd) 73 51h mode 555¢55555545411 m z. ...... :.:.:.:.:.:.:. e b :1. e . h ll C rII(Itltlldltrtllllttlllllfl. n ( o 1 a O u m I .1 t m a C . 2&35‘55‘ m a u a n 4 Any nouum Ill-Illlllllllllla szazszazazezazoxazgz ...-....o........o.-o.o-o..o III/’I’IIIIII’IZ. {...-i.e...e.’ ....o.o.a...o...nuou-.o.p.oqo.ouo.oco.e.ou Aav uouum =9 nelx Location along the beam :6 nelx (22-I) on: node on nouum 11 Location along the beam nelx: (21-/) 7th mode AOV hOHHm Location along the beam nelx=13 (cont’d) Figure 7. 71+ 4.5.3 Node Relocation : An Improved Grid In this section we discuss the construction of an improved grid by node relocation within a given number of degrees of freedonL ii node moving algorithm that can produce improved grids is also presented. Since it was shown that the interpolation type 2 indicator is more accurate for error distributions than the interpolation type 1 in the previous section, the interpolation type 2 indicator is used in node relocation. . I In the node relocation method, the variables consist of the nodal positions as well as the nodal variables (displacements and slopes). The objective is to arrange the nodal positions so that the best possible approximation for the given total number of degrees of freedom is achieved. This has been formulated in an explicit form by Diaz et a1. [14-16]. An optimal relocation problem for eigenvalue problems in a 1 dimensional case is described below. Find the vector of element lengths h-(hl’h2"""hN) such that . . . 2 20 2 minimizes B (u£,h) - K21 hK lu£|k+l,K (39.a) subject toxglhK-l , hK 2 o The alternative form of the objective function B2(u2,h) is 2 2a 2 B (u£,h) - Mflx { hK Iufilk+l,K } (39.b) The 'true' optimality condition associated with (39.a or 39.b) is 75 2 2 2 fK- hKalu£|k+l,K = constant for K=l,2 ..... N (40) A Equation (40) can be used in computations using a known function u instead of u, i e., ‘1 2a A 2 2 fK - hK 'u£|k+l,K - "K - constant for K=l,2,...,N (41) where n: denotes the error indicator of the fi-th mode in the K-th element. Equation (41) will be referred to a 'near’ optimality condition. It indicates that the near optimality condition canbe achieved by relocating nodes such that error indicators are same for all elements. To achieve the 'near’ optimal (improved) grid, however, an appropriate node-moving algorithm must be proposed first. A simple node-moving algorithm was proposed by Diaz et.al. [14,15,16] . The ‘basic idea.is that a 'mass-like' quantity is assigned to each element and is assumed to be at the element's geometric center. The mass assigned to an element is proportional to the element error indicator. The mass 'attracts’ the nodes surrounding it. thereby reducing the size of elements with high error. In this work, a slightly modified version is used. We assign a weighting factor to each mass srufli that nodes cannot move too far from their prevhmuspmsitions. The algorithm used in this*work is presented below for a 1 dimensional case . 76 Let Ki denote the i-th element and Xi denote the i-th node in l-D finite element model. The motion of the node Xi is described by * x + *2 fl new ”Ki—l i-l "Ki 1 * (”Ki) Xi = >2 + * ”Ki: —————L. (42) "K. "K. 1 i-l i new . . . . where Xi is the new location of the node xi. Li 18 the length of ’H1 is the measured error indicator of element K, X1 is i element Ki’ the geometric center of element K and B iseaxfiflgmt to be assigned. If 6 >1, node movements are accelerated, i 63.. nodes nmnna far from their previous positions. If 6 <1, nodermnmmmnts are deccelerated, i.e., nodes cannot move too far from their previous positiorm; The Choice of the value of B may depend on the type of problem and thus, some preliminary tests with the different values of fizmgrbe necessary. In our problem, the weight 6 between 0.3 and 0.5 was proved to be a proper value. This process is repeated.unti1 the measured error indicators satisfy the near optimality condition (41) , i.e., "K — constant for K=1,2,....,N We now turn our attention to the effect of finitmzealement nodal positions to the approximation error of eigenvalues and eigenfunctions. A simple test can illustrate tfliis. The rnxflalem is stated as follows Suppose that a finite element beam model of total length 1 has only 2 cubic elements and 3 nodes. The first and third nodes are fixed but the second node is free to move along the beam. Our concern is the variation of the relative eigenvalue error as the second node moves along the beam and also the error in eigenfunctions in the energy norm. For this purpose, two kinds of plots are shown in Figure 8. One depicts eigenvalue errors and the other eigenfunction errors in the energy norm: (A? - A2) Eigenvalue error (%) - ————:——————- x 100 where A2 , A; are the 2nd mode eigenvalue and its approximation and |El— E2| Eigenfunction error (%) - '7I777I73_ x 100 1 2 where I .| denotes the absolute value and EK (K-l,2) denotes the exact error in the energy norm of the K«th element. At the optimal location, the error in elements 1 and 2 is the same, i e., E1- E2. At this point, the eigenfunction error (%) is 0. Vertical lines in Figure 8 denote the point where the minimum occurs. In this test problem, symmetric cases (SS-SS and CL-CL) are special since a uniform grid of lements is actually optimal for the 2nd mode. In general, however,the optimal node location does not coincide with the uniform mesh. 78 33-53 2114 M“ 55-55 2nd mode ‘0 4O 1 j s-O.3 ['03 A 1 o ‘ A 4 4' 33-4 V 30., v 4 c ‘ ‘ O 4 0 , -—4 ,5 3. 4 U . Ir a 304 .5 g 20-4 H > ‘ a 4 a s , gm 4 g 33 I L‘ 5 I 1 " ‘ 1 2° ‘r'r*"l’"*T*" O fvrrvvv Yivvvvv OJ <1‘ 05 ¢10 OJ' OJ 414 ca (16 DJ “04“ position nodal position CL-CL 2nd mode CL-CL 2nd mode 230 604 j A 4 8.0.3 4 8.... . ‘ 4 v 4 2 ‘90‘ c ‘54 v 4 O 4 c e ‘ U l r' 3 1 5 g 1 ad “ § "°1 .. .3 ”i 3 c 4 o c ‘ L4 0 4 H 0 ‘ 0 no 4 t; f ‘3 4 ... d a 110‘ 0 1 4 4 4 4 70‘ o‘ -.. "for. flew” -..,n. 0.3 034 0.: OT 0.7 0.3 0.4 0.: 0.3 0.7 nodal position nodal position Figure 8. Error in eigenvalue and eigenfunction vs. Nodal Position 79 2nd mode CL 2nd mode CL .i 7 a O 1am u .o ‘. v ... .s 3 IO. ‘ T. fio L .3 O1i‘oq1‘40‘1j‘114d1411 O. o o o a 6 ‘ 2 o on ::_uo::u=owuo Cu HOHHO 7. o r we“ 0 .0 ‘. v 0 :3 . In I '0 v .4 v0 til .3 jIJIJII- 2 9 6 J Onv on 03~o>comue cu nouns nodal position nodal position CL-SS 2nd mode CL-SS 2nd mode 0. o v u .5 4 mm a . ‘ .I4 '9. 11.11... .3. 4 0 w w w x on e:~e>:euae :«.uouue 7. V.9 v 7‘ Wu «. ms .4 .o I r VA .0 1 4414-4444d441444414q4444 M w w w n m c on co«UOCSHCOUuo ca Houuo nodal position nodal position Fig. 8. (cont'd) 80 Other results related to node relocation are illustrated in Figure 9, where numerical results associated with the improved grid with respect Us a single mode are presented. There, two kinds of plots are shown. One depicts the error distrflmnjLWIof the unifonn and improved grid based only on the interpolation tvpe 2 indicgtor. The other depicts the error distribution of the met error in the energy norm using the same grids. The dotted line is the error distribution of the uniform grid and the solhiljmm:the error distribution of the improved grid. A 10 element mesh is used in all cases except in the CL-CL case, where only 7 elements are used. In symmetric cases (SS-SS and CL-CL cases), modes up to 4th are tested. In unsymmetric cases (CL and CL-SS cases), modes up to 3rd are tested. In Figure 10, the nodal location of the improved grids anf improved eigenvalue with respect to a single mode is illustrated in two cases, 88-88 and CL cases. Table 1 lists eigenvalues and relative errors for the uniform grid in those cases. Symmetric cases (SS-SS and CL-CL cases) need N = (22 +tfl elements for accurate estimates of error distributions, whereas unsymmetric cases (CL and CL-SS cases) need N = (32 + b) elements. In all cases except the 4th mode in 58-85 case, the similar reductions in both the exact error and in the indicator can be achieved. For example, suppose that after node relocations the maximum indicator is reduced by 50 %. Ifmniiwa can expect that the maximum exact error will be also reduced by about 50 %. 81 improved grid - - - - uniform grid 41b anode 88-88 4th mode 'el\\‘“ ifihflliol .0 Vine... 'udfla ill-.... 17 ... . I: ..-...n-HH r. a ’1‘ _ 'H‘ihiill r2 _ nuiuoa _ 111.6 r1. 4:43.940 illiiiiillliro m m m m m QOH x gouge oomxo J? I ‘6 1W . IIIIIII OIH llll m. lllllllllllll ....O 1. 'HH iiiiiii 17 it: 0. #0 ii L .5 OIHHHHI v‘ u l I lIlIlIIiih-‘O 13 ” ...... ....... ... . it 6 T... OH x ocumowch hOth element # element 5 3rd mode 3rd mode SS-SS 'C —i s|“t 1w IIIIII “|II Oils: 1. " " !!!!! ”He. 17 iiiii IIIII ‘III 1. . ."-' I liltililse fiC h iiiii . ..i\01 r) — onHHH rel ._- - ---. T 1'4 1 4i( 4 4 i Jl‘liiro n u a a a n o OH x nouooqpcfi uouuo element 4} element f 2nd mode 88-88 2nd mode _ --\. W \i‘ iiiiii r. er\I\o\\ \\ I litreaunv / r7 iiiiii T . iiiIVAH To \0 iiiiii r4 _ ......,- «x .. :15 111111 I/ ..2 _ ..... 1. i Tia-1.. .. - .5 oo~ x Louuo uoexe Il|i0 v” \\e iiiiii 1. Inn‘s 1. ---.-- T _ ...... a r. _ iiiie I: _ ...... .- . R--- r. . III . Id lllll Ia “ lllll 0 T.- .r -11 lilIIIIlier u u u u N o“ x woodcuts“ uouuo element # element- { Figure 9. Error distribution of uniform and improved grid 82 4th mode 83—58 4th mode SS—SS s, l I '- I \ I l C 10 ‘1 1'2 13 fillllllil lll‘S 2 ii iii 11 ’NHHI 1” tile :40. . ilie II OIHHIII '0 I III. v7 _ Ilia x6 .YHHII ..5 IIIII e v4 . ill. VJ 04.111111: 2 II 'I‘ 1| 111144441111 0 s o 5 c O Y] J mod 7. going—or; Congo element # element # 4th mode SS-SS 4th mode SS—SS V. W m qu x uouuo uomXo wmw N CH X HOUQOHQCH .AOHHO 1° :9 \\ \\ on... ‘8 JUN. LI l‘l I|||| OHHH 1‘ llllll l. 15 'nnli i4 IHUI. l3 enlii 2 I, III I a 4 A! element # element. # 3rd mode CL 3rd mode .~O._.~.u u Umwxw .: x ooH x uoumowch uonuo 9.1!!! Age-W III lllilili.‘ 1' V. lllll ... \ Fl!!! 17 I Ill.“ 1‘ e In I“. Is \9 |||||||| v4 \ a. llllll 1.. III!!! 1 4 4 4 AV) m w w m element # element # Figure 9. (cont'd) 2nd mode CL 2nd mode CL V w H o- uuuuu g. Tm T I . T \3\J r7 \i\\\ r. \ O\I f3 0 III .f/ : 74 11‘, as It) r2 41 4 4 4 a I: 2 .- C J 5 ooH x .aoumoHUCH uouuo element # element. # 1st mode CL— CL 1st mode CL-CL OH X HOHHG UUNXQ element # N OH X HOufiowUCH Mowkm element # Figure 9. (cont'd) 84, 3rd mode CL- SS 3rd mode CL—SS I I .li‘i‘ Illiu l li‘iilIIIIIII n ‘1}... ‘14. i - ‘ II- " |‘-|‘ I! " “ ' {T1311 m m m m .. OOH X HOHHm uUme H III:IIe 1” OflHHH iiiii .v T. IIIIII .701 T I 1:2». .7 ‘HMMH llll c #8 IIIIII 47/ T3 IIATIIHHJI 14 8w“ IIIIIIII A. f} I If. ..... : f2 A ------- .1 (i.. 44111111141110 CH x mounowpcw uouuo G element 5 element # 2nd mode CL-SS 2nd mode CL—SS ‘ | O 6 ion 5 l l l i l l l l I I I 6 IL»..- T. r4 1101!! r3 iiiiiiii x .. ‘ ( T1 JIiIi 4 4+4 1*14 {lira w w u w m o ‘I 00H X HOHHG UUQXG OH x HOUmoHUCw nouuo element- if element it Figure 9. (cont'd) 85 {whxnode CL-CL 4th anode CL—CL ed 2 nonuo uooxo o~ x |||i0 1. - “““ in. “ “ ‘l’ ’80 Ill”. 13 ‘| W 1. v. H. vs I‘ Leumowozfi QOHHU element { element 5 3rdlnode CL-CL 3rdxnode CL—CL Vi““‘¢l‘l‘1¢{‘ li‘.‘“‘ii‘ oo~ x uouuo uooxo OH x aoumoaoca uouuo element 5 element 5 2ndxnode CL-CL 2ndlnode CL-CL ' 4' . l-ml ml--.- oo~ x uouuo uomxo H IIIII .6 r7 1 ........ + T. u Fain: I: II I... lllllll ..z r‘ e iiiiiiiiii .3 a .- ----- + r. H IIIIIII 6 r1 m-i-m-l-m--I-mIl- Nod x wouooHoCL nouuo element I element 5 Figure 9. (cont'd) 86 pined-pined (SS-SS) uniform grid - A ; A e - v 2nd mode optimized grid A2 - 1558.79439 3rd mode optimized grid A3 - 7895.98101 A a A A v v A f 4th mode optimized grid A“ f 24996.22791 cantilever (CL) 2nd mode optimized grid A2 - 485.53422 3rd mode optimized grid A3 - 3807.76813 Figure 10. Nodal location of the improved grid PINED-PINED (SS-SS) CASE NUMBER OF ELEMENTS - CANTILEVER (CL) CASE NUMBER OF ELEMENTS- ClAMPED-PINED (CL-SS) CASE 10 10 NUMBER OF ELEMENTS- 10 CLAPMED-CLAMPED (CLvCL) CASE NUMBER OF ELEMENTS- 7 87 MODE PLO-3N MODE N MODE WNH WNW EIGENVALUE 1558.87662 7898-54283 25019-26421 EIGENVALUE 485.54857 3808.46023 EIGENVALUE 237-72772 2497.32518 10883.44947 EIGENVALUE 500.69958 3811.58909 14733.62586 REL ERROR(§) 0.0212481848 0.1065438590 0.3309853171 REL ERROR(§) 0.0061282991 0.0502809047 REL ERROR(\) 0.0027999471 0.0335568255 0.1460053260 REL ERROR(‘) 0.0271058632 0.2116979730 0.7935330701 Table 1. Eigenvalues and relative errors for the uniform grid CHAPTER V. CONCLUSIONS AND RECOMMENDATIONS 5.1 Conclusions We established error estimates for elliptic eigenvalue problems for each mode based on interpolation error theory. The error indicator was implemented using two different techniques, the smoothing technique and the direct substitution technique. The accuracy of error indicators was checked using the eigenvalue problem of a uniform, Bernoulli-Euler beam. Based on the accuracy on maximum error, the following conclusions were drawn : (l) The maximum error of the i-th mode eigenfunction is the same for all modes 2 when the number of elements N - (22 + b), (32 + b), etc. is used. The constant b depends on the boundary conditions. It equals 0 in the 85-85 case, -1 in the CL case, and l in the CL-SS and the CL-CL cases. This is true for both cubic and quintic elements and is true for the maximum error estimator and maximum error indicators. (2) Both error indicators can estimate accurately the maximum 88 error in the I-th mode eigenfunction when more than N=- (32 + b) elements are used. Based on the accuracy on error distributions. the following conclusions were drawn: (1) The interp type 2 indicator can estimate more accurately the distribution of the true error than the interpolation type 1. When higher order elements are used. the interpolation type 2 estimator can estimate error distributions very accurately. (2) Symmetric boundary conditions have an important effect on the number of elements needed to estimate accurately the error distribution. Using cubic beam elements. the interpolation type 2 indicator can estimate accurately the distribution of the true error when N-(22 i 1) elements are used in symmetric cases and N-(32 t 1) elements are used in unsymmetric cases. Using these results,improved grids were constructed by node relocation. It was observed that the similar reductions in both the exact error and in the indicator were achieved. 5.2 Discussions and Future Reasearch The results from the improved grids (e.g.improved eigenvalues and eigenfunctions) may not be so impressive. The reason is that, in our model problem (eigenvalue problem of the uniform, Bernoulli-Euler 9O beam), the uniform grid itself can estimate lower mode eigenvalues with enough accuracy. An extension of l-D beam problems to 2-D, 4th order problems (e.g. vibration of plates) may be another candidate. Main issues in 2-D problems are somewhat different from those in l-D problems. In l-D problems, the issue is 'how many elements’ and 'what type of element' are needed to estimate eigenfunctions with enough accuracy. In 2-D problems, however, an important issue is 'where we should put additional elements'. Accuracy estimates of eigenvalue problems may be applied to many engineering problems. As a familiar example of a stability problem from solid mechanics, the buckling of a column may be one candidate. Here, the buckling load is proportional to the smallest eigenvalue for the corresponding differential operator of beam theory and buckling mode is the associated eigenfunction. Accuracy estimates for l-D finite element eigenvalue problems may also be applied to modelling of beam-like structures for control design problems. In the development of feedback control theory for distributed parameter systems (DPS) which are described by partial differential equations, it is important to design an implementable finite dimensional controller. The implementation is usually done with on-line digital computers. The dimension of the controller is directly related to the memory capacity and the access time for retrieval of information from the computer memory. Therefore, the controller of the least possible dimension is preferred from a practical point of View. In order to design the finite dimensional 9l controller often a reduced order model (ROM) is sought. If the actual modes (eigenfuntions) are known, the dimension of the DPS may be reduced using projection into a modal space. In general, however. these modes are never known exactly and some other reasonable approximation procedure must be used. [One popular choice is the finite element method [35,36]. The feedback control force or moment is computed based on approximated mode shapes and mode slopes, and is applied to atabilize the motion of the DPS. The accuracy of approximated.solutions affect the stabilty of the DPS. To design the proper ROM, following two issues may be resolved : 1. Model selection problem, i.e., 'how many' and 'which' modes should be retained in the ROM. A modal cost analysis [38,39] may resolve this issue. 2. Finite element modeling problem, i.e., 'what type' auui "how many' elements should be used in the finite element ROM to produce stability of the original DPS. Accuracy estimates may resolve this. [10] [11] LIST OF REFERENCES Babuska, I. and Rheinbolt, V.C.,"Adaptive Approaches and Reliability Estimations in Finite Element Analysis" Comp. Meth. Appl. Mech. Engrg. (1979) pp 519-540 ,"A Posteriori Error Estimates for the Finite Element Method" Int. J. Num. Meth. Engrg. 12 (1978) pp 1597~1615 ,"Error Estimates for Adaptive Finite Element Computations" SIAM J. Numer. Anal. 15(4) (1978) pp 736-754 ,"Reliable Error Estimation and Mesh Adaptation for the Finite Element Method" in Oden, J.T. ed., Computational Methods in Nonlinear Mechanics pp 67-108 North- Holland Publishing Company (1980) I Daflandcz, L.,De1voo, Ph. and Oden, J.T.,"On an h-type Mesh- Refinement Strategy Based on Minimization of Interpolation Errors" Comp. Meth. Appl. Mech. Engrg. 53 (1985) pp 67-89 Babuska, I. and Miller, A.,"The Post-Processing Approach in the Finite Element Method" Part I: Calculation of Displacements, Stresses and.0ther Higher Derivatives of the Displacements Part II: The Calculation of Stress Intensity Factor Int. J. Num. Meth. Engrg. 20 (1984) pp 1085-1121 Babuska, I.,Szabo, B.A. and Katz, I.N.,"The p-version of the FEM" SIAM J. Numer. Anal. 18 (1981) pp 515-545 Zienkiewicz, O.C.,Gago, J. and Kelly, D.W.,"The Hierarchical Concept in Finite Element Analysis" Computers and Structures 16 (1983) pp 53-65 Zienkiewicz, 0.6.,Ke11y, D.W.,Gago, J. and Babuska,I., "Hierarchical Finite Element Approaches, Adaptive Refinement and Error Estimates" Proc. MAFELAP 1981 Zienkievicz, 0.0. and Craig, A., "Adaptive Refinement, Error Estimates, Multigrid Solution, and Hierarchic Finite Element Method Concepts " in Accuracy Estimates and Adaptive Refinements in Finite Element Computations (1986), John Wiley & Sons Ltd. Szabo, B.A~, "Estimation and Control of Error Based on p- convergence" in Accuracy Estimates and Adaptive Refinements in Finite Element Computations (1986), John Wiley & Sons Ltd. 92 [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] 93 Babuska, I. and Szabo, B.,"On the Rate of Convergence of FEM" Int. J. Num. Meth. Engrg. 18 (1982) pp 323-341 Babuska, I. and Dorr, ILR. ,"Error Estimates for the Combined h and p Versions of the FEM" Numer. Math. 37 (1981) pp 257-277 Diaz, A.R. ,Kikuchi, N. and Taylor, J.E.,"A Method of Grid Optimization for FE Methods" Comp. Meth. Appl. bkufli. Engrg. 41 (1983) pp 29-45 ' ,"Optimal Design Formulations for Finite Element Grid Adaptation" in Komkov V. eds, Sensitivity of Functionals with Application to Engineering Science, Lecture Notes in Mathematics 1086 pp 56-76 Diaz, A.R.,Kikuchi, N.,Papa1anbros, P. and Taylor, J.E.,"Design of an Optimal Grid for FE Methods" J. Struct. Mech. 11(2) (1983) ' pp 215-230 . Kikuchi, N.,"Adaptive Grid Design Methods for FE Analysis" Comp. Meth. Appl. Mech. Engrg. 55 (1986) pp 129-160 Demkowicz, L. and Oden, J.T. ,"On a Mesh Optimization Method based on a Minimization of Interpolation Error" Iru:..J. Engrg. Sci. 24(1) (1986) pp 55-68 Shephard, M.S.,Yerry, ILA. and Bachnann, P.L.,"Automatic Mesh Generation Allowing for Efficient a prior and a posteriori Mesh Refinement" Comp. Meth. Appl. Mech. Engrg. 55 (1986) pp 161-180 Bennighof, J.K. and Heirovitch, L. , "Eigenvalue Convergence in the Finite Element Method" Int. J. Num. Meth. in Engrg. 23 (1986) pp 2153-2165 Friberg, P.O.,"An Error Indicator for the Generalized Eigenvalue Problem using Hierarchical FEM” Int. J. Num. Meth. in Engrg. 23 (1986) pp 91-98 Priberg, P.O. and.Hb11er, P. Hakovicka, D. and Wieberg, N.E., "An Adaptive Procedure for Eigenvalue Problems using Hierarchical FEM” 24 (1987) pp 319-335 Killer, K. and Miller, R.,"Moving Finite Elements Part I" SIAM J. Numer. Anal. (1981) Miller, K.,"Moving Finite Elements Part II" SIAM J. Numer. Anal. (1981) pp 1019-1057 Miller, K.,"Recent Results on FEM with Moving Nodes" in Accuracy [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] Estimates and Adaptive Refinements in Finite Element Computations (1986), John Wiley & Sons Ltd. Adjerid, S. and Flaherty, J.E.,"A Moving Mesh FEM with Local Refinement for Parabolic PDEs" Comp. Meth. Appl. Mech. Engrg. 55 (1986) pp 3-26 Oden, J.T.,Denkowicz, L.,Strouboulis, T. and Delvoo, P.,"Adaptive Methods for Problems in Solid and Fluid Mechanics" in Accuracy Estimates and Adaptive Refinements in Finite Element Computations (1986), John Wiley & Sons Ltd. Demkowicz, L. and Oden, J.T.,"An Adaptive Characteristic Petrov- Galerkin FEM for Convection-Dominated Linear and Nonlinear Parabolic Problems in two space variables" Comp. Meth. Appl. Mech. Engrg. 59 (1986) pp 63-87 Oden, J.T.,Strouboulis, T. and Delvoo, P.,"Adaptive FEM for the analysis of Invscid Compressible Flow : Part I,Fast Refinement/Unrefinement and Moving Mesh Methods for Unstructured. Meshes" Comp. Meth. Appl. Mech. Engrg. 59 (1986) pp 327-362 Strouboulis, T. and Delvoo, P. and Oden,J.T., "A Moving-Grid Finite Element Algorithm for Supersonic Flow Interaction between Moving Bodies" Comp. Meth. Appl. Mech. Engrg. 59 (1986) pp 235- 255 Oden, J.T.,Stroubou1is, T.,De1voo, Ph. and Hove, M.,"Recent Advances in Error Estimation and Adaptive Improvement of FE Calculation" Hughes, T.J.R. and Hulbert, G.M., "Space-Time Finite Element Methods for Elastodynamics : Formulations and Error Estimates" Comp. Meth. Appl. Mech. Engrg. 66 (1988) pp 339-363 Babuska, I.,Zienk1evicz, O.C.,Gago, J. and Oliveira, A. eds., Accuracy Estimates and Adaptive Refinements in Finite Element Comptations (1986), John Wiley & Sons Ltd. Babuska, I.,Chandra, J. and Flaherty, J. eds., Adaptive Computational Methods for PDEs (1983) Balas, H.J. ,"The Galerkin Method and Feedback Control of Linear Distributed Parameter Systems" J. Math. Anal. Appl. 91 (1983) pp 527-546 ,"Toward a practical control theory for distributed parameter systems in control and dynamic systems : Advances in Theory and Applications" Leondes ed. vol.18 Academic Press (1982) pp 361-421 [37] [38] [39] [40] [41] [42] [43] [44] [45] [47] [48] [49] 95 Hughes, P.C.,"Space Structure Vibration Modes : How many exist? Which modes are important?" IEEE Control Systems Magazine (1987) Hu, A. and Skelton, R.E. ,"Convergence Properties of Modal Costs for Certain Distributed Parameter Systems" ASME Vibration Conference Proc., Boston (1987) Hu, A.,Ske1ton, R.E. and Yang, T.,"Modeling and Control of Beam- like Structure" J. Sound and Vibration (1987) Diaz, A.R., Jayasuria, S. and Kim, Y.J., "On the design of Nonlinear Controllers for Distributed Parameter Systems using Adaptive Finite Elements" Midwest Mechanics Conference (1987) Kikuchi, Noboru, Finite Element Methods in Mechanics Cambridge Univ Press (1986) The Texas Finite Element Series Prentice Hall Vol 1 Finite Elements : An Introduction Vol 2 Finite Elements : A Second Course Vol 3 Finite Elements : Computational Aspects Hughes, T.J.R., The Finite Element Method : Linear Static and Dynamic Finite Element Analysis Prentice-Hall (1987) Strang, G. and Fix, 6., An Analysis of the Finite Element Method Prentice-Hall (1973) Ciarlet, P.G., The Finite Element Method for Elliptic Problems North-Holland (1978) Melosh, R.J., "Finite Element Approximations in Transient Analysis" The Journal of the Astronautical Sciences, Vol XXXI, No. 3, PP 343-358 (1983) Sun, C.T. and Huang, S.N., "Transverse Input Problems by Higher order Beam Element" Computers and Structures 5, pp 297-303 (1975) Meirovitch, Leonard , Analytical Methods in Vibration The Macmillan Company (1967) APPENDICES Appendix A . Mathematical Background [42,43] A.1 Real Linear Space A real linear space is a collection of objects for which the operations of addition and scalar multiplication are defined and behave as follows: if x and y are members of the linear space and a and B are scalars, then ax + fly is also a member of the linear space. Linear spaces may be endowed with important structures, most significantly, inner products and norms. Definitions An inner product (.,.) on a real linear space A is a map (.,.) : AxA 4 R with the following properties: Let x,y,zEA and a E R, then 1. (x,y)-(y,x) (symmetry) 2. (x+y)- (y+x) (linearity) 3. (x+y,z)-(x+z)+(y+z) 4. (x,x) 2 0 (positive definiteness) (x,x)-0 if and only if x-O 96 A QQLE ||. || on a linear space A is a map || .||: A + m, with the following properties: Let x,yEA and aE R, then 1. ||x|| Z O and [|x|| - 0 if and only if x - 0 (positive definiteness) 2- IIOIXII ' lal ||X|| 3. [|x+y|| S ||x|| + ||y|| (triangular inequality) A semi-norm |.| on a linear space A is a map |.|:A * R with the following properties Let x,y 6 A and aE R, then 1. |x] 2 0 (positive semi-definiteness) 2. [0er ' lallxl 3. |x+y| 5 |x| + |y] (triangular inequality) A.2 The Continuity Class Cm(0) Suppose that 0 is a bounded region in 'R% and that u is a given real-valued function of position in 0. Then u is said to be the class of Cm(0) on 0 if u and all of its partial derivatives up to order m (m is nonnegative integer) are continuous at every point in O, i e., r Cm(0)-{ v | —2—%———-— is continuous in Q C.R21 6x 6yj i,j z 0 , i+j - r , r s m The class Cm(0) is a linear space of functions. 98 A.3 The L2(0) Class . . . . 2 . . . The functuxni f is said to be in L class if its derivatives can be defined onLy in a square integral sense. The function f in L2 class will have the property of f0 fzdn < w . A.4 The Sobolev Class Hm(0) A function u is in Hmm) if u and all of its partial derivatives of up to order m (nonnegative integer) are members of L2(0), i.e. , m m av 8v 3 v '2 Y The space Hm(0) is a Hilbert space. The Hm(0) inner product, norm, and semi-norm are defined by 8au aflv 2 h m . (u,v)m- ( In a§BSm[ a xa a yfl 1 dx } (H inner product) 8 m [lullm' (U,u)m (H norm) (a+fi) a u 2 dX )5 (Hm semi-norm) lulm- { f0 agfl-m [ a xa a yfi ] where a and 6 are nonnegative integers. A.5 Equivalence of Two Norms 99 (l) I(2) Two norms, [ |.]| and |].| , on a linear space A, are said c such that to be equivalent if there exist constants c 2 l' (l) IXI|(1)S lenm for all x A cll SCZIIXII 100 Appendix B . Interpolation Functions for Quintic Beam Elements 3.1 Quintic Element Type 1 (3 nodes, 2 dof/node) We list here interpolation functions for the quintic beam element with 3 nodes and 2 degrees of freedom per each node. A model of N quintic elements has a total 2(2N+1) degrees of freedom. This element has continuity in displacements and slopes (lst derivative) over the domain. Moments (2nd derivative) and shear forces (3rd derivative) are, however, discontinuous at element boundaries. The 5th order polynomial interpolation functions of this quintic elementcuuibe constructed using continuity conditions of displacement and slope at: three nodes. Interpolation functions for this element are as follows 1 - 23 n2 + 66 "3. 68 "4+ 24 "5 N1(n) 2 3 4 5 N2(n) - LK( n - 6 n + 13 n - 12 n + 4 n ) N3(n) 16 "2 - 32 "3. 16 n“ 2 3 4 5 N4(n) - LK( - 8 n + 32 n - 40 n + 16 n ) N5(n) - 7 "2 - 36 "3+ 52 n“- 24 "5 101 2 3 4 N6(n) - LK( - n + 5 n - 8 n + 4 n ). U! 3.2 Quintic Element Type 2 (2 nodes, 3 def/node) Another type of quintic element can also be defined using 2 nodes and 3 degrees of freedom per node. A model of N quintic elements has a total 3(N+l) degrees of freedom. The difference between the quintic element type 1 and type 2 is the following: lkithe quintic element type 1, the displacement and slope (up to lst derivative) are continuous across elements, and the moment (curvature, 2nd derivative) is discontinuous at element boundaries. However, in the quintic element type 2, the moment (curvature, up to 2nd derivatives) is also continuous at element boundaries. This type of quintic element has proven to be very efficient in the dynamical analysis of beam-like structures [48]. Interpolation functions for this type of element are as follows: N1(n) - 1 - 10 n3 + 15 "a - 6 "5 3 4 5 N2(n) - LK( n - 6 n + 8 n - 3 n ) 2 2 3 3 5 N3(n) - LK( 0.5 n - 1.5 n + 1.5 n - 0.5 n ) 10 n3 - 15 "4 + 6 n5 Na(n) N5(n) - LK< -4 n3 + 7 n“ - 3 n5) 102 N6(n) - L§< 0.5 n - n“ + 0.5 as) 103 Appendix C . Exact Eigenpairs of a Uniform, Bernoulli-Euler Bean The equation of motion of a uniform, Bernoulli-Euler beam of length 2 can be written as 84(xt) 62(xt) EI 3" +m y' -o t>0,0