83_

Coraputt'rs & Structure~ Voi 16. Nu I--4. pp. 53--65, 1983 0045-7949!83{010053--13503.00/0 Printed in Great Britain Pergamon Press Ltd THE HIERARCHICAL CONCEPT IN FINITE ELEMENT ANALYSIS O. C. ZIENKIEWlCZ, J, P. DE S. R. GAGO and D. W. KELLY Department of Civil Engineering, University College of Swansea, Singleton Park, Swansea SA2 8PP, Wales Abstract--The hierarchical concept for finite element shap e functions was introduced many years ago as a convenient device for mixed order interpolation. Its full advantages have not been realized until a much later time--and these include in addition (a) improved conditioning; (b) ease of introducing error indicators if successive refinement is sought. Further, it is possible to use the ideas to construct a range of error estimators which compare well with alternatives and are ideally suited for adaptive refinements of analysis. As the hierarchical elements are equally simple to implement as "standard" fixed order elements it is felt that more programs will. in the furore, turn to utilize their advantages. This is especially true in the field of nonlinear analysis where even today computational economies are necessary. INTRODUCTION Despite the wide usage of finite elements today little attention has been given to the problem of estimating errors occurring in the analysis and to correcting these errors by automatic refinement. These aspects of the problems are of obvious importance if intelligent use of programs is to be made in engineering practice. In typical application today the experienced analyst decides on the mesh of elements to be used and generates the information necessary automatically[l-3]. In distributing his nodes he is guided by (1) Previous solutions which gave satisfactory ans- wers when tested against benchmark problems. (2) Guidelines generated by previous investigations into mesh "optimality" to which much attention has been given in the past [4--8]. When presented with the computer output a quick check on some obvious errors such as static force balances, stress discontinuities between elements reas- sures him that all is well and he then proceeds to draw necessary design conclusions. Generally if all is NOT well, a new mesh is generated from scratch and the problem resolved with much of the costly human inter- vention duplicated. In regenerating the mesh the user has often the choice of saving some data preparation efforts by retaining the original mesh and refining locally, either by (a) Introducing new elements of the type used ori- ginally but of smaller size (h) or (b) Using same element definitions but increasing the order of polynomials used (p) involving new modes placed on such elements. In Fig. 1 we illustrate the two possibilities, often called h or p refinements respectively. Once the new subdivision is accomplished a com- pletely new solution is generally performed and con- vergence observed. The computer at this stage repeats all the computations and makes no use of the previously generated solution. The reader can contrast this procedure with a more classical Ritz (or Galerkin) approach in which a generally simple geometry may be solved utilizing a trigonometric Fourier series expansion. Here as new terms are called the previous terms remain unaltered and a single term is introduced at a time, the refinement ceasing when the accuracy of solution is deemed sufficient. Why the difference between the two approaches when we are all aware that both the finite element process AND the series solutions are mathematically precisely of / / I OriginoL mesh h- Refinemen+~ (subdivision) 53 _2/f p-Refinement (increose of polynomiaL, order) Fig. 1. Possible refinements of an inaccurate mesh. 54 ~). C. ZiENK,.'EWiCZ et aL the same kind? The answer lies in the fact that generally the sequence of meshes derived in subdivision is NOT hierarchical while the Fourier series expansion is. By hierarchical we mean the following: if u is the unknown function (e.g. the structural displacement) and its ap- proximation is u=~ =~] Nia, (l) then approximation is hierarchical if an increase of n to n + 1 does not alter the shape functions N~(i = l to n). Clearly this is not the case in a standard finite element formulation where a subdivision of an element is ac- companied by introducing entirely new functions (Fig. 2a). This could of course be remedied by introducing the new approximation hierarchically as shown in Fig. 2(b), but here the conventional concept of elements, each formulated by an immutable procedure would have to be partially abandoned. The only "element" in such a form would remain as the original one. Such an idea is of course possible although rather special integration rules would have to be designed for the element to allow for the interaction of "senior" hierarchical functions (such as Nt) with "juniors" (typi- cal of N~). Hierarchical functions are introduced more simply in a "p" type context and here their simplicity is evident. In Fig. 3 we show for instance the generation of hierar- chical and non-hierarchical functions for a quadratic isoparametric element. It was indeed in this environment that hierarchical functions were initially introduced by Zienkiewicz et aL[9] as early as 1971. At that time the objective was simply the introduction of p-graded meshes in an a priori chosen manner. Since that time some use of the hierarchical subdivision has taken place[10] Fig. 4, and indeed new and useful families of hierarchical p type elements were introduced by Peano et aL[l 1-13]. We discuss some such functions in more detail elsewhere [20], but would like to stress that the hierar- 3 Ns N~ (o1 Standard [~i Nier~rc~J¢ Fig. 3. Increase of polynomials in element (p-refinement). Stan- dard (a) and hierarchical (b) shape functions. chical forms have other advantages additional to the ease of introducing local refinements and that their use could well become much more widespread than it is at present. 7. mmrrs oF ~ mmumcmc~ coNcv.vr 2.1 Structure o/equations Let us consider a typical linear problem for which a discrete approximation is sought. Such a problem could well be the solution of the elasticity equations in a typical stress analysis situation but for generality we shall state it is satisfying a set of differential equations and boundary conditions Lu + q = 0 in fl (domain) Mn + S = 0 on F (boundary). (2) 8 6 5 N6 Ns In the above L and M are linear differential operators. When an approximate solution is sought in the form n = ~ N~a~ = Na t"' (3) i=t where a (") are undetermined parameters and a Galerkin type weighting is used in a standard approach we obtain a set of algebraic equations /n NrL(Na~"') dll + fn Nrq dfI + b = 0 (4) where b is dependent on boundary conditions, or (on integration by parts) ( a ) S tandard ( b ) H ierarch ic Fig. 2. Subdivision of an element (h-refinement) standard (a) and hierarchical (b) shape functions. I~.@" = f"' ~5) where I~,) is an appropriate stiffness matrix involving the derivates of the shape functions N and f is the "force" vector. Details of the computation are well known and need not be repeated here [14]. The hierarchical concept in finite element analysis DE GENfRA'rE ~ l ELEMEN1 la~ h~ERA~J, eC.~ E~.EME~ 55 Fig. 4. Three dimensional analysis of a pressure vessel indicating advantages in idealization provided by hierarchical elements. When the mesh is "subdivided", or new degrees of freedom introduced, the number of parameters is in- creased to a ~'+m~, and identical discretization process results in algebraic equations If the refinement is made hierarchically the original stiffness coefficients reappear and eqn (6) could be writ- ten as KI,.,,~/~ i, = (7) / LK ... . . g,,,, J I.a"~'J [ f " ' J Immediately we observe the first merit, i.e. that the matrix K,,~ is preserved, possibly saving some coefficient computations on refinement. This by itself may be significant if elements are complex but more important perhaps is the possible use of the original solution of the unrefined mesh {eqn 5) which we could denote as a-,~. = (Kl,,~)-'.f I'> to start the next solution. Peano[19] discusses the pos- sibility of storing the triangularized form of (K,,>) -~ in attempting solution of eqn (7) and such a procedure has much to commend it. More importantly however a~,,* could be used as a start of an iterative solution giving the first approximation for a '"> as a I~'* =K. , ,~( f -K l~ . , .a ). 2.2 Equation conditioning With the hierarchical degrees of freedom appearing as perturbations on the original solution rather than its substitute we would expect eqn (7) to have a more dominantly diagonal form than that obtainable in a direct (6) approximation involving the identical number of non- hierarchic degrees of freedom. This has important con- sequences of ensuring an improved conditioning of the matrix K¢,+,,~ of eqn (6) and a faster rate of iteration convergence than would be possible with non-hierar- chical forms. In Fig. 5 we show in detail this improvement in con- ditioning on the basis of a single cubic element and an assembly of four such elements. In both cases the con- dition number of the stiffness matrix is improved by an order of magnitude and in larger assemblies yet more dramatic improvements may be expected. This clearly is of importance in the equation solution if a limited number of digits is carried. With the increasing use of micro-computers this aspect alone could well "sell" the concept to many users. In Table 1 and Fig. 6 (8) we show a very simple example of a slender cantilever and the roundoff errors occurring with different slender- ness ratios for both hierarchical and non-hierarchical elements which being of the same order should give precisely identical results. At this point it is useful to digress and mention that precisely this point of conditioning was advanced early by others. Wilson in the U.S. and Japan Seminar of 1975 and elsewhere[15, 16] has addressed the problem show- (9i ing the usefulness of introducing variables involving ~6 (). (7. ZIENK;EWICZ et aL {:7 :]. ~!07 : : - Ic ELEMENTWISE I :lbBIC ELEMENT f 25 0 25 0 i o AFTER NORMALISATION OF K C0nd K'~ = .Lmc,,, = 390 27 ~m,n Cond Km' = 36 28 PATCH OF ~I~MF,NTS ~ CUBlC ELEMENTS; / I>~ z ~ : ,, Ill J 25 0 21 666 Cond K~%I6t~3 ~2 Cond K m'=12378 ~ tC> ,:=, 25 0 21 666 Fig. 5. Improvement of conditioning associated with hierarchical formulation. ~ame hierarchical structure and characteristics ~Flg. ~). A hierarchical form was introduced ti~o by Wachspress[18] in which local element subdivision can easily be achieved. 2.3 Error measures The perturbation nature of hierarchical forms has a further merit of providing an immediate estimate of the error in the solution. This error can be defined as e = u-Li i 0) where u and 6 are the exact and approximate solutions respectively. If two successive solutions are used with n and (n + m) parameters this error can be estimated as the difference between successive solutions e=tV"~' -6 .... Ill) and this indeed is precisely what an intelligent analyst would generally do. However, now we shall introduce a possibility of estimating the error without the recourse to the complete second solution. If fi(~) is determined as the first solution (viz. eqn (8)) @~"' = Nl"'a '"'* (12) then the first approximation to the error can be obtained as the first iteration for a ('~ given by eqn (9). Thus N" 'a " ) * ~N")K -j/f"~ K a"* ) e ~ - ( , .~ - . , , . . ) . (13) This approximation is of crucial importance to our further discussions. _L L i i i i , i _ t= 15.0 h Fig. 6. Cantilever beam for conditioning studies. displacement differences in thin beam or joint problems. Th is indeed can be cast as a special case of hierarchical form illustrated in Fig. 7. Indeed the concept of global and local approximations introduced by Mote[H] has the 2.4 Summary In the preceding we have noted that the hierarchical forms have the merits of (1) Utilizing previous solutions and computation when attempting a refinement: (2) Permitting a simple iteration in solving for suc- cessive refinements; (3) Always resulting in improved equation condition- ing; and (4) Giving an immediate error measure. While items (1) and (2) are only relevant if successive refinements are to be attempted (3) and (4) are present even if a single solution only is to be used. Table 1.1. Pure traction example on cantilever beam with varying thickness t. Stress evaluation at point (7.7.0.50 ~,,~, = 1.0000, tip displacement 8oxac, -- 15.0 I £ i! I" t A B A B 1, 0.01000 0.0% 0.0% 0.0% 0.0% I I 0.00100 1.6% 0.1% I 0.27, 0.0% O. 00010 15.7% 7.0% ] 17.4% O. 5% 0. 00001 FALLS FALLS I. CASE A - I0 e igh~-noded isopara~etr ie etements CASE 8 - i0 e ight-noded h ierarchic e lomen:s CASE I - relative error in compu=ed Cxx s=ress CASE IT - relative error in computed d isp lacement The hierarchical concept in finite element analysis 57 Table 1.2. Applied shear at the free end with varying thickness: stress evaluation at point (8.25, 0.3873t) ,r+x~ct = 31.371t, tip displacement 6~x~c~ = 0.135 (theory of beams) t I.OO O. I0 0.01 I 0.O% 3.9% FALLS O. 0% 3.3% FAILS A 0.02% 4.58% FALLS B 0.02% 0.07% FALLS CASE A - iO e ight-noded isoparametr ic elements CASE B - IO eight-noded hierarchic elements CASE I - relative error in acomputed stesses CASE Ii - relative error in computed displacement h ' [ l i l l l l i l t l i l l l i l l l l~ l l l l l l l ] l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l~ ,/,///,///,/,-, ,:~ ~ , ,////:~ ~ "//: @j or~+A¢, Fig. 7. A joint element in a hierarchical form. N¢; , . . . . l~obo l function f / f f ~ Finite element suoerposeO Fig. 8. Local-global finite element method, a hierarchical form. What then are the drawbacks which prevent the widespread use. The answer to this is less clear. With subdivided elements in hierarchic form, the drawback of integration on subregions was already men- tioned. With p refinement the above objection disappears but another one remains. This concerns the traditional concept of the unknown parameters a corresponding to nodal displacements. Now these are the differences of such displacements or indeed less easily identifiable quantities. Both defects can however be readily remedied. Firstly special integration rules can be introduced. Secondly by a very simple transformation (implicit in eqn (12)) the parameters a can be replaced by displacements (or values of the unknown function) or specified points. We can not therefore see any reason for not using hierarchical variables in standard programs. CAS 16:1/4 - E 3. ERROR INDICATORS AND FA'TIMATORS 3.1 Definitions If we are to adopt the procedures of refinement two types of question will arise. First, we shall need an INDICATION WHERE additional refinement is most effective and second a MEASURE OF ERROR deter- mining whether the refinement is necessary. In the previous section we have defined the "error" by eqn (10) as e=u-fl. (14) Now this is a purely local definition and inconvenient if the whole problem is to be considered. Further this does not tell us how welt the original (equilibrium) equations such as eqn (2) are satisfied. For this reason we introduce a more global error measure defined as the energy of the error [ie[l~: = fn erLe dfl. 05) This measure is of considerable importance as it can be easily related to the residual (r) which tell us how well the original equations are satisfied. (We neglect here any error in boundary conditions which can be treated similarly.) This residual is obtained as r=Lf i+q. (16) Using eqn (14) we note that Le = Lu - Lfi 58 O.C. ZIENKIEWICZ et ai. and inserting eqn (2) and (16) we have r = -Le. Equation (15) can thus be rewritten I]ei]~" = - ( er r df~ Jfl giving one of the basic relationships useful in later analysis. 3.2 Error indicators Consider now a hierarchical introduction of a single DOF a . . , into an n degree of freedom system for which we have determined a solution a ~"). We shall attempt to find the approximate energy error i le°-,tl~: ~ '7~+, which would be CORRECTED by the introduction of the new variable and which will provide therefore an error INDICATOR. By previous arguments (viz. eqn 13) we have already succeeded in obtaining the approximation to the value of e.+l as e,,+) --- gn . la , ,+t and a . . l = (./,,+,- K.+I.,,#"))IK. +,..., (21) if a,÷, is a single scalar. Now we shall examine in more detail the expression in parenthesis of eqn (21). We note that this has been obtained by substituting ~ = Na t') into the governing equation and "weighting" this with N.+ z. Thus 1".. t - K,,+,.,,a(") = ( N . . , r dl). 111 3.3 Error estimation In principle error indicators evaluated for ali possibte (17) new degrees of freedom capable of introduction and summed should provide an accurate measure to the total energy error occurring in the solution at the refinement stage reached. Unfortunately this is impracticabie and only the NEXT degree of refinement is generally studied (18) If it so happens that the next shape function is ortho- gonal (or nearly so) to the residual (i.e. j'fl N,,_, r dl~ ~ 0). and the bulk of the error lies in a higher order refinement, then the error estimate will be poor--ai- though the value of expression (24) as an indicator continues telling us simply that little change of the error would occur by introducing the refinement in question. Much effort has therefore gone into the matter ,:)f est- ablishing reliable error estimators which will give a more accurate error value possibly exceeding the vaiue from above. Here classical work by Babuska et at. pioneered (19) suitable estimates for linear quadrilateral elements and more recently[20.26,27] more elaborate expressions have been proposed for higher order elements. In the present paper we shall introduce an alternative estimator based directly on the hierarchical concept and indeed on the indicator given by the expression I24a). The performance of this new estimator appears satis- factory and at all stages as shown in later examples, overestimates the true error by a reasonable margin. It (20) must be mentioned that rigorous proofs of its ultimate convergence to the true errors are still lacking and that it is now proposed heuristically. If we consider expression (24a), we note that the sign of both the shape functions and of the residual may well change within an element. However it is possible to say that We now can write using definitions (18) and (19) (22) or L rl.+, = - e . . i r dfZ. (23) Inserting (20)-(22) we can identify v/~,.t as O'.+l = N . . l rd ,,.-i,,,+, (24a) r l~,*, = 0c .+, - K . . , . , ,a ( " ) ) : /K .+, . .+~. (24b) While the form (24b) may appear more direct and has been/dent/fled as an energy error indicator by Peano[19] we find that the alternative form given by (24a) and introduced by the present authors [20, 271 has a special sign/ficance since we can avoid the computation of some stiffness terms of the additional solution. We stress once again that no additional solution is required to evaluate the indicators and this can be done by expression (24) separately for each possible new degree of freedom INDICATING where corrections are most desirable. ,2,, and that the quantity on the r.