FINITE ELEMENT ANALYSIS AS COMPUTATION

FINITE ELEMENT ANALYSIS AS COMPUTATION

What the textbooks don't teach you about finite element analysis Gangan Prathap

Preface
These course notes are specially prepared for instruction into the principles underlying good finite element analysis practice. You might ask: Why another book on FEA (Finite Element Analysis) and the FEM (Finite Element Method) when there are hundreds of good textbooks in the market and perhaps thousands of manuals, handbooks, course notes, documentation manuals and primers available on the subject. Some time earlier, I had written a very advanced monograph on the analytical and philosophical bases for the design and development of the finite elements which go into structural analysis software packages. Perhaps a hundred persons work on element development in the world and out of these, maybe a dozen would probe deep enough to want to explore the foundational bases for the method. This was therefore addressed at a very small audience. In the world at large are several hundred thousand engineers, technicians, teachers and students who routinely use these packages in design, analysis, teaching or study environments. Their needs are well served by the basic textbooks and the documentation manuals that accompany installations of FEM software. However, my experience with students and teachers and technicians and engineers over two decades of interaction is that many if not most are oblivious to the basic principles that drive the method. All of them can understand the superficial basis of the FEM - what goes into the packages; what comes out; how to interpret results and so on. But few could put a finger on why the method does what it does; this deeper examination of the basis of the method is available to a very few. The FEM is now formally over thirty-five years old (the terminology "finite element" having been coined in 1960). It is an approximate method of solving problems that arise in science and engineering. It in fact originated and grew as such, by intuition and inspired guess, by hard work and trial and error. Its origins can be traced to aeronautical and civil engineering practice, mainly from the point of view of structural engineering. Today, it can be used, with clever variations, to solve a wide variety of problems in science and engineering. There are billions of dollars worth of installed software and hardware dedicated to finite element analysis all over the world and perhaps billions of dollars are spent on analysis costs using this software every year. It is therefore a remarkable commercial success of human imagination, skill and craft. The science of structural mechanics is well known. However, FEA, being an approximate method has a uniqueness of its own and these curious features are taken for closer examination in this course. The initial chapters are short ones that examine how FEM is the art of numerical and computational analysis applied to problems arising from structural engineering science. The remaining chapters elucidate the first principles, which govern the discretisation process by which "exact theory becomes numerical approximations. We shall close these notes with a telescopic look into the present scenario and the shape of things to come in the future.

iii

Acknowledgements: I am extremely grateful to Dr T S Prahlad, Director, N.A.L. for his encouragement and support. Many colleagues were associated at various stages of writing this book. They include Dr B P Naganarayana, Dr S Rajendran and Mr Vinayak Ugru, Dr Sudhakar Marur and Dr Somenath Mukherjee. I owe a great deal to them. Mr Arjun provided secretarial help and also prepared the line drawings and I am thankful to him. I am also thankful to Latha and Rahul for their patience and understanding during the preparation of this manuscript.

Gangan Prathap Bangalore 2001

iv

Contents Preface 1. Introduction: From Science to Computation 2. Paradigms and some approximate solutions 3. Completeness and continuity: How to choose shape functions? 4. Convergence and Errors 5. The locking phenomena 6. Shear locking 7. Membrane locking, parasitic shear and incompressible locking 8. Stress consistency 9. Conclusion iii 1 11

28 33 45 57

79 97 111

Chapter 1 Introduction: From Science to Computation

1.1 From Art through Science to Computation
It will do very nicely to start with a description of what an engineer does, before we proceed to define briefly how and where structural mechanics plays a crucial role in the life of the engineer. The structural engineer's appointed office in society is to enclose space for activity and living, and sometimes does so giving time as well - the ship builder, the railway engineer and the aerospace engineer enable travel in enclosed spaces that provide safety with speed in travel. He did this first, by imitating the structural forms already present in Nature. From Art imitating Life, one is led to codification of the accumulated wisdom as science - the laws of mechanics, elasticity, theories of the various structural elements like beams, plates and shells etc. From Archimedes' use of the Principle of Virtual Work to derive the law of the lever, through Galileo and Hooke to Euler, Lagrange, Love, Kirchhoff, Rayleigh, etc. we see the theoretical and mathematical foundations being laid. These where then copiously used by engineers to fabricate structural forms for civil and military functions. Solid and Structural Mechanics is therefore the scientific basis for the design, testing, evaluation and certification of structural forms made from material bodies to ensure proper function, safety, reliability and efficiency. Today, analytical methods of solution, which are too restricted in application, have been replaced by computational schemes ideal for implementation on the digital computer. By far the most popular method in Computational Structural Mechanics is that called the Finite Element Method. People who use computational devices (hardware and software) with hardly any reference to experiment or theory now perform design and Analysis. From advertising claims made by major fem and cad software vendors (5000 sites or installations of one package, 32,000 seats of another industry standard vendor, etc.); it is possible to estimate the number of fem computationalists or analysts as lying between a hundred thousand and two hundred thousand. It is to them that this book is addressed.

1.2 Structural Mechanics

A structure is "any assemblage of materials which is intended to sustain loads". Every plant and animal, every artifact made by man or beast has to sustain greater or less forces without breaking. The examination of how the laws of nature operate in structural form, the theory of structural form, is what we shall call the field of structural or solid mechanics. Thus, the body of knowledge related to construction of structural form is collected, collated, refined, tested and verified to emerge as a science. We can therefore think of solid and structural mechanics as the scientific basis for the design, testing, evaluation and certification of structural forms made from material bodies to ensure proper function, safety, reliability and efficiency. To understand the success of structural mechanics, one must understand its place in the larger canvas of classical mechanics and mathematical analysis,

especially the triumph of the infinitesimal calculus of continuum behavior in physics. The development and application of complex mathematical tools led to the growth of the branch of mathematical physics. This therefore encouraged the study of the properties of elastic materials and of elasticity - the description of deformation and internal forces in an elastic body under the action of external forces using the same mathematical equipment that was been used in other classical branches of physics. Thus, from great mathematicians such as Cauchy, Navier, Poisson, Lagrange, Euler, Sophie Germain, came the formulation of basic governing differential equations which originated from the application of the infinitesimal calculus to the behavior of structural bodies. The study is thus complete only if the solutions to such equations could be found. For nearly a century, these solutions were made by analytical techniques however, these were possible only for a very few situations where by clever conspiracy, the loads and geometries were so simplified that the problem became tractable. However, by ingenuity, the engineer could use this limited library of simple solutions to construct meaningful pictures of more complicated situations. It is the role of the structural designer to ensure that the artifacts he designs serve the desired functions with maximum efficiency, but putting only as much flesh as is needed on the bare skeleton. Often, this is measured in terms of economy of cost and/or weight. Thus for vehicles,low weight is of the essence, as it relates directly to speed of movement and cost of operation. The efficiency of such a structure will depend on how every bit of material that is used in components and joints is put to maximum stress without failing. Strength is therefore the primary driver of design. Other design drivers are stiffness (parts of the structures must not have excessive deflections), buckling (members should not deflect catastrophically under certain loading conditions), etc. It is not possible to do this check for structural integrity without having sophisticated tools for analysis. FEM packages are therefore of crucial relevance here. In fact, the modern trend is to integrate fem analysis tools with solid modeling and CAD/CAM software in a single seamless chain or cycle all the way from concept and design to preparation of tooling instructions for manufacture using numerically controlled machines.

1.3 From analytical solutions to computational procedures 1.3.1 Introduction

We have seen earlier in this chapter how a body of knowledge that governs the behavior of materials, solids and structures came to exist. Gifted mathematicians and engineers were able to formulate precise laws that governed the behavior of such systems and could apply these laws through mathematical models that described the behavior accurately. As these mathematical models had to take into account the continuum nature of the structural bodies, the description often resulted in what are called differential equations. Depending on the sophistication of the description, these differential equations could be very complex. In a few cases, one could simplify the behavior to gross relationships - for a bar, or a spring, it was possible to replace the differential equation description with a simple relation between forces and displacements; we shall in fact see that this is a very simple and direct act of the discretization process that is the very essence of the finite element process. For a more general continuum structure, such discretizations or simplifications were not obvious at first. Therefore, there was no recourse

Fig. 1.1 A simple bar problem under axial load

except to find by mathematical and analytical techniques, viable solutions to the governing differential equations. This can become very formidable in cases of complex geometry, loading and material behavior, and often were intractable; however with ingenuity, solutions could be worked out for a useful variety of structural forms. Thus a vast body of solutions to the elastic, plastic, viscoelastic behavior of bars, beams, plates and shells has been built over the last two centuries. It was recognized quite early that where analytical techniques fail or were difficult to implement, approximate techniques could be devised to compute the answers. In fact, much before algebra and analysis arrived, in contexts more general than solid or structural mechanics, simple approximation and computational schemes were used to solve engineering and scientific problems. In this chapter, we shall review what we mean by an analytical solution for a very simple problem and then proceed to examine some computational solutions. This will reveal to us how the discretization process is logically set up for such a problem.

1.3.2 The simple bar problem

Perhaps the earliest application of the discretization technique as it appeared in civil engineering practice was to the bar problem. The bar is a prismatic one-dimensional structure which can resist only axial loads, in tension or in compression, and cannot take any bending loads. Figure 1.1 shows its simplest configuration - a cantilever bar. This is a statically determinate problem, meaning that one can obtain a complete solution, displacements as well as strains and stresses from considerations of equilibrium alone. We shall now compute the analytical solution to the problem depicted in Fig.1.1 using very elementary engineering mathematics.

1.3.2.1 Analytical solution

The problem is categorized as a boundary value problem. We presume here that for this problem, the reader is able to understand that the governing differential equation describing the situation is EA u,xx = 0 (1.1)

where E is the Youngs modulus, A is the area of cross-section, u(x) is the function describing the variation of displacement of a point, and ,x denotes differentiation with respect to coordinate x. This is the equation of equilibrium describing the rate at which axial force varies along the length of

the bar. In general, governing differential equations belong to a category called partial differential equations but here we have a simpler form known as an ordinary differential equation. Equation (1.1) is further classified as a field equation or condition as one must find a solution to the variable u which must satisfy the equation over the entire field or domain, in this case, for x ranging from 0 to L. A solution is obviously

u = ax + b where a and b are as yet to determine a and b in conditions shown in Fig. conditions. Two boundary

(1.2)

undetermined constants. To complete the solution, i.e. a unique sense for the boundary (support) and loading 1.1, we must now introduce what are called the boundary conditions are needed here and we state them as

u = 0 at x = 0 EA u,x = P at x = L

(1.3a) (1.3b)

The first relates to the fact that the bar is fixed at the left end and the second denotes that a force P is applied at the free end. Taking these two conditions into account, we can show that the following description completely takes stock of the situation:

u(x) = Px/EA (x) = P/EA (x)= E(x)= P/A P(x) = A(x) = P

(1.4a) (1.4b) (1.4c) (1.4d)

where (x), (x) and P(x) are the strain, stress and axial force (stress resultants) along the length of the bar. These are the principal quantities that an engineer is interested in while performing stress analysis to check the integrity of any proposed structural design.

1.3.2.2 Approximate solutions

It is not always possible to determine exact analytical solutions to most engineering problems. We saw in the example above that an analytical solution provides a unique mathematical expression describing in complete detail the value of the field variable at every location in the body. From this expression, other expressions can be derived which describe further quantities of practical interest, e.g. strains, stresses, stress resultants, etc. In most problems where analytical solutions appear infeasible, it is necessary to resort to some approximate or numerical methods to obtain these values of practical interest. One obvious method at this stage is to start with the equations we have already derived, namely Equations (1.1) and (1.3). This can be achieved using a technique known as the finite difference method. We shall briefly describe this next.

Fig. 1.2 Grid for finite difference scheme

Other approximate solutions are based on what are called functional descriptions or variational descriptions of the problem. These describe the problem at a level which is more fundamental in the laws of physics than the more specialized descriptions in terms of governing differential equations as seen in Equations (1.1) and (1.3) above. We deal directly with energy or work quantities and discretization or approximation is applied at that level. If we can appreciate that the governing differential equations were derived analytically from such energy or work principles using what is called a variational or virtual work approach, then it would become easier to understand that an approximation applied directly to the energy or work quantities and a variation carried out subsequently on these terms will preserve both the physics as well as the approximating or numerical solution in one single process. This will be the basis for the finite element method as we shall see in subsequent subsections.

1.3.2.3 Finite difference approximations

For a problem where the differential equations are already known, the technique known as the finite difference method can be used to obtain an approximate or numerical solution. In fact, before the finite element method was established, it was this method which was used to obtain solutions for very complicated problems in structural mechanics, fluid dynamics and heat transfer. Figure 1.2 shows a very simple uniformly spaced mesh or grid of nodal points that can be used to discretize the problem described by Equations (1.1) and (1.3). The larger the number of nodal points used (i.e. the smaller the grid spacing h), the greater will be the accuracy involved. The field variable is now taken to be known (or very strictly, unknown, at this stage of the computation) at these nodal points. Thus, u1,..., ui,..., un, are the unknown variables or degrees of freedom. The next step is to rewrite the governing differential equations and boundary conditions in finite difference form. We see that the finite difference forms for u,x and u,xx are required. It is easy to show (again, details are omitted, as these can be found in most books on numerical analysis or fem) that at a grid point i. u,x = (ui+1 ui)/h u,xx = (ui+1 2ui + ui-1)/h
2

(1.5a) (1.5b)

We now attempt a solution in which four points are used in the finite difference grid. The governing differential equations and boundary conditions are replaced by the following discrete set of linear, algebraic equations: u1 = 0 u3 2u2 + u1 = 0 u4 - 2u3 + u2 = 0 = Ph/EA u4 - u3 (1.6a) (1.6b) (1.6c) (1.6d)

where h=L/3 and Equations (1.6a) and (1.6d) represent the boundary conditions at x = 0 and L respectively. The reader can easily work out that the solution obtained from this set of simultaneous algebraic equations is, u2 = PL/3EA, u3 = 2PL/3EA and u4 = PL/EA. Comparison with Equation (1.4a) will confirm that we have obtained the exact answers at the grid points. Other quantities of interest like strain, stress, etc. can be computed from the grid point values using the finite difference forms of these quantities. The above illustrates the solution to a very simple one-dimensional problem for which the finite difference procedure yielded the exact answer. Generalizations to two and three dimensions for more complex continuous problems can be made, especially if the meshes are of regular form and the boundary lines coincide with lines in such regular grids. For more complex shapes, the finite difference approach becomes difficult to use. It is in such applications that the finite element method proves to be more versatile than the finite difference scheme.

1.3.2.4

Variational formulation

We began this chapter with an analytical solution to the differential equations governing the problem. We did not ask how these equations originated. There are two main ways in which it can be done. The first, and the one that is most often introduced at the earliest stage of the study of structural mechanics, is to formulate the equations of equilibrium in terms of force or stress quantities. Then, depending on considerations of static determinacy, make the problem determinate by introducing stress-strain and strain-displacement relations until the final, solvable, set of governing equations is obtained. The second method is one that originates from a more fundamental statement of the principles involved. It recognizes that one of the most basic and elegant principles known to man is the law of least action, or minimum potential energy, or virtual work. It is a statement that Nature always chooses the path of minimum economy. Thus, the law of least action or minimum total potential is as axiomatic as the basic laws of equilibrium are. One can start with one axiom and derive the other or vice-versa for all of mechanics or most of classical physics. In dynamics, the equivalent principle is known as Hamilton's principle. In structural mechanics, we start by measuring the total energy stored in a structural system in the form of elastic strain energy and potential due to external loads. We then derive the position of equilibrium as that state where this energy is an extremum (usually a minimum if this is a stable state of equilibrium). This statement in terms of minimum energy is strictly true only for what are called conservative systems. For non-conservative systems, a more

general statement known as the principle of virtual work would apply. From these brief general philosophical reflections we shall now proceed to the special case at hand to see how the energy principle applies. The total energy is stated in the form of a functional. In fact, most of the problems dealt with in structural mechanics can be stated in a variational form as the search for a minimal or stationary point of a functional. The functional, , is an integral function of functions which are defined over the domain of the structure. Thus, for our simple bar problem, we must define terms such as the strain energy, U, and the potential due to the load, V, in terms of the field variable u(x). This can be written as

U =

L 1/ 2 EA u 2 ,x dx 0

(1.7a) (1.7b)

V = P ux=L Then,

( (x )) =

0 1/ 2 EA

u 2 ,x dx - P ux = L

(1.8)

Using standard variational procedures from the calculus of variations, we can show that the extremum or stationary value of the functional is obtained as = 0 (1.9)

If this variation is carried out, after integrating by parts and by regrouping terms, we will obtain the same equation of equilibrium and boundary conditions that appeared in Equations (1.1) and (1.3) above. Our interest is now not to show what can be done with these differential equations but to demonstrate that the approximation or discretization operation can be implemented at the stage where the total potential or functional is defined.

1.3.2.5

Functional approximation

There are several ways in which an approximation can be applied at the level of energy, virtual work or weighted residual statements, e.g. the Rayleigh-Ritz (R-R) procedure, the Galerkin procedure, the least-squares procedure, etc. We first choose a set of independent functions which satisfy the geometric or essential (also called kinematic) boundary conditions of the problem. These functions are called admissible functions. In this case, the condition specified in Equation (1.3a) is the geometric boundary condition. The other condition, Equation (1.3b) is called a nonessential or natural or force boundary condition. The admissible functions need not satisfy the natural boundary conditions (in fact, it can be shown that if they are made to do so, there can even be a deterioration in performance). The approximate admissible configuration for the field variable, say u (x ) , is obtained by a linear combination of these functions using unknown coefficients or parameters, also known as degrees of freedom or generalized coordinates. These unknown constants are then determined to satisfy the extremum or stationary statement. Methods like the Rayleigh-Ritz procedure

(the one going to be used now), the Galerkin method, collocation methods, least square methods, etc. all proceed with approximating fields chosen in this way. For our problem, let us choose only one simple linear function, using a polynomial form for the purpose. The advantages of using polynomial functions will become clear later in the book. Thus, only one constant is required, i.e. u (x ) = a1x (1.10)

Note that this configuration is an admissible one satisfying the geometric condition u (0 ) = 0 . By substituting in Equation (1.8) and carrying out the variation prescribed in Equation (1.9), we are trying to determine the value of a1 which provides the best functional approximation to the variational problem. It is left to the reader to show that the approximate solution obtained is u (x ) = Px/EA (1.11a) (1.11b) (1.11c) (1.11d)

(x ) = P/EA (x ) = E(x) = P/A

P (x ) = A(x) = P

Thus, for this problem, the approximate solution coincides with the exact solution. Obviously, we cannot draw any conclusions here as to the errors involved in approximation by the R-R procedure. In more complicated problems, approximate solutions will be built up from several constants and admissible functions and convergence to the exact solution will depend on the number of terms used and the type of functions used. The R-R procedure, when applied in a piecewise manner, over the element domains that constitute the entire structure, becomes the finite element method. Thus, an understanding of how the R-R method works will be crucial to our understanding of the finite element method.

1.3.2.6 Finite element approximation

We shall now describe an approximate solution by the finite element method. We shall use a form known as the displacement-type formulation (also called the stiffness formulation). The simplest element known for our purpose is a two noded line element such as shown in Fig. 1.3. This element is also known as the bar, truss or rod element. It has two nodes at which nodal displacements u1 and u2 are prescribed. These are also called the nodal degrees of freedom. Also indicated in Fig. 1.3 are the forces P1 and P2 acting at these nodes. Thus, any one-dimensional placed from end to end. The element is that stored as nodal now see how this piecewise Ritz structure can be replaced by contiguous elements only information communicated from element to degrees of freedom and the nodal forces. We can approximation is performed.

Fig. 1.3 A two-node bar element We first derive an expression for the approximate representation of the field variable, u (x ) , within the region of the element by relating it to the nodal degrees of freedom. It can be shown that,

u (x ) = u1N 1 + u2 N 2

(1.12)

where N1 =(1-)/2 and N2 =(1+)/2,where =x/l. N1 and N2 are called the shape functions. Note that this form is exactly the same as the linear function used in our R-R approximation earlier, Equation (1.10), with generalized coordinates ai being now replaced by nodal quantities ui. It is possible therefore to compute the energy stored in a beam element as

{u1 u2 }

EA 2l

1 - 1 u1 - 1 1 u2

(1.13)

and the potential due to the applied loads as

{u1 u2 } P1 P
2
Thus, if a variation is taken over u1 and equilibrium as EA 1 - 1 u1 = 1 u2 2l - 1

(1.14)

u2, we can write the equation of

P1 P2

(1.15)

This matrix representation is typical of the way finite element equations are assembled and manipulated automatically on the digital computer. We shall now attempt to solve our bar problem using a finite element model comprising just one element. Node 1 is therefore placed at the fixed end and node 2 coincides with the free end where the load P2 = P is applied. The fixity condition at node 1 is simply u1 = 0 and P1 indicates the reaction at this end. It is very simple to show from these that u2 = PL/EA. We can see that the exact solution has been obtained. What we have seen above is one of the simplest demonstrations of the finite element method that is possible. Generalizations of this approach to two

and three dimensional problems permit complex plate, shell and three dimensional elasticity problems to be routinely handled by ready made software packages.

1.4 Concluding remarks

In this chapter, we have briefly investigated how classical analytical techniques, which are usually very limited in scope, can give way to computational methods of obtaining acceptable solutions. We have seen that the finite element method is one such approximate procedure. By using a very large number of elements, it is possible to obtain solutions of greater accuracy. Over the years, finite element theorists and analysts have produced a very large body of work that shows how such solutions are managed for a great variety of problems. This is in the first-order tradition of the method. The remaining part of this book hopes to address questions like: How good is the approximate solution? What kinds of situation affect these approximations in a trivial and in not-so-trivial ways? Since solutions depend very much on the number of elements used and the type of shape function patterns used within each element, we must ask, what patterns are the best to assume? All these questions are linked to one very fundamental question: In what manner does the finite element method actually compute the solution? Does it make its first approximation of the displacement field directly or on the stress field? The exercises and exposition in this book are aimed at throwing more light on these questions.

10

Chapter 2 Paradigms and some approximate solutions

2.1 Introduction
We saw in the preceding chapter that a continuum problem in structural or solid mechanics can either be described by a set of partial differential equations and boundary conditions or as a functional based on the energy principle whose extremum describes the equilibrium state of the problem. Now, a continuum has infinitely many material points and therefore has infinitely many degrees of freedom. Thus, a solution is complete only if analytical functions can be found for the displacement and stress fields, which describe these states, exactly everywhere in the domain of the problem. It is not difficult to imagine that such solutions can be found only for a few problems. We also saw earlier that the Rayleigh-Ritz (RR) and finite element (fem) approaches offer ways in which approximate solutions can be achieved without the need to solve the differential equations or boundary conditions. This is managed by performing the discretization operation directly on the functional. Thus, a real problem with an infinitely large number of degrees of freedom is replaced with a computational model having a finite number of degrees of freedom. In the RR procedure, the solution is approximated by using a finite number of admissible functions i and a finite number of unknown constants ai so that the approximate displacement field is represented by a linear combination of these functions using the unknown constants. In the fem process, this is done in a piecewise manner - over each sub-region (element) of the structure, the displacement field is approximated by using shape functions Ni within each subregion and nodal degrees of freedom ui at nodes strategically located so that they connect the elements together without generating gaps or overlaps. The functional now becomes a function of the degrees of freedom (ai or ui as the case may be). The equilibrium configuration is obtained by applying the criterion that must be stationary with respect to the degrees of freedom. It is assumed that this solution process of seeking the stationary or extremum point of the discretized functional will determine the unknown constants such that these will combine together with the admissible or shape functions to represent some aspect of the problem to some "best" advantage. Which aspect this actually is has been a matter of some intellectual speculation. Three competing paradigms present themselves. It is possible to believe that by "best" we mean that the functions tend to satisfy the differential equations of equilibrium and the stress boundary conditions more and more closely as more terms are added to the RR series or more elements are added to the structural mesh. The second school of thought believes that it is the displacement field, which is approximated to greater accuracy with improved idealization. The "aliasing" paradigm, which will be critically discussed later, belongs to this school. It follows from this that stresses, which are computed as derivatives of the approximate displacement fields, will be less accurate. In this book however, we will seek to establish the currency of a third paradigm - that the RR or fem process actually seeks to find to best advantage, the state of stress or strain in the structure. In this

11

scenario, the displacement fields are computed from these "best-fit" stresses as a consequence. Before we enter into a detailed examination of the merits or faults of each of these paradigms, we shall briefly introduce a short statement on what is meant by the use of the term "paradigm" in the present context. We shall follow this by examining a series of simple approximations to the cantilever bar problem but with more and more complex loading schemes to see how the overall picture emerges.

2.2 What is a "paradigm"?

Before we proceed further it may be worthwhile to state what we mean by a paradigm here. This is a word that is uncommon to the vocabulary of a trained engineer. The dictionary meaning of paradigm is pattern or model or example. This does not convey much in our context. Here, we use the word in the greatly enlarged sense in which the philosopher T. S. Kuhn introduced it in his classic study on scientific progress [2.1]. In this sense, a paradigm is a "framework of suppositions as to what constitutes problems, theories and solutions. It can be a collection of metaphysical assumptions, heuristic models, commitments, values, hunches, which are all shared by a scientific community and which provide the conceptual framework within which they can recognize problems and solve them [2.2]. The aliasing and best-fit paradigms can be thought of as two competing scenarios, which attempt to explain how the finite element method computes strains and stresses. Our task will therefore be to establish which paradigm has greater explanatory power and range of application.

Fig. 2.1 Bar under uniformly distributed axial load

12

2.3 Bar under uniformly distributed axial load

Consider a cantilever bar subjected to a uniformly distributed axial load of intensity q 0 per unit length (Fig. 2.1). Starting with the differential equation of equilibrium, it is easy to show that the analytical solution to the problem is

u (x ) = (q 0 EA ) Lx x 2 2

(2.1a)

(x ) = (q 0 A ) (L x )

(2.1b)

Consider a one-term RR solution based on ur = a1x , where the subscript r denotes the use of the RR approach. It is left to the reader to show that the approximate solution obtained is ur (x ) = (q 0 EA ) (Lx 2 ) (2.2a)

r (x ) = q 0 A

) (L 2 )

(2.2b)

We now compare this with an fem solution based on a two-noded linear element. The solution obtained will be

uf (x ) = (q 0 EA ) (Lx 2 )

(2.3a)

f (x ) = q 0 A

) (L 2 )

(2.3b)

We see that the RR and fem solutions are identical. This is to be expected because the fem solution is effectively an RR solution. We can note also the curious coincidence where all three solutions predict the same displacement at the tip. However, from Fig. 2.1 we can see that the ur and uf are approximate solutions to u. It is also clear from Fig. 2.1 that r = f = at the mid-point of the beam. It is possible to speculate that r and f bear some unique relationship to the true variation . Consider next what will happen if two equal length linear (i.e., twonoded) bar elements are used to model the bar. It is left to the reader to work out that the solution described in Fig. 2.2 will be obtained. First, we must note that the distributed axial load is consistently lumped at the nodes. Thus the physical load system that the fem equations are solving is not that described in Fig. 2.1 or Fig. 2.2 as . Instead, we must think of a substitute stairstep distribution f, produced by the consistent load lumping process (see Fig. 2.2) which is sensed by the fem stiffness matrix. Now, a solution of the set of algebraic equations will result in f and uf as the fem solution. We see once again that the nodal predictions are exact. This is only a coincidence for this particular type of problem and nothing more can be read into this fact. More striking is the observation that the stresses computed by

13

the finite element stairstep fashion.

system

now

approximates

the

original

true

stress

in

It also seems reasonable to conclude that within each element, the true state of stress is captured by the finite element stress in a "best-fit" sense. In other words, we can generalize from Figs. 2.1 and 2.2, that the finite element stress magnitudes are being computed according to some precise rule. Also, there is the promise that by carefully understanding what this rule is, it will be possible to derive some unequivocal guidelines as to where the stresses are accurate. In this example, where an element capable of yielding constant stresses is used to model a problem where the true stresses vary linearly, the centroid of the element yields exact predictions. As we take up further examples later, this will become more firmly established. A cursory comparison of Figs. 2.1 and 2.2 also indicates that in a general sense, the approximate displacements are more accurate than the approximate stresses. It seems compelling now to argue that this is so because the approximate displacements emerge as "discretized" integrals of the stresses or strains and for that reason, are more accurate than the stresses.

2.4 Bar under linearly varying distributed axial load

We now take up a slightly more difficult problem. The cantilever bar has a load distributed in a linearly varying fashion (Fig. 2.3). It is easy to establish that the exact stress distribution in this case will be quadratic in nature. We can write it as

(x ) = q 0 L2 8A 4 3 2 - 1 - 3 2 3

(2.4)

Some interesting features about this equation can be noted. A dimensionless coordinate, = 2x/L-1 has been chosen so that it will also serve as a natural coordinate system taking on values -1 and 1 at the ends of the single three-node bar element shown as modeling the entire bar in Fig.2.3. We have also very 2 curiously expanded the quadratic variation using the terms 1, , (1-3 ). These can be identified with the Legendre polynomials and its relevance to the treatment here will become more obvious as we proceed further.

Fig. 2.2 Two element model of bar under uniformly distributed axial load.

14

Fig. 2.3 Bar under linearly varying axial load

We shall postpone the first obvious approximation, that of using a oneterm series ur = a1x till later. For now, we shall consider a two-term series

ur = a1x + a2 x 2 . This is chosen so that the essential boundary conditions at x=0 is satisfied. No attempt is made to satisfy the force boundary condition at x=L. The reader can verify by carrying out the necessary algebra associated with the RR process that the solution obtained will yield an approximate stress pattern given by

r (x ) = q 0 L2 8A (4 3 2 )

(2.5)

This is plotted on Fig. 2.3 as the dashed line. A comparison of Equations (2.4) and (2.5) reveals an interesting fact - only the first two Legendre polynomial terms are retained. Taking into account the fact that the Legendre polynomials are orthogonal, what this means is that in this problem, we have obtained r in a manner that seems to satisfy the following integral condition: (2.6) That is, the RR procedure has determined a r that is a "best-fit" of the true state of stress in the sense described by the orthogonality condition in Equation (2.6). This is a result anticipated from our emerging results to the various exercises we have conducted so far. We have not been able to derive it from any general principle that this must be so for stronger reasons than shown here up till now.
1 -1 r ( r ) d = 0

15

Let us now proceed to an fem solution. It is logical to start here with the three-node element that uses the quadratic shape functions, N = 1 2 ,

N2 = 1

) and

N 3 = + 1 2. We first compute the consistent loads to be placed

at the nodes due to the distributed load using Pi =

N iq dx . This results in the

2 following scheme of loads at the three nodes identified in Fig. 2.3: P1 = q 0 L 6 ,

P2 = q 0 L2 3 and P3 = 0 . The resulting load configuration can be represented in the form of a stress system shown as f, represented by the dashed-dotted lines in Fig. 2.3. Thus, any fem discretization automatically replaces a smoothly varying stress system by a step-wise system as shown by f in Figs. 2.2 and 2.3. It is this step-wise system that the finite element solution f responds to. If the finite element computations are actually performed using the stiffness matrix for the three-node element and the consistent load vector, it turns out, as the reader can assure himself, that the computed fem stress, pattern will be
(2.7) (4 3 2 ) This turns out to be exactly the same as the r computed by the RR process in Equation (2.5). At first sight, this does not seem to be entirely unexpected. Both the RR process and the fem process here have started out with quadratic admissible functions for the displacement fields. This implies that both have the capability to represent linear stress fields exactly or more complicated stress fields by an approximate linear field that is in some sense a best approximation. On second thought however, there is some more subtlety to be taken care of. In the RR process, the computed r was responding to a
2 f x = q 0 L 8A

( )

quadratically varying system (see Fig. 2.3). We could easily establish through Equation (2.6) that r responded to in a "best-fit" manner. However, in the fem process, the load system that is being used is the f system, which varies in the stairstep fashion. The question confronting us now is, in what manner did f respond to f - is it also consistent with the "best-fit" paradigm? Let us now assume an unknown field = c0 + c1 which is a "best-fit" of

the stairstep field given by f = 2q 0 L2 3 in 0 <x <L/2 and f = 0 in L/2 <x <L. We shall determine the constants c0 and c1 so that the "best-fit" variation shown below is satisfied:

2q 0 L2 3 d +

0 0 d = 0

(2.8)

The reader can work out that this leads to

(x ) = q 0 L2 8A

) (4 3 2 ) ,
16

(2.9)

which is identical to the result obtained in Equation (2.7) by carrying out the finite element process. In other words, the fem process follows exactly the "best-fit" description of computing stress fields. Another important lesson to be learnt from this exercise is that the consistent lumping process preserves the "best-fit" nature of the stress representation and subsequent prediction. Thus, f is a "best-fit" of both and f! Later, we shall see how the best-fit nature gets disturbed if the loads are not consistently derived. It again seems reasonable to argue that the nodal displacements computed directly from the stiffness equations from which the stress field f has been processed can actually be thought of as being "integrated" from the "best-fit" stress approximation. It seems possible now to be optimistic about our ability to ask and answer questions like: What kind of approximating fields are best to start with? and, How good are the computed results? Before we pass on, it is useful for future reference to note that the approximate solutions r or f intersect the exact solution at two points. A comparison of Equation (2.4) with Equations (2.5) and (2.7) indicates that these are the points where the quadratic Legendre polynomial,

(1 3 )
2

vanishes, i.e.,

at = 1 3 . Such points are well known in the literature of the finite element method as optimal stress points or Barlow points, named after the person who first detected such points. Our presentation shows clearly why such points exist, and why in this problem, where a quadratic stress state is sought to be approximated by a linear stress state, these points are at = 1 3 . Curiously, these are also the points used in what is called the two-point rule for GaussLegendre numerical integration. Very often, the two issues are confused. Our interpretation so far shows that the fact that such points are sampling points of the Gauss-Legendre rule has nothing to do with the fact that similar points arise as points of optimal stress recovery. The link between them is that in both, the Legendre polynomials play a decisive role - the zeros of the Legendre polynomials define the optimal points for numerical integration and these points also help determine where the "best-fit" approximation coincides with a polynomial field which is exactly one order higher. We shall now return to the linear Ritz admissible function, ur = a1x , to see if it operates in the best-fit sense. This would be identical to using a single two-node bar element to perform the same function. Such a field is capable of representing only a constant stress. This must now approximate the quadratically varying stress field (x) given by Equation (2.4). This gives us an opportunity to observe what happens to the optimal stress point, whether one exists, and whether it can be easily identified to coincide with a Gauss integration point. Again, the algebra is very simple and is omitted here. One can show that the one-term approximate solution would lead to the following computed stress:

r (x ) = q 0 L2 8A

) (4 3 )

(2.10)

17

What becomes obvious by comparing this with the true stress (x) in Equation (2.4) and the computed stress from the two-term solution, r (x ) in Equation (2.5) is that the one-term solution corresponds to only the constant part of the Legendre polynomial expansion! Thus, given the orthogonal nature of the Legendre polynomials, we can conclude that we have obtained the "best-fit" state of stress even here. Also, it is clear that the optimal stress point is not easy to identify to coincide with any of the points corresponding to the various Gauss integration rules. The optimal point here is given by = 1 4 3 .

2.5 The aliasing paradigm

We have made out a very strong case for the best-fit paradigm. Let us now examine the merits, if any, of the competing paradigms. The argument that fem procedures look to satisfy the differential equations and boundary conditions does not seem compelling enough to warrant further discussion. However, the belief that finite elements seek to determine nodal displacements accurately was the basis for Barlow's original derivation of optimal points [2.3] and is also the basis for what is called the "aliasing" paradigm [2.4]. The term aliasing is borrowed from sample data theory where it is used to describe the misinterpretation of a time signal by a sampling device. An original sine wave is represented in the output of a sampling device by an altered sine wave of lower frequency - this is called the alias of the true signal. This concept can be extended to finite element discretization - the sample data points are now the values of the displacements at the nodes and the alias is the function, which interpolates the displacements within the element from the nodal displacements. We can recognize at once that Barlow [2.3] developed his theory of optimal points using an identical idea - the term substitute function is used instead of alias. Let us now use the aliasing concept to derive the location of the optimal points, as Barlow did earlier [2.3], or as MacNeal did more recently [2.4]. We assume here that the finite element method seeks discretized displacement fields, which are substitutes or aliases of the true displacement fields by sensing the nodal displacements directly. We can compare this with the "best-fit" interpretation where the fem is seen to seek discretized strain/stress fields, which are the substitutes/aliases of the true strain/stress fields in a "best fit" or "best approximation" sense. It is instructive now to see how the alternative paradigm, the "displacement aliasing" approach leads to subtle differences in interpreting the relationship between the Barlow points and the Gauss points. Table 2.1 Barlow and Gauss points for one-dimensional case p 1 2 3 Nodes at

2 2

Gauss points

Barlow points best-fit aliasing

1 0, 1

3
2

0 1/3

0 1/3 0, 3 5

0 1/3
0, 5/3

3 1/3, 1 4 3 2 0, 3 5 4 1, .., indicate polynomial orders

18

2.5.1 A one dimensional problem

We again take up the simplest problem, a bar under axial loading. We shall assume that the bar is replaced by a single element of varying polynomial order for its basis function (i.e. having varying no. of equally spaced nodes). Thus, from Table 2.1, we see that p=1,2,3 correspond to basis functions of linear, quadratic and cubic order, implying that the corresponding elements have 2,3,4 nodes respectively. These elements are therefore capable of representing a constant, linear and quadratic state of strain/stress, where strain is taken to be the first derivative of the displacement field. We shall adopt the following notation: The true displacement, strain and stress field will be designated by u, and . The discretized displacement, strain and stress field will be designated by u , and . The aliased displacement, strain and stress field will be designated by u ui.
a

, a and a. Nodal displacements will be represented by

We shall examine three scenarios. In the simplest, Scenario a, the true displacement field u is exactly one polynomial order higher than what the finite element is capable of representing - we shall see that the Barlow points can be determined exactly in terms of the Gauss points only for this case. In Scenario b, we consider the case where the true field u is two orders higher than the discretized field u . The definition of an identifiable optimal point becomes difficult. In both Scenarios a and b, we assume that the rigidity of the bar is a constant, i.e., = D . In Scenario c, we take up a case where the rigidity can vary, i.e., = D() . We shall see that once again, it becomes difficult to identify the optimal points by any simple rule. We shall now take for granted that the best-fit rule operates according to the orthogonality condition expressed in Equation (2.6) and that it can be used interchangeably for stresses and strains. We shall designate the optimal points determined by the aliasing algorithm as a, the Barlow points (aliasing), and the points determined by the best-fit algorithm as b, the Barlow points (bestfit). Note that a are the points established by Barlow [2.3] and MacNeal [2.4], while b will correspond to the points given in References 2.6 and 2.7. Note that the natural coordinate system is being used here for convenience. Thus, for Scenarios a and b, this leads to
T ( ) dV = 0

(2.11)

whereas for Scenario c, it becomes

D ( ) ( ) dV = 0

(2.12)

Note that we can consider Equation (2.11) as a special case of the orthogonality condition in Equation (2.12) with the weight function D()=1. This case corresponds to one in which a straightforward application of Legendre polynomials can be made. This point was observed very early by Moan [2.5]. In this case, one can determine the points where = as those corresponding to points, which are the zeros of the Legendre polynomials. See Table 2.2 for a list of unnormalised Legendre polynomials. We shall show below that in Equation

19

(2.11), the points of minimum error are the sampling points of the Gauss Legendre integration rule only if is exactly one polynomial order lower than

. And in Equation (2.12), the optimal points no longer depend on the nature of the Legendre polynomials, making it difficult to identify the optimal points.
Scenario a
We shall consider fem solutions using a linear (two-node), a quadratic (threenode) and a cubic (four-node) element. The true displacement field is taken to be one order higher than the discretized field in each case. Linear element (p = 1)

u = quadratic = bo + b1 + b2 = linear = u, = b1 + 2b2 =

2
p =1 s=0

(2.13)

s Ps ( )

(2.14)

Note that we have written in terms of the Legendre polynomials for future convenience. Note also that we have simplified the algebra by assuming that strains can be written as derivatives in the natural co-ordinate system. It is now necessary to work out how the algebra differs for the aliasing and best-fit approaches.
a Aliasing: At i = 1 , ui = ui ; then points where a = are given by a=0. (The algebra is very elementary and is left to the reader to work out). Thus, the Barlow point (aliasing) is a = 0, for this case.

Best-fit: u =linear, is undetermined at first. Let = c0 , as the element is capable of representing only a constant strain. Equation (2.11) will now give = c0 = b1 . Thus, the optimal point is b = 0, the point where the Legendre
polynomial P1() = vanishes. Therefore, the Barlow point (best-fit) for this example is b = 0.

Table 2.2

The Legendre polynomials Pi Polynomial Pi 1

Order of polynomial i 0 1 2 3 4

(1-3 )
2

(3-5 )
3

(3-30 +35
2

20

Quadratic element (p = 2)

u = cubic = b0 + b1 + b2

+ b3

= quadratic = u,

(2.15)

= (b1 + b3 ) + 2b2 b3 1 3 2 =

s=0

sPs ( )

(2.16)

a Aliasing: At i = 0, 1, ui = ui ; ; then points where a = are given by a = 1 3 . (Again, the algebra is left to the reader to work out). Thus, the Barlow points

(aliasing) are a = 1

3 , for this case.

Let

= c0 + c1 as this element is capable of representing a linear strain. Equation (2.11) will now give = (b 1 + b 3 ) + 2b 2 .
Best-fit: u
=

quadratic.

Thus,

the

optimal

polynomial P2 ( ) = 1 3 2 this example are b = 1

points

are

b = 1

3,

the

points

where

the

Legendre

vanishes. Therefore, the Barlow points (best-fit) for

3.

Note that in these two examples, i.e. for the linear and quadratic elements, the Barlow points from both schemes coincide with the Gauss points (the points where the corresponding Legendre polynomials vanish). In our next example we will find that this will not be so. Cubic element (p = 3)

u = quartic = b0 + b1 + b2 2 + b3 3 + b4 4

(2.17)

= cubic = u,
= (b1 + b3 ) + (2b2 + 12b4 5 ) b3 1 3 2 4b4 5 3 5 3
p =3

)
(2.18)

s=0

sPs ( )
At
a i = 1 3 , 1, ui = ui ;

Aliasing:

then

points

where

a =

are

given

by

a = 0, 5 3 . Thus, the Barlow points (aliasing) are a = 0, 5 3 for this case.

Note that the points where the Legendre polynomial P3 ( ) = 3 5 3

vanishes are

0 = 0, 3 5 !
Best-fit: u = cubic. Let = c0 + c1 + c2 1 3 2 , as this element is capable of representing a quadratic strain. Equation (2.11) will now give

= (b1 + b3 ) + (2b2 + 12b4 5 ) b3 1 3 2 . Thus, the Barlow points (best-fit) for

21

this

example

P3 ( ) = 3 5

are

b = 0, 3 5 ,

the

points

where

the

Legendre

polynomial

vanishes.

Therefore, we have an example where the aliasing paradigm does not give the correct picture about the way the finite element process computes strains. However, the best-fit paradigm shows that as long as the discretized strain is one order lower than the true strain, the corresponding Gauss points are the optimal points. Table 2.1 summarizes the results obtained so far. The experience of this writer and many of his colleagues is that the best -fit model, is the one that corresponds to reality. If one were to actually solve a problem where the true strain varies cubically using a 4-noded element, which offers a discretized strain which is of quadratic order, the points of optimal strain actually coincide with the Gauss points.

Scenario b:
So far, we have examined simple scenarios where the true strain is exactly one polynomial order higher than the discretized strain with which it is replaced. If P p ( ) , denoting the Legendre polynomial of order p, describes the order by which the true strain exceeds the discretized strain, the simple rule is that he optimal points are obtained as P p ( b ) = 0 . These are therefore the set of p Gauss points at which the Legendre polynomial of order p vanishes. Consider now a case where the true strain is two orders higher than the discretized strain - e.g. a quadratic element (p = 2) modeling a region where the strain and stress field vary cubically. Thus, we have,

= (b1 + b3 ) + (2b2 + 12b4 5 ) b3 1 3 2 4b4 5 3 5 3 = c0 + c1

)
(2.19)

Equation (2.11) allows us to determine the coefficient ci in terms of bi; it turns out that

= (b1 + b3 ) + (2b2 + 12b4 5 ) ,

(2.20)

a representation made very easy by the fact that the Legendre polynomials area orthogonal and that therefore can be obtained from by simple inspection. It is not however a simple matter to determine whether the optimal points coincide with other well-known points like the Gauss points. In this example, we have to seek the zeros of

b3 1 3 2

) + 4b 5 (3 5 )
3 4

(2.21)

Since b3 and b4 are arbitrary, depending on the problem, it is not possible to seek universally valid points where this would vanish, unlike in the case of Scenario a earlier. Therefore, in such cases, it is not worthwhile to seek points of optimal accuracy. It is sufficient to acknowledge that the finite

22

element procedure yields strains obtain in the circumstances.

, which are the most reasonable one can

Scenario c
So far, we have confined attention to problems where is related to through a simple constant rigidity term. Consider an exercise where (the one-dimensional bar again) the rigidity varies because the cross-sectional area varies or because the modulus of elasticity varies or both, i.e. = D(). The orthogonality condition that governs this case is given by Equation (2.12). Thus, it may not be possible to determine universally valid Barlow points a priori if D() varies considerably.

2.6 The `best-fit' rule from a variational theorem

Our investigation here will be complete in all respects if the basis for the best-fit paradigm can be derived logically from a basic principle. In fact, some recent work [2.6,2.7] shows how by taking an enlarged view of the variational basis for the displacement type fem approach we will be actually led to the conclusion that strains or stresses are always sought in the best-fit manner. The `best-fit' manner in which finite elements compute strains can be shown to follow from an interpretation uses the Hu-Washizu theorem. To see how we progress from the continuum domain to the discretized domain, we will find it most convenient to develop the theory from the generalized Hu-Washizu theorem [2.8] rather than the minimum potential theorem. As we have seen earlier, these theorems are the most basic statements that can be made about the laws of nature. The minimum potential theorem corresponds to the conventional energy theorem. However, for applications to problems in structural and solid mechanics, Hu proposed a generalized theorem, which had somewhat more flexibility [2.8]. Its usefulness came to be recognized when one had to grapple with some of the problems raised by finite element modeling. One such puzzle is the rationale for the "best-fit" paradigm. Let the continuum linear elastic problem have an exact solution described by the displacement field u, strain field and stress field . (We project that the strain field is derived from the displacement field through the straindisplacement gradient operators of the theory of elasticity and that the stress field is derived from the strain field through the constitutive laws.) Let us now replace the continuum domain by a discretized domain and describe the computed state to be defined by the quantities u , and , where again, we take that the strain fields and stress fields are computed from the straindisplacement and constitutive relationships. It is clear that is an approximation of the true strain field . What the Hu-Washizu theorem does, following the interpretation given by Fraejis de Veubeke [2.9], is to introduce a dislocation potential to augment the usual total potential. This dislocation which can be potential is based on a third independent stress field considered to be the Lagrange multiplier removing the lack of compatibility appearing between the true strain field and the discretized strain field .

23

Note that, is now an approximation of . The three-field Hu-Washizu theorem can be stated as

= 0
where

(2.22)

T 2 +

( )

+ P dv

(2.23)

and P is the potential energy of the prescribed loads. In the simpler minimum total potential principle, which is the basis for the derivation of the displacement type finite element formulation in most textbooks, only one field (i.e., the displacement field u) is subject to variation. However, in this more general three field approach, all the three fields are subject to variation and leads to three sets of equations which can be grouped and classified as follows:

Variation on

Nature Equilibrium

Equation

+ terms from P = 0

(2.24a)

Orthogonality (Compatibility )

) dv = 0

(2.24b)

Orthogonality (Equilibrium )

dv = 0

(2.24c)

Equation (2.24a) shows that the variation on the displacement field u requires that the independent stress field must satisfy the equilibrium equations ( signifies the operators that describe the equilibrium condition. Equation (2.24c) is a variational condition to restore the equilibrium imbalance between

and . In the displacement type the orthogonality condition seen in with the orthogonality condition in tries to restore the compatibility and the discretized strain field can be stated as,

formulation, we choose = . This satisfies Equation (2.24c) identically. This leaves us Equation (2.24b). We can now argue that this imbalance between the exact strain field . In the displacement type formulation this

( ) dV

=0

(2.25)

Thus we see very clearly, that the strains computed by the finite element procedure are a variationally correct (in a sense, a least squares correct) `best approximation' of the true state of strain.

24

2.7 What does the finite element method do?

To understand the motivation for the aliasing paradigm, it will help to remember that it was widely believed that the finite element method sought approximations to the displacement fields and that the strains/stresses were computed by differentiating these fields. Thus, elements were believed to be "capable of representing the nodal displacements in the field to a good degree of accuracy." Each finite element samples the displacements at the nodes, and internally, within the element, the displacement field is interpolated using the basis functions. The strain fields are computed from these using a process that involves differentiation. It is argued further that as a result, displacements are more accurately computed than the strain and stress field. This follows from the generally accepted axiom that derivatives of functions are less accurate than the original functions. It is also argued that strains/stresses are usually most inaccurate at the nodes and that they are of greater accuracy near the element centers - this, it is thought, is a consequence of the mean value theorem for derivatives. However, we have demonstrated convincingly that in fact, the Ritz approximation process, and the displacement type fem which can be interpreted as a piecewise Ritz procedure, do exactly the opposite - it is the strain fields which are computed, almost independently as it were within each element. This can be derived in a formal way. Many attempts have been made to give expression to this idea, (e.g., Barlow [2.1] and Moan [2.7]), but it appears that the most intellectually satisfying proof can be arrived at by starting with a mixed principle known as the Hu-Washizu theorem [2.8]. This proof has been taken up in greater detail in Reference 2.3 and was reviewed only briefly in the preceding section. Having said that the Ritz type procedures determine strains, it follows that the displacement fields are then constructed from this in an integral sense. The stiffness equation actually reflecting this integration process and the continuity of fields across element boundaries and suppression of the field values at domain edges being reflected by the imposition of boundary conditions. It must therefore be argued that displacements are more accurate than strains because integrals of smooth functions are generally more accurate than the original data. We have thus turned the whole argument on its head.

2.8 Conclusions
In this chapter, we postulated a few models to explain how the displacements type fem works. We worked out a series of simple problems of increasing complexity to establish that our conjecture that strains and stresses appear in a "best-fit" sense could be verified (falsified, in the Popperian sense) by carefully designed numerical experiments. An important part of this exercise depended on our careful choice and use of various stress terms. Thus terms like and f was actual physical states that were sought to be modeled. The stress terms r and f were the quantities that emerged in what we can call the "first-order tradition" analysis in the language of Sir Karl Popper [2.10] - where the RR or fem operations are mechanically carried out using functional approximation and finite element stiffness equations respectively. We noticed certain features that seemed to relate these computed stresses to the true system they were modeling in a predictable or repeatable manner. We then proposed a mechanism to explain how

25

this could have taken place. Our bold conjecture, after examining these numerical experiments, was to propose that it is effectively seeking a best-fit state. To confirm that this conjecture is scientifically coherent and complete, we had to enter into a "second-order tradition" exercise. We assumed that this is indeed the mechanism that is operating behind the scenes and derived predicted quantities r that will result from the best-fit paradigm when this was applied to the true state of stress. These predicted quantities turned out to be exactly the same as the quantities computed by the RR and fem procedures. In this manner, we could satisfy ourselves that the "best-fit" paradigm had successfully survived a falsification test. Another important step we took was to prove that the "best-fit" paradigm was neither gratuitous nor fortuitous. In fact, we could also establish that this could be derived from more basic principles - in this regard, the generalized theorem of Hu was found valuable to determine that the best-fit paradigm had a rational basis. In subsequent chapters, we shall use this foundation to examine other features of the finite element method. One important conclusion we can derive from the best-fit paradigm is that, an interpolation field for the stresses (or stress resultants as the case may be) which is of higher order than the strain fields . On which it must `do work' in the energy or virtual work principle is actually self-defeating because the higher order terms cannot be `sensed'. This is precisely the basis for de Veubeke's famous limitation principle [2.9], that it is useless to look for a better solution by injecting additional degrees of freedom in the stresses. We can see that one cannot get stresses, which are of higher order than are reflected in the strain expressions.

2.9 References
2.1. T. S. Kuhn, The Structure of Scientific Revolution, University of Chicago Press, 1962. 2.2. S. Dasgupta, Understanding design: Artificial explanatory paradigm, SADHANA, 19, 5-21, 1994. intelligence as an

2.3. J. Barlow, Optimal stress locations in finite element models, Int. J. Num. Meth. Engng. 10, 243-251 (1976). 2.4. R. H. MacNeal, Finite Dekker, NY, 1994.

Elements:

Their

Design

and

Performance,

Marcel

2.5

T. Moan, On the local distribution of errors by finite element approximations, Theory and Practice in Finite Element Structural Analysis. Proceedings of the 1973 Tokyo Seminar on Finite Element Analysis, Tokyo, 43-60, 1973.

2.6. G. Prathap, The Finite Element Method in Structural Mechanics, Kluwer Academic Press, Dordrecht, 1993.

26

2.7. G. Prathap, A variational basis for the Barlow points, Comput. Struct. 49, 381-383, 1993. 2.8. H. C. Hu, On some variational principles in the theory of elasticity and the theory of plasticity, Scientia Sinica, 4, 33-54 (1955). 2.9. B. F. de Veubeke, Displacement and equilibrium models in the finite element method, in Stress Analysis, Ellis Horwood, England, 1980. 2.10. B. Magee, Popper, Fontana Press, London, 1988.

27

Chapter 3 Completeness and continuity: How to choose shape functions?

3.1 Introduction
We can see from our study so far that the quality of approximation we can achieve by the RR or fem approach depends on the admissible assumed trial, field or shape functions that we use. These functions can be chosen in many ways. An obvious choice, and the one most universally preferred is the use of simple polynomials. It is possible to use other functions like trigonometric functions but we can see that the least squares or best-fit nature of stress prediction that the finite element process seeks in an instinctive way motivates us to prefer the use of polynomials. When this is done, interpretations using Legendre polynomials become very convenient. We shall see this happy state of affairs appearing frequently in our study here. An important question that will immediately come to mind is - How does the accuracy or efficiency of the approximate solution depends on our choice of shape functions. Accuracy is represented by a quality called convergence. By convergence, we mean that as we add more terms to the RR series, or as we add more nodes and elements into the mesh that replaces the original structure, the sequence of trial solutions must approach the exact solution. We want quantities such as displacements, stresses and strain energies to be exactly recovered, surely, and if possible, quickly. This leads to the questions: What kind of assumed pattern is best for our trial functions? What are the minimum qualities, or essences or requirements that the finite element shape functions must show or meet so that convergence is assured. Two have been accepted for a long time: continuity (or compatibility) and completeness.

3.2 Continuity
This is very easy to visualize and therefore very easy to understand. A structure is sub-divided into sub-regions or elements. The overall deformation of the structure is built-up from the values of the displacements at the nodes that form the net or grid and the shape functions within elements. Within each element, compatibility of deformation is assured by the fact that the choice of simple polynomial functions for interpolation allows continuous representation of the displacement field. However, this does not ensure that the displacements are compatible between element edges. So special care must be taken otherwise, in the process of representation, gaps or overlaps will develop. The specification of continuity also depends on how strains are defined in terms of derivatives of the displacement fields. We know that a physical problem can be described by the stationary condition =0, where =() is the functional. If contains derivatives of the field variable to the order m, then we obviously require that within each element, , which is the approximate field chosen as trial function, must contain a complete polynomial of degree m. So that is continuous within elements and the completeness requirements (see below) are also met. However, the more important requirement now is that continuity of field variable must be maintained across element boundaries -

28

this requires that the trial function and its derivatives through order m-1 must be continuous across element edges. In most natural formulations in solid mechanics, strains are defined by first derivatives of the displacement fields. In such cases, a simple continuity of the displacement fields across element 0 edges suffices - this is called C continuity. Compatibility between adjacent 0 elements for problems which require only C -continuity can be easily assured if the displacements along the side of an element depend only on the displacements specified at all that nodes placed on that edge. There are problems, as in the classical Kirchhoff-Love theories of plates and shells, where strains are based on second derivatives of displacement fields. In this case, continuity of first derivatives of displacement fields 1 across element edges is demanded; this is known as C -continuity. Elements, which satisfy the continuity conditions, are called conforming or compatible elements. We shall however find a large class of problems (Timoshenko beams, Mindlin plates and shells, plane stress/strain flexure, incompressible elasticity) where this simplistic view of continuity does not assure reasonable (i.e. practical) rates of convergence. It has been noticed that formulations that take liberties, e.g. the non-conforming or incompatible approaches, significantly improve convergence. This phenomenon will be taken up for close examination in a later chapter.

3.3 Completeness
We have understood so far that in the finite element method, or for that matter, in any approximate method, we are trying to replace an unknown function (x), which is the exact solution to a boundary value problem over a domain enclosed by a boundary by an approximate function (x ) which is constituted from a set of trial, shape or basis functions. We desire a trial function set that will ensure that the approximation approaches the exact solution as the number of trial functions is increased. It can be argued that the convergence of the trial function set to the exact solution will take place if (x ) will be sufficient to represent any well behaved function such as (x) as closely as possible as the number of functions used becomes indefinitely large. This is called the completeness requirement. In the finite element context, where the total domain is sub-divided into smaller sub-regions, completeness must be assured for the shape functions used within each domain. The continuity requirements then provide the compatibility of the functions across element edges. We have also seen that what we seek is a best-fit arrangement in some sense between (x ) and (x). From Chapter 2, we also have the insight that this best-fit can be gainfully interpreted as taking place between strain or stress quantities. This has important implications in further narrowing the focus of the completeness requirement for finite element applications in particular. By bringing in the new perspective of a best-fit strain or stress paradigm, we are able to look at the completeness requirements entirely from

29

what we desire the strains or stresses to be like. Of paramount importance now is the idea that the approximate strains derived from the set of trial functions chosen must be capable in the limit of approximation (i.e. as the number of terms becomes very large) to describe the true strain or stress fields exactly. This becomes very clear from the orthogonality condition derived in Equation 2.25. From the foregoing discussion, we are now in a position to make a more meaningful statement about what we mean by completeness. We would like displacement functions to be so chosen that no straining within an element takes place when nodal displacements equivalent to a rigid body motion of the whole element are applied. This is called the strain free rigid body motion condition. In addition, it will be necessary that each element must be able to reproduce a state of constant strain; i.e. if nodal displacements applied to the element are compatible with a constant strain state, this must be reflected in the strains computed within the element. There are simple rules that allow these conditions to be met and these are called the completeness requirements. If polynomial trial functions are used, then a simple assurance that the polynomial functions 0 contain the constant and linear terms, etc. (e.g. 1, x in a one-dimensional C 0 problem;) 1, x, y in a two-dimensional C problem so that each element is certain to be able to recover a constant state of strain, will meet this requirement.

3.3.1 Some properties of shape functions for C0 elements

The displacement field is approximated by using shape functions Ni within each element and these are linked to nodal degrees of freedom ui. It is assumed that these shape functions allow elements to be connected together without generating gaps or overlaps, i.e. the shape functions satisfy the continuity requirement described above. Most finite element formulations in practical use in general 0 purpose application require only what is called C continuity - i.e. a field quantity must have interelement continuity but does not require the continuity of all of its derivatives across element boundaries.

C shape functions have some interesting properties that derive from the completeness requirements imposed by rigid body motion considerations. Consider a field (usually a displacement) which is interpolated for an n-noded element according to

i=1

N ii

(3.1) natural co-ordinates (e.g. , argue that one property the the distribution of within at node i has unit value and

where the Ni are assumed to be interpolated using for a 2-dimensional problem). It is simple to shape functions must have is that Ni must define the element domain when the degree of freedom i all other nodal s are zero. Thus

i) Ni = 1 when x = xi (or = i) and Ni = 0 when x = xj (or = j) for i j.

30

Consider now, the case of the element being subjected to a rigid body motion so that at each node a displacement ui=1 is applied. It is obvious from physical considerations that every point in the element must have displacements u=1. Thus from Eq 3.1 we have ii) Ni = 1. This offers a very simple check for shape functions. Another useful check comes again from considerations of rigid body motion - that a condition of zero strain must be produced when rigid body displacements are applied to the nodes. The reader is encouraged to show that this entails conditions such as iii) Ni, = 0, Ni, = 0 for a 2-dimensional problem. Note that for C elements, where derivatives of are also used as nodal degrees of freedom, these rules apply only to those Ni which are associated with the translational degrees of freedom.
1

3.4 A numerical experiment regarding completeness

Let us now return to the uniform bar carrying an axial load P as shown in Fig. 2.1. The exact solution for this problem is given by a state of constant strain. Thus, the polynomial series u (x ) = a0 + a1x (see Section 2.3) will contain enough terms to ensure completeness. We could at the outset arrange for a0 to be zero to satisfy the essential boundary condition. With u1 (x ) = a1x we were able to obtain the exact solution to the problem. We can now investigate what will happen if we omit the a1x term and start out with the trial function set

u2 (x ) = a2 x 2 . It is left to the reader to show that this results in a computed

stress variation 2 = 3Px 2EAL ; i.e. a linear variation is seen instead of the actual constant state of stress. The reader is now encouraged to experiment with RR solutions
2

based
3

on

other

incomplete

sets

like

u2 +3 (x ) = a2 x

u3 (x ) = a3 x 3

and

+ a3 x . Figure 3.1 shows how erroneous the computed stresses 2, 3

and 2+3 are. The reason for this is now apparent. The term that we have omitted is the one needed to cover the state of constant strain. In this case, this very term also provides the exact answer. The omission of this term has led to the loss of completeness. In fact, the reader can verify that as long as this term is left out of any trial set containing polynomials of higher order, there will not be any sort of clear convergence to the true state of strain. In the present example, because of the absence of the linear term in u (the constant strain term therefore vanishes), the computed strain vanishes at the origin!

31

Fig. 3.1 Stresses from RR solutions based on incomplete polynomials

3.5 Concluding remarks

We have now laid down some specific qualities (continuity and completeness) that trial functions must satisfy to ensure that approximate solutions converge to correct results. It is also hoped that if these requirements are met, solutions from coarse meshes will be reasonably accurate and that convergence will be acceptably rapid with mesh refinement. This was the received wisdom or reigning paradigm around 1977. In the next chapter, we shall look at the twin issues of errors and convergence from this point of view. Unfortunately, elements derived rigorously from only these basic paradigms can behave in unreasonably erratic ways in many important situations. These difficulties are most pronounced in the lowest order finite elements. The problems encountered were called locking, parasitic shear, etc. Some of the problems may have gone unrecorded with no adequate framework or terminology to classify them. As a very good example, for a very long time, it was believed that curved beam and shell elements performed poorly because they could not meet the strain-free rigid body motion condition. However, more recently, the correct error inducing mechanism has been discovered and these problems have come to be called membrane locking. An enlargement of the paradigm is now inevitable to take these strange features into account. This shall be taken up from Chapter 5.

32

Chapter 4 Convergence and Errors

4.1 Introduction
Error analysis is that aspect of fem knowledge, which can be said to belong to the "second-order tradition" of knowledge. To understand its proper place, let us briefly review what has gone before in this text. We have seen structural engineering appear first as art and technology, starting as the imitation of structural form by human fabrication. The science of solid and structural mechanics codified this knowledge in more precise terms. Mathematical models were then created from this to describe the structural behavior of simple and complex systems analytically and quantitatively. With the emergence of cheap computational power, techniques for computational simulation emerged. The finite element method is one such device. It too grew first as art, by systematic practice and refinement. In epistemological terms, we can see this as a first-order tradition or level or inner loop of fem art and practice. Here, we know what to do and how to do it. However, the need for an outer loop or second-order tradition of enquiry becomes obvious. How do we know why we should do it this way? If the first stage is the stage of action, this second stage is now of reflection. Error analysis therefore belongs to this tradition of epistemological enquiry. Let us now translate what we mean by error analysis into simpler terms in the present context. We must first understand that the fem starts with the basic premises of its paradigms, its definitions and its operational procedures to provide numerical results to physical problems which are already described by mathematical models that make analytical quantification possible. The answer may be wrong because the mathematical model wrongly described the physical model. We must take care to understand that this is not the issue here. If the mathematical model does not explain the actual physical system as observed through carefully designed empirical investigations, then the models must be refined or revised so that closer agreement with the experiment is obtained. This is one stage of the learning process and it is now assumed that over the last few centuries we have gone through this stage successfully enough to accept our mathematical models without further doubt or uncertainty. Since our computational models are now created and manipulated using digital computers, there are errors which occur due to the fact that information in the form of numbers can be stored only to a finite precision (word length as it is called) at every stage of the computation. These are called round-off errors. We shall assume here that in most problems we deal with, word length is sufficient so that round-off error is not a major headache. The real issue that is left for us to grapple with is that the computational model prepared to simulate the mathematical model may be faulty and can lead to errors. In the process of replacing the continuum region by finite elements, errors originate in many ways. From physical intuition, we can argue that this will depend on the type and shape of elements we use, the number of elements used and the grading or density of the mesh used, the way distributed loads are assigned to nodes, the manner in which boundary

33

conditions are modeled by specification of nodal degrees of freedom, etc. These are the discretization errors that can occur. Most of such errors are difficult to quantify analytically or determine in a logically coherent way. We can only rely on heuristic judgement to understand how best to minimize errors. However, we shall now look only at that category of discretization error that appears because the computational or discretized model uses trial functions, which are an approximation of the true solution to the mathematical model. It seems possible that to some extent, analytical quantification of these errors is possible. We can recognize two kinds of discretization error belonging to this category. The first kind is that which arises because a model replaces a problem with an infinitely large number of degrees of freedom with a finite number of degrees of freedom. Therefore, except in very rare cases, the governing differential equations and boundary conditions are satisfied only approximately. The second kind of error appears due to the fact that by overlooking certain essential requirements beyond that specified by continuity and completeness, the mathematical model can alter the physics of the problem. In both cases, we must be able to satisfy ourselves that the discretization process, which led to the computational model, has introduced a certain predictable degree of error and that this converges at a predictable rate, i.e. the error is removed in a predictable manner as the discretization is improved in terms of mesh refinement. Error analysis is the attempt to make such predictions a priori, or rationalize the errors in a logical way, a posteriori, after the errors are found. To carry out error analysis, new procedures have to be invented. These must be set apart from the first-order tradition procedures that carry out the discretization (creating the computational model from the mathematical model) and solution (computing the approximate results). Thus, we must design auxiliary procedures that can trace errors in an a priori fashion from basic paradigms (conjectures or guesses). These error estimates or predictions can be seen as consequences computed from our guesses about how the fem works. These errors must now be verified by constructing simple digital computation exercises. This is what we seek to do now. If this cycle can be completed, then we can assure ourselves that we have carved out a scientific basis for error analysis. This is a very crucial element of our study. The fem, or for that matter, any body of engineering knowledge, or engineering methodology, can be said to have acquired a scientific basis only when it has incorporated within itself, these auxiliary procedures that permit its own self-criticism. Therefore, error analysis, instead of being only a posteriori or post mortem studies, as it is usually practised, must ideally be founded on a priori projections computed from intelligent paradigms which can be verified (falsified) by digital computation. In this chapter, we shall first take stock of the conventional wisdom regarding convergence. This is based on the old paradigm that fem seeks to approximate displacements accurately. We next take note of the newly established paradigm that the Ritz-type and fem approximations seek strains/stresses in a `best-fit' manner. From such an interpretation, we examine if it is possible to argue that errors, whether in displacements, stresses or energies, due to finite element discretization must diminish 2 rapidly, at least in a (l/L) manner or better, where a large structure (domain) of dimension L is sub-divided into elements (sub-domains) of dimension l. Thus,

34

with ten elements in a one-dimensional structure, errors must not be more than a few percent. This is the usual range of problems where the continuity and completeness paradigms explain completely the performance of finite elements. In subsequent chapters, we shall discover however that a class of problems exist where errors are much larger - the discretization fail in a dramatic fashion. Convergence and error analysis must now be founded on a more complex conceptual framework - new paradigms need to be introduced and falsified. This will be postponed to subsequent chapters.

4.2 Traditional order of error analysis

The order of error analysis approach that is traditionally used is inherited from finite difference approaches. This seems reasonable because quite often simple finite difference and finite element approximations result in identical equations. It is also possible with suitable interpretations to cast finite difference approximations as particular cases of weighted residual approximations using finite element type trial functions. However, there are some inherent limitations in this approach that will become clear later. It is usual to describe the magnitude of error in terms of the mesh size. Thus, if a series of approximate solutions using grids whose mesh sizes are uniformly reduced is available, it may be possible to obtain more information about the exact answer by some form of extrapolation, provided there is some means to establish the rate at which errors are removed with mesh size. The order of error analysis proceeds from the understanding that the finite element method seeks approximations for the displacement fields. The errors are therefore interpreted using Taylor Series expansions for the true displacement fields and truncations of these to represent the discretized fields. The simple example below will highlight the essential features of this approach.

4.2.1 Error analysis of the axially loaded bar problem

Let us now go back to the case of the cantilever bar subjected to a uniformly distributed axial load of intensity q0 (Section 2.3). The equilibrium equation for this problem is

u,xx + q 0 AE = 0

(4.1)

We shall use the two-node linear bar elements to model this problem. We have seen that this leads to a solution, which gives exact displacements at the nodes. It was also clear to us that this did not mean that an exact solution had been obtained; in fact while the true solution required a quadratic variation of the displacement field u(x), the finite element solution u (x) was piecewise linear. Thus, within each element, there is some error at locations between the nodes. Fig. 4.1 shows the displacement and strain error in the region of an element of length 2l placed with its centroid at xi in a cantilever bar of length L. If for convenience we choose q 0 AE = 1 , then we can show that

u(x) = Lx x 2 2 u (x) = Lx (l
2

(4.2a)
2 xi )

+ 2xxi

(4.2b)

35

(x ) (x )

Fig. 4.1 Displacement and strain error in a uniform bar under distributed axial load q 0

(x)=Lx

(4.2c) (4.2d)

(x) = L xi

If we denote the errors in the displacement field and strain field by e(x) and e(x) respectively, we can show that

e (x ) = (x xi + l ) (x xi l ) 2 e (x ) = (x xi )

(4.3a) (4.3b)

From this, we can argue that in this case, the strain error vanishes at the element centroid, x=xi, and that it is a maximum at the element nodes. This is clear from Fig. 4.1. It is also clear that for this problem, the displacement error is a maximum at x=xi. What is more important to us is to derive measures for these errors in terms of the element size l, or more usefully, in terms of the parameter h = l L , which is the dimensionless quantity that indicates the mesh division relative to the size of the structure. It will also be useful to have these errors normalized with respect to typical displacements and strains occurring in the problem. From Equations (4.2) and (4.3), we can estimate the maximum normalized errors to be of the following orders of magnitude, where O stands for "order",

e u max = O h 2

( )
36

(4.4a)

max

= O (h )

(4.4b)

Note that the order of error we have estimated for the displacement field is that for a location inside the element as in this problem, the nodal deflections are exact. The strain error is O(h) while the displacement error is 2 O(h ). At element centroids, the error in strain vanishes. It is not clear to us in this analysis that the discretized strain (x ) is a "best-fit" of the true strain (x ) . If this is accepted, then it is possible to show that the error in
2

the strain energy stored in such a situation is O(h ) as well. It is possible to generalize these findings in an approximate or tentative way for more complex problems discretized using elements of greater precision. Thus, if an element with q nodes is used, the trial functions are of degree q-1. If the exact solution requires a polynomial field of degree q at q least, then the computed displacement field will have an error O(h ). If strains th for the problem are obtained as the r derivative of the displacement field, q-r then the error in the strain or stress is O(h ). Of course, these measures are tentative, for errors will also depend on how the loads are lumped at the nodes and so on.

4.3 Errors and convergence from "best-fit" paradigm

We have argued earlier that the finite element method actually proceeds to compute strains/stresses in a `best-fit' manner within each element. We shall now use this argument to show how convergence and errors in a typical finite element model of a simple problem can be interpreted. We shall see that a much better insight into error and convergence analysis emerges if we base our treatment on the `best-fit' paradigm. We must now choose a suitable example to demonstrate in a very simple fashion how the convergence of the solution can be related to the best-fit paradigm. The bar under axial load, an example we have used quite frequently so far, is not suitable for this purpose. A case such as the uniform bar with linearly varying axial load modeled with two-node elements gives nodal deflections, which are exact, even though the true displacement field is cubic but between nodes, in each element, only a linear approximate trial function is possible. It can actually be proved that this arises from the special mathematical nature of this problem. We should therefore look for an example where nodal displacements are not exact. The example that is ideally suited for this purpose is a uniform beam with a tip load as shown in Fig. 4.2. We shall model it with linear (two-node) Timoshenko beam elements which represent the bending moment within each element by a constant. Since the bending moment varies linearly over the beam for this problem, the finite element will replace this with a stairstep approximation. Thus, with increase in number of elements, the stress pattern will approach the true solution more closely and therefore the computed strain energy due to bending will also converge. Since the applied load is at the tip, it is very easy to associate this load with the deflection under the load using Castigliano's theorem. It will then be possible to discern the convergence for the tip deflection. Our challenge is therefore to see if the best-fit paradigm can be used to predict the convergence rate for this example from first principles.

37

Fig. 4.2 Cantilever beam under tip load

4.3.1 Cantilever beam idealized with linear Timoshenko beam elements

The dimensions of a cantilever under tip load (Fig. 4.2) are chosen such that the tip deflection under the load will be w=4.0. The example chosen represents a thin beam so that the influence of shear deformation and shear strain energy is negligible. We shall now discretize the beam using equal length linear elements based on Timoshenko theory. We use this element instead of the classical beam element for several reasons. This element serves to demonstrate the features of shear locking which arise from an inconsistent definition of the shear strains (It is therefore useful to take up this element for study later in Chapter 6). After correcting for the inconsistent shear strain, this element permits constant bending and shears strain accuracy within each element - the simplest representation possible under the circumstances and therefore an advantage in seeing how it works in this problem. We shall restrict attention to the bending moment variation as we assume that the potential energy stored is mainly due to bending strain and that we can neglect the transverse shear strain energy for the dimensions chosen. Figure 4.3 shows the bending moment diagrams for a 1, 2 and 4 element idealizations of the present problem using the linear element. The true bending moment (shown by the solid line) varies linearly. The computed (i.e. discretized) bending moments are distributed in a piecewise constant manner as shown by the broken lines. In each case, the elements pick up the bending moment at the centroid correctly - i.e. it is doing so in a `best-fit' manner. What we shall now attempt to do for this problem is to relate this to the accuracy of results. We shall now interpret accuracy in the conventional sense, as that of the deflection at the tip under the load. Table 4.1 shows the

Table 4.1 Normalized tip deflections of a thin cantilever beam, L/t = 100

No. of elements Predicted rate Element without RBF Element with RBF

1 .750 .750 1.000

2 .938 .938 1.000

3 .984 .984 1.000

38

normalized tip deflection with increasing idealizations (neglecting a very small amount due to shear deformation). An interesting pattern emerges. If error is measured by the norm {w w (fem ) } w , it turns out that this is given exactly by the formula 1 4N 2 where N is the number of elements. It can now be seen that this relationship can be established by arguing that this feature emerges from the fact that strains are sought in the `best-fit' manner shown in Fig. 4.3. Consider a beam element of length 2l (Fig. 4.4). Let the moment and shear force at the centroid be M and V. Thus the true bending moment over the element region for our problem can be taken to vary as M+Vx. The discretized bending moment sensed by our linear element would therefore be M. We shall now compute the actual bending energy in the element region (i.e. from a continuum analysis) and that given by the finite element (discretized) model. We can show that

Fig. 4.3 Bending moment diagrams for a one-, two- and four-element dealizations of a cantilever beam under tip load.

39

Energy in continuum model Energy in discretized model

= =

(l EI) (M 2 + V 2l 2 3) (l EI) (M 2 )

(4.5) (4.6)

Thus, as a result of the discretization process involved in replacing each continuum segment of length 2l by a linear Timoshenko beam element which can give only a constant value M for the bending moment, there is a reduction 2 2 (error) in energy in each element equal to (l/EI) (V l /3). It is simple now to show from this that for the cantilever beam of length L with a tip load P, the 2 total reduction in strain energy of the discretized model for the beam is U/4N 2 3 where U=P L /6EI is the energy of the beam under tip load. We are interested now to discover how this error in strain energy translates into an error in the deflection under load P. This can be very easily deduced using Castigliano's second theorem. It is left to the reader to show that the tip deflections of the continuum and discretized model will differ as

{w w (fem ) } w

= 1 4N 2 .

Table 4.1 shows this predicted rate of convergence. Our foregoing analysis shows that this follows from the fact that if any linear variation is approximated in a piecewise manner by constant values as seen in Fig. 4.3, this is the manner in which the square of the error in the stresses/strains (or, equivalently, the difference in work or energy) will converge with idealization. Of course, in a problem where the bending moment is constant, the rate of convergence will be better than this (in fact, exact) and in the case where the bending moment is varying quadratically or at a higher rate, the rate of convergence will be decidedly less.

Fig. 4.4 Bending moment variation in a linear beam element.

40

We also notice that convergence in this instance is from below. This can be deduced from the fact that the discretized potential energy U is less than the actual potential energy U for this problem. It is frequently believed that the finite element displacement approach always underestimates the potential energy and a displacement solution is consequently described as a lower bound solution. However, this is not a universally valid generalization. We can see briefly below (the reader is in fact encouraged to work out the case in detail) where the cantilever beam has a uniformly distributed load acting on it using the same linear Timoshenko beam element for discretization that this is not the case. It turns out that tip rotations converge from above (in a 1 2N 2 rate) while the tip deflections are fortuitously exact. The lower bound solution nature has been disturbed because of the necessity of altering the load system at the nodes of the finite element mesh under the `consistent load' lumping procedure. As promised above, we now extend the concept of "best-fit" and variationally correct rate of convergence to the case of uniformly distributed load of intensity q on the cantilever with a little more effort. Now, when a finite element model is made, two levels of discretization error are introduced. Firstly, the uniformly distributed load is replaced by consistent loads, which are concentrated at element nodes. Thus, the first level of discretization error is due to the replacement of the quadratically varying bending moment in the actual beam with a linear bending moment within each beam element. Over the entire beam model, this variation is piecewise linear. The next level of error is due to the approximation implied in developing the stiffness matrix which we had considered above this effectively senses a "bestapproximated" constant value of the bending moment within each element of the linear bending moment appearing to act after load discretization. With these assumptions, it is a simple exercise using Castigliano's theorem and fictitious tip force and moment P and M respectively to demonstrate that the finite element model of such a problem using two-noded beam elements will yield a fortuitously correct tip deflection w = qL4 8EI for all idealizations (i.e. even with one element!) and a tip rotation that converges at the rate 1 2N 2 from above to the exact value = qL3 6EI . Thus, as far as tip deflections are concerned, the two levels of discretization errors have nicely cancelled each other to give correct answers. This can deceive an unwary analyst into believing that an exact solution has been reached. Inspection of the tip rotation confirms that the solution is approximate and converging. We see from the foregoing analysis that using linear Timoshenko beam elements for the tip loaded cantilever, the energy for this problem converges 2 as O(h ) where h=2l=L/N is the element size. We also see that this order of convergence carries over to the estimate of the tip deflections for this problem. Many text-books are confused over such relationships, especially those that proceed on the order of error analysis. These approaches arrive at conclusions such as strain error is proportional to element size, i.e. O(h) and 2 displacement error proportional to the square of the element size, i.e. O(h ) for this problem. We can see that for this problem (see Fig. 4.3) this estimate is meaningful if we consider the maximum error in strain to occur at element nodes (at centroids the errors are zero as these are optimal strain points). We also see that with element discretization, these errors in strain vanish as

41

O(h). We can also see that the strain energies are now converging at the rate 2 of O(h ) and this emerges directly from the consideration that the discretized strains are `best-fits' of the actual strain. This conclusion is not so readily arrived at in the order of error analysis methods, which often argue that the 2 strains are accurate to O(h), then strain energies are accurate to O(h ) because strain energy expressions contain squares of the strain. This conclusion is valid only for cases where the discretized strains are `best-fit' approximations of the actual strains, as observed in the present example. If the `best-fit' paradigm did not apply, the only valid conclusion that could be drawn is that the strains that have O(h) error will produce errors in strain energy that are O(2h).

4.4 The variationally correct rate of convergence

It is possible now to postulate that if finite elements are developed in a clean and strictly "legal" manner, without violating the basic energy or variational principles, there is a certain rate at which solutions will converge reflecting the fidelity that the approximate solution maintains with the exact solution. We can call this the variationally correct rate of convergence. However, this prescription of strictly legal formulation is not always followed. It is not uncommon to encounter extra-variational devices that are brought in to enhance performance.

4.5 Residual bending flexibility correction

The residual bending flexibility (RBF) correction [4.1,4.2] is a device used to enhance the performance of 2-node beam and 4-node rectangular plate elements (these are linear elements) so that they achieve results equivalent to that 0 obtained by elements of quadratic order. In these C elements (Timoshenko theory for beam and Mindlin theory for plate) there is provision for transverse shear strain energy to be computed in addition to the bending energy (which is the 1 only energy present in the C Euler-Bernoulli beam and Kirchhoff plate formulations). The RBF correction is a deceptively simple device that enhances performance of the linear elements by modifying the shear energy term to compensate for the deficient bending energy term so that superior performance is obtained. To understand that this is strictly an extra-variational trick (or "crime"), it is necessary to understand that a variationally correct procedure will yield linear elements that have a predictable and well-defined rate of convergence. This has been carried out earlier in this chapter. It is possible to improve performance beyond this limit only by resorting to some kind of extra-variational manipulation. By a variationally correct procedure, we mean that the stiffness matrix is derived strictly using a BTDB triple product, where B and D are the strain-displacement and stress-strain matrices respectively.

4.5.1 The mechanics of the RBF correction.

MacNeal [4.2] describes the mathematics of the RBF correction using an original cubic lateral displacement field (corresponding to a linearly varying bending moment field) and what is called a aliased linearly interpolated field (giving a constant bending moment variation for the discretized field). We might recall that we examined the aliasing paradigm in Section 2.5 earlier. We find that a quicker and more physically insightful picture can be obtained by using Equations (4.5) and (4.6) above. Consider a case where only one element of

42

length 2l is used to model a cantilever beam of length L with a tip load P. Applying Castigliano's theorem, we can show very easily that the continuum and discretized solutions will differ by,

w L = P L2 3EI + 1 kAG
2

(4.8) Note the difference in bending flexibilities. This describes the inherent approximation involved in the discretization process if all variational norms are strictly followed. The RBF correction proposes to manipulate the shear flexibility in the discretized case so that it compensates for the deficiency * in the bending flexibility for this problem. Thus if k is the compensated shear correction factor, from Eq (4.7) and (4.8), we have

( w L = P (L

) 4EI + 1 kAG )

(4.7)

L2 3EI + 1 kAG = L2 4EI + 1 k* AG 1 k* AG = 1 kAG + l 2 3EI

(4.9)

It is not accurate to say here that the correction term is derived from the bending flexibility. The bending flexibility of an element that can represent only a constant bending moment (e.g. the two-noded beam element used here) is

L2 4EI . The missing L2 12EI

(or l 2 3EI ) is now brought in as a compensation

*

through the shear flexibility, i.e. k is changed to k . This "fudge factor" therefore enhances the performance of the element by giving it a rate of convergence that is not variationally permissible. Two wrongs make a right here; or do they? Table 4.1 shows how the convergence in the problem shown in Fig. 4.2 is changed by this procedure.

4.6 Concluding remarks.

We have shown that the best-fit paradigm is a useful starting point for deriving estimates about errors and convergence. Using a simple example and this simple concept, we could establish that there is a variationally correct rate of convergence for each element. This can be improved only by taking liberties with the formulation, i.e. by introducing extra-variational steps. The RBF is one such extra-variational device. It is also possible to argue that errors, whether in displacements, stresses or energies, due to finite element discretization must converge rapidly, at least 2 in a O(h ) manner or better. If a large structure (domain) of dimension L is sub-divided into elements (sub-domains) of dimension l, one expects errors of the order of (l L )2 . Thus, with ten elements in a one-dimensional structure, errors must not be more than a few percent. We shall discover however that a class of problems exist where errors are much larger - the discretizations fail in a dramatic fashion, and this cannot be resolved by the classical (pre-1977) understanding of the finite element method. A preliminary study of the issues involved will be taken up in the next chapter; the linear Timoshenko beam element serves to expose the factors clearly. Subsequent chapters will undertake a comprehensive review of the various manifestations of such errors. It is felt that this detailed treatment is justified, as an understanding of

43

such errors has been one of the most challenging problems that the finite element procedure has faced in its history. Most of the early techniques to overcome these difficulties were ad hoc, more in the form of `art' or `black magic'. In the subsequent chapters, our task will be to identify the principles that establish the methodology underlying the finite element procedure using critical, rational scientific criteria.

4.7 References
4.1. R. H. MacNeal, A simple quadrilateral plate element, Comput. Struct., 8, 175-183, 1978. 4.2. R. H. MacNeal, Finite Elements: Their Design and Performance, Marcel Dekker, NY, 1994.

44

Chapter 5 The locking phenomena

5.1 Introduction
A pivotal area in finite element methodology is the design of robust and accurate elements for applications in general purpose packages (GPPs) in structural analysis. The displacement type approach we have been examining so far is the technique that is overwhelmingly favored in general purpose software. We have seen earlier that the finite element method can be interpreted as a piecewise variation of the Rayleigh-Ritz method and therefore that it seeks strains/stresses in a `best-fit' manner. From such an interpretation, it is possible to argue that errors, whether in displacements, stresses or energies, due to finite element discretization must manner or better, where a large structure into elements (sub-domains) of dimension dimensional structure, errors must not be converge rapidly, at least in a (l L )2 (domain) of dimension L is sub-divided l. Thus, with ten elements in a onemore than a few percent.

By and large, the elements work without difficulty. However, there were spectacular failures as well. These are what are now called the locking 0 problems in C finite elements - as the element libraries of GPPs stabilized these elements came to be favored for reasons we shall discuss shortly. By locking, we mean that finite element solutions vanish quickly to zero (errors saturating quickly to nearly 100%!) as certain parameters (the penalty multipliers) become very large. It was not clear why the displacement type method, as it was understood around 1977, should produce for such problems, answers that were only a fraction of a per cent of the correct answer with a practical level of discretization. Studies in recent years have established that an aspect known as consistency must be taken into account. The consistency paradigm requires that the interpolation functions chosen to initiate the discretization process must also ensure that any special constraints that are anticipated must be allowed for in a consistent way. Failure to do so causes solutions to lock to erroneous answers. The paradigm showed how elements can be designed to be free of these errors. It also enabled error-analysis procedures that allowed errors to be traced to the inconsistencies in the representation to be developed. We can now develop a family of such error-free robust elements for applications in structural mechanics. In this chapter, we shall first introduce the basic concepts needed to understand why such discretized descriptions fail while others succeed. We 0 compare the equations of the Timoshenko beam theory (a C theory) to the 1 classical beam theory (a C theory) to show how constraints are generated in such a model. This permits us to discuss the concept of consistency and the nature of errors that appear during a Ritz type approximation. These same errors are responsible for the locking seen in the displacement type finite element models of similar problems.

45

5.2 From C1 to C0 elements

As the second generation of GPPs started evolving around the late 70s and early 1 80s, their element libraries replaced what were called the C elements with what 0 were known as the C elements. The former were based on well known classical theories of beams, plates and shells, reflecting the confidence structural analysts had in such theories for over two centuries - namely the EulerBernoulli beam theory, the Kirchhoff-Love plate theory and the equivalent shell theories. These theories did not allow for transverse shear strain and permitted the modeling of such structures by defining deformation in terms of a single field, w, the transverse deflection of a point on what is called the neutral axis (in a beam) and neutral surface of a plate or shell. The strains could then be computed quite simply from the assumption that normal to the neutral surface remained normal after deformation. One single governing differential equation resulted, although of a higher order (in comparison to other theories we shall discuss shortly), and this was considered to be an advantage. There were some consequences arising from such an assumption both for the mathematical modeling aspect as well as for the finite element (discretization) aspect. In the former, it turned out that to capture the physics of deformation of thick or moderately thick structures, or the behavior of plates and shells made of newly emerging materials such as high performance laminated composites, it was necessary to turn to more general theories accounting or transverse shear deformation as well - these required the definition of rotations of normals which were different from the slopes of the neutral surface. Some of the 1 contradictions that arose as a result of the old C theories - e.g. the use of the fiction of the Kirchhoff effective shear reactions, could now be removed, restoring the more physically meaningful set of three boundary conditions on the edge of a plate or shell (the Poisson boundary conditions as they are called) to be used. The orders of the governing equations were correspondingly reduced. A salutary effect that carried over to finite element modeling was that the elements could be designed to have nodal degrees of freedom which were the six basic engineering degrees of freedom - the three translations and the three rotations at a point. This was ideal from the point of view of the organization of a general purpose package. Also, elements needed only simple basis functions requiring only the continuity of the fields across element boundaries - these 0 1 are called the C requirements. In the older C formulations, continuity of slope was also required and to achieve this in arbitrarily oriented edges, as would be found in triangular or quadrilateral planforms of a plate bending element, it was necessary to retain curvature degrees of freedom (w,xx, w,xy, w,yy) at the nodes and rely on quintic polynomials for the element shape or basis functions. So, as general purpose packages ideal for production run 0 analyses and design increasingly found favour in industry, the C beam, plate 1 and shell elements slowly began to replace the older C equivalents. It may be instructive to note that the general two-dimensional (i.e. plane stress, plane strain and axisymmetric) elements and three-dimensional (solid or brick as they 0 are called) elements were in any case based on C shape functions - thus this 0 development was welcome in that universally valid C shape functions and their derivatives could be used for a very wide range of structural applications.

46

However, life was not very simple - surprisingly dramatic failures came to be noticed and the greater part of academic activity in the late seventies, most of the eighties and even in the nineties was spent in understanding and eliminating what were called the locking problems.

5.3 Locking, rank and singularity of penalty-linked stiffness matrix, and consistency of strain-field
When locking was first encountered, efforts were made to associate it with the rank or non-singularity of the stiffness matrix linked to the penalty term (e.g. the shear stiffness matrix in a Timoshenko beam element which becomes very large as the beam becomes very thin, see the discussion below). However, on reflection, it is obvious that these are symptoms of the problem and not the cause. The high rank and non-singularity is the outcome of certain assumptions made (or not made, i.e. leaving certain unanticipated requirements unsatisfied) during the discretization process. It is therefore necessary to trace this to the origin. The consistency approach argues that it is necessary in such problems to discretize the penalty-linked strain fields in a consistent way so that only physically meaningful constraints appear. In this section, we would not enter into a formal finite element discretization (which would be taken up in the next section) but instead, illustrate the concepts involved using a simple Ritz-type variational method of approximation of the beam problem via both classical and Timoshenko beam theory [5.1]. It is possible to show how the Timoshenko beam theory can be reduced to the classical thin beam theory by using a penalty function interpretation and in doing so, show how the Ritz approximate solution is very sensitive to the way in which its terms are chosen. An `inconsistent' choice of parameters in a low order approximation leads to a full-rank (non-singular) penalty stiffness matrix that causes the approximate solution to lock. By making it `consistent', locking can be eliminated. In higher order approximations, `inconsistency' does not lead to locked solutions but instead, produces poorer convergence than would otherwise be expected of the higher order of approximation involved. It is again demonstrated that a Ritz approximation that ensures ab initio consistent definition will produce the expected rate of convergence - a simple example will illustrate this.

5.3.1 The classical beam theory

Consider the transverse deflection of a thin cantilever beam of length L under an uniformly distributed transverse load of intensity q per unit length of the beam. This should produce a linear shear force distribution increasing from 0 at the free end to qL at the fixed end and correspondingly, a bending moment that varies from 0 to qL2 2 . Using what is called the classical or Euler-Bernoulli theory, we can state this problem in a weak form with a quadratic functional given by,

2 (1 2 EI w,xx qw ) dx 0 L

(5.1)

This theory presupposes a constraint condition, assuming zero transverse shear strain, and this allows the deformation of a beam to be described entirely

47

in terms of a single dependent variable, the transverse deflection w of points on the neutral axis of the beam. An application of the principle of minimum total potential allows the governing differential equations and boundary conditions to be recovered, but this will not be entered into here. A simple exercise will establish that the exact solution satisfying the governing differential equations and boundary conditions is,

w (x ) = qL2 4EI x 2 (qL 6EI ) x 3 (q 24EI ) x 4 EI w,xx = (q 2 ) (L x )2 EI w,xxx = q (L x )

(5.2a) (5.2b) (5.2c)

5.3.1.1 A two-term Ritz approximation

Let us
2

now
3

consider

an

approximate

Ritz

solution

based

on

two

terms,

w = b2 x + b3 x . Note that the constant and linear terms are dropped, anticipating the boundary conditions at the fixed end. One can easily work out that the approximate solution will emerge as, w (x ) = 5qL2 24EI x 2 (qL 12EI ) x 3

(5.3a)

so that the approximate bending moment and shear force determined in this Ritz process are,

EI w ,xx = 5qL2 12 (qL 2 ) x EI w ,xxx = qL 2

(5.3b) (5.3c)

If the expressions in Equations (5.2) and (5.3) are written in terms of the natural co-ordinate system , where x = (1 + ) L 2 so that the ends of the beam are represented by = -1 and +1, the exact and approximate solutions can be expanded as,

EI w,xx = qL2 8 4 3 2 1 3 (1 3 2 ) EI w ,xx = qL2 8

)(

(5.4) (5.5)

) (4 3 2 )

The approximate solution for the bending moment is seen to be a `best-fit' or `least-squares fit' of the exact solution, with the points = 1 3 , which are the points where the second order Legendre polynomial vanishes, emerging as points where the approximate solution coincides with the exact solution. From Equations (5.2c) and (5.3c), we see that the shear force predicted by the Ritz approximation is a 'best-fit' of the actual shear force variation. Once again we confirm that the Ritz method seeks the 'best-approximation' of the actual state of stress in the region being studied.

48

5.3.2 The Timoshenko beam theory and its penalty function form
The Timoshenko beam theory [5.1] offers a physically more general formulation of beam flexure by taking into account the transverse shear deformation. The description of beam behavior is improved by introducing two quantities at each point on the neutral axis, the transverse deflection w and the face rotation so that the shear strain at each point is given by = w,x , the difference between the face rotation and the slope of the neutral axis. The total strain energy functional is now constructed from the two independent functions for w(x) and (x), and it will now account for the bending (flexural) energy and an energy of shear deformation.

L 2 2 (1 2 EI ,x + 1 2 ( w,x ) q w ) dx 0

(5.6)

where the curvature =,x, the shear strain =-w,x and =kGA is the shear rigidity. The factor k accounts for the averaged correction made for the shear strain distribution through the thickness. For example, for a beam of rectangular cross-section, this is usually taken as 5/6. The Timoshenko beam theory will asymptotically recover the elementary beam theory as the beam becomes very thin, or as the shear rigidity becomes very large, i.e. . This requires that the Kirchhoff constraint -w,x0 must emerge in the limit. For a very large , these equations lead directly to the simple fourth order differential equation for w of elementary beam theory. Thus, this is secured very easily in the infinitesimal theory but it is this very same point that poses difficulties when a simple Ritz type approximation is made.

5.3.2.1 A two-term Ritz approximation

Consider now a two term Ritz approximation based on = a1x and w = b1x . This permits a constant bending strain (moment) approximation to be made. The shear strain is now given by

= a1x b1

(5.7)

and it would seem that this can represent a linearly varying shear force. The Ritz variational procedure now leads to the following set of equations,

L3 3 L2 2 a1 0 L 0 = 2 + EI 2 L 2 L b1 qL 2 0 0
Solving, we get

(5.8)

a1 =

3qL2 12EI + L2

; b1 =

qL 1.5qL3 2 12EI + L2

(5.9)

49

As , both a1 and b10. This is a clear case of a solution `locking'. This could have been anticipated from a careful examination of Equations (5.6) and (5.7). The penalty limit in Equation (5.6) introduces a penalty condition on the shear strain and this requires that the shear strain must vanish in the Ritz approximation process - from Equation (5.7), the conditions emerging are a10 and b10. Clearly, a10 imposes a zero bending strain ,x0 as well this spurious constraint produces the locking action on the solution. Thus, a meaningful approximation can be made for a penalty function formulation only if the penalty linked approximation field is consistently modeled. We shall see how this is done next.

5.3.2.2 A three-term consistent Ritz approximation

Consider now a three term approximation chosen to satisfy what we shall call the consistency condition. Choose = a1x and w = b1x + b2 x 2 . Again, only a constant bending strain is permitted. But we now have shear strain of linear order as,

= b1 + (a1 2b2 )x

(5.10)

Note that as , the constraint 1-2b20 is consistently balanced and will not result in a spurious constraint on the bending field. The Ritz variational procedure leads to the following equation structure:

L3 3 L2 2 L 0 0 2 L EI 0 0 0 + L 2 3 0 0 0 2L 3 L2

2L3 3 L2 4L3 3

0 a1 2 b1 = qL 2 3 b qL 3 2

(5.11)

It can be easily worked out from this that the approximate solution is given by,

qL2 x 6EI

(5.12a)

w =

qL q 2 qL2 x x + x2 2 12EI EI ,x = qL2 6

(5.12b)

(5.12c) (5.12d)

w ,x = q (L x )

There is no locking seen at all - the bending moment is now a least squares correct constant approximation of the exact quadratic variation (this can be seen by comparing Equation (5.12c) with Equations (5.4) and (5.5) earlier). The shear force is now correctly captured as a linear variation - the consistently represented field in Equation (5.12) being able to recover this even as a !

50

A comparison of the penalty linked matrices in Equations (5.8) and (5.11) shows that while in the former, we have a non-singular matrix of rank 2, in the latter we have a singular matrix of rank 2 for a matrix of order 3. It is clear also that the singularity (or reduced rank) is a result of the consistent condition represented by (1-2b2) in the linear part of the shear strain definition in Equation (5.10) - as a direct consequence, the third row of the penalty linked matrix in Equation (5.11) is exactly twice the first row. It is this aspect that led to a lot of speculation on the role the singularity or rank of the matrix plays in such problems. We can further show that non-singularity of penalty linked matrix arises in an inconsistent formulation only when the order of approximation is low, as seen for the two-term Ritz approximation. We can go for a quadratic inconsistent approximation (with linear bending strain variation) in the next section to show that there is no locking' of the solution and that the penalty linked matrix is not non-singular - however the effect of inconsistency is to reduce the performance of the approximation to a sub-optimal level.

5.3.2.3 A four-term inconsistent Ritz approximation

We now take up a four term approximation which provides theoretically for a linear variation in the approximation for bending strain, i.e. = a1x + a2 x 2 and

w = b1x + b2 x 2 so that = ,x = a1 + 2a2 x and the shear strain is

= b1 + (a1 2b2 )x + a2 x 2

(5.13)

Note now that the condition forces a20 this becomes a spurious constraint on the bending strain field. We shall now see what effect this spurious constraint has on the approximation process. The Ritz variational procedure leads to the following set of equations:

L 2 L EI 0 0

L2 4L3 3 0 0

0 0 0 0

0 0 + 0 0
(5.14)

L3 3 L4 4 L2 2 3 2 L 3

L4 4 L5 5 L3 3 L4 2

L2 2 L3 3 L L2

2L3 3 L4 2 L2 4L3 3

0 a1 a2 0 = 2 b1 qL 2 b2 3 qL 3

It can be seen that the penalty linked 44 matrix is singular as the fourth row is exactly twice the second row - this arises from the consistent representation (a1 2b2 ) of the linear part of shear strain in Equation (5.13). The rank of the

51

matrix is therefore 3 and the solution should be free of "locking" - however the inconsistent constraint forces a2 0 and this means that the computed bending strain a1 ; i.e. it will have only a constant bending moment prediction capability. What it means is that this four term inconsistent approach will produce answers only as efficiently as the three term consistent Ritz formulation. Indeed, one can work out from Equation (5.14) that,

a1 =

qL2 15qL2 6EI 60EI + L2 15qL2 60EI + L2

(5.15a)

a2 =

(5.15b)

b1 =

qL 5qL3 2 60EI + L2

(5.15c)

b2 =

qL2 q + 6EI 2

(5.15d)

As , the bending moment and shear force variation are given by,

EI ,x = qL2 6

(5.16a)

w ,x = q (L x ) + 2.5qL (6 (x L )2 6 (x L ) + 1

)
x L = 1 2 1+1

(5.16b)

The solution can sense only a constant bending moment in the thin beam limit. There are now violent quadratic oscillations in the shear force and these oscillations

1 2 1 1 3 , or = 1 3 , which are the points corresponding to the 2 point Gauss-Legendre integration rule. The effect of the inconsistent representation has been to reduce the effectiveness of the approximation. We shall next see how the effectiveness can be improved by making the approximation consistent before the variational process is carried out.

can

be

shown

to

vanish

at

the

points

and

5.3.2.4 A four-term consistent Ritz approximation

Let us now take up a Ritz approximation with a consistently represented function

. This can be achieved by noting that in Equation (5.13), the offending term is the quadratic term associated with a2. We
for the shear strain defined as also see from Equation (5.16b) that this leads to a spuriously excited quadratic form

(6x

6Lx + L2 . Our experience with consistent finite element formulations

2

[5.2] allows us to replace the x

term in Equation (5.13) with Lx L2 6 so that,

52

 = b1 + (a1 2b2 )x + a2 Lx L2 6
= b1 a2 L2 6

)
(5.17)

)+

(a1 + a2 L - 2b2 )x

In fact, it can be proved using the generalized (mixed) variational theorem such as the Hellinger-Reissner and the Hu-Washizu principles [5.2] that the variationally correct way to determine from the inconsistent is,

T dx = 0

and this will yield precisely the form represented in Equation (5.17). The Ritz variational procedure now leads to,

L 2 L EI 0 0

L2 4L3 3 0 0

0 0 0 0

0 0 + 0 0
(5.18)

L3 3 L4 4 L2 2 3 2 L 3

L4 4 7L5 36 L3 3 L4 2

L2 2 L3 3 L L2

2L3 3 L4 2 L2 4L3 3

0 a1 0 a2 = 2 qL 2 b1 3 b2 qL 3

It is important to recognize now that since the penalty linked matrix emerges from the terms - b1 + a2 L2 6 and (a1 + a2 L - 2b2 ) there would only be two linearly independent rows and therefore the rank of the matrix is now 2. The approximate solution is then given by,

EI =

5 qL qL2 x 12 2

(5.19a) (5.19b)

= q (L - x )

Comparing Equation (5.19a) with Equation (5.3b) we see that the approximate solution to the Timoshenko equation for this problem with a consistent shear strain assumption gives exactly the same bending moment as the Ritz solution to the classical beam equation could provide. We also see that Equation (5.19b) is identical to Equation (5.12b) so that the shear force is now being exactly computed. In other words, as , the terms a1, a2, b1 and b1 all yield physically meaningful conditions representing the state of equilibrium correctly.

53

5.4 Consistency and C0 displacement type finite elements

We shall see in the chapters to follow that locking, poor convergence and 0 displacement type finite element violent stress oscillations seen in C formulations are due to a lack of consistent definition of the critical strain fields when the discretization is made - i.e. of the strain-fields that are constrained in the penalty regime. The foregoing analysis showed how the lack of consistency translates into a non-singular matrix of full rank that causes locking in low-order Ritz approximations of such problems. It is also seen that in higher order approximations, the situation is not as dramatic as to be described as locking, but is damaging as to produce poor convergence and stress oscillations. It is easy to predict all this by examining the constrained strain-field terms from the consistency point of view rather than performing a post-mortem examination of the penalty-linked stiffness matrix from rank or singularity considerations as is mostly advocated in the literature. We shall now attempt a very preliminary heuristic definition of these requirements as seen from the point of view of the developer of a finite element for an application in a constrained media problem. We shall see later that if a simple finite element model of the Timoshenko beam is made), the results are in very great error and that these errors grow without limit as the beam becomes very thin. This is so even when the shape functions for the w and have been chosen to satisfy the completeness and continuity conditions. We saw in our Ritz approximation of the Timoshenko beam theory in this section that as the beam becomes physically thin, the shear strains must vanish and it must begin to enforce the Kirchhoff constraint and that this is not possible unless the approximate field can correctly anticipate this. In the finite element statement of this problem, the shape functions chosen for the displacement fields cannot do this in a meaningful manner - spurious constraints are generated which cause locking. The consistency condition demands that the discretized strain field interpolations must be so constituted that it will enforce only physically true constraints when the discretized functionals for the strain energy of a finite element are constrained. We can elaborate on this definition in the following descriptive way. In the development of a finite element, the field variables are interpolated using interpolations of a certain order. The number of constants used will depend on the number of nodal variables and any additional nodeless variables (those corresponding to bubble functions). From these definitions, one can compute the strain fields using the strain-displacement relations. These are obtained as interpolations associated with the constants that were introduced in the field variable interpolations. Depending on the order of the derivatives of each field variable appearing in the definition of that strain field (e.g. the shear strain in a Timoshenko theory will have and the first derivative of w), the coefficients of the strain field interpolations may have constants from all the contributing field variable interpolations or from only one or some of these. In some limiting cases of physical behavior, these strain fields can be constrained to be zero values, e.g. the vanishing shear strain in a thin Timoshenko beam. Where the strain-field is such that all the terms in it (i.e. constant, linear, quadratic, etc.) have, associated with it, coefficients from all the independent interpolations of the field variables that appear in the definition of that

54

strain-field, the constraint that appears in the limit can be correctly enforced. We shall call such a representation field-consistent. The constraints thus enforced are true constraints. Where the strain-field has coefficients in which the contributions from some of the field variables are absent, the constraints may incorrectly constrain some of these terms. This fieldinconsistent formulation is said to enforce additional spurious constraints. For simple low order elements, these constraints are severe enough to produce solutions that rapidly vanish - causing what is often described as locking.

5.5 Concluding remarks

These exercises show us why it is important to maintain consistency of the basis functions chosen for terms in a functional, which are linked, to penalty multipliers. The same conditions are true for the various finite element formulations where locking, poor convergence and stress oscillations are known to appear. It is also clear why the imposition of the consistency condition into the formulation allows the correct rank or singularity of the penalty linked stiffness matrix to be maintained so that the system is free of locking or suboptimal convergence. Again, it is worthwhile to observe that non-singularity of the penalty linked matrix occurs only when the approximate fields are of very low order as for the two-term Ritz approximation. In higher order inconsistent formulations, as for the four-term inconsistent Ritz approximation, solutions are obtained which are sub-optimal to solutions that are possible if the formulation is carried out with the consistency condition imposed a priori. We shall see later that the use of devices such as reduced integration permits the consistency requirement to be introduced when the penalty linked matrix is computed so that the correct rank which ensures the imposition of the true constraints only is maintained. In the next chapter, we shall go to the locking and other over-stiffening phenomena found commonly in displacement finite elements. The phenomenon of shear locking is the most well known - we shall investigate this closely. The Timoshenko beam element allows the problem to be exposed and permits a mathematically rigorous error model to be devised. The consistency paradigm is introduced to show why it is essential to remove shear locking. The correctness concept is then brought in to ensure that the devices used to achieve a consistent strain interpolation are also variationally correct. This example serves as the simplest case that could be employed to demonstrate all the principles involved in the consistency paradigm. It sets the scene for the need for new paradigms (consistency and correctness) to complement the existing ones (completeness and continuity) so that the displacement model can be scientifically founded. The application of these concepts to Mindlin plate elements is also reviewed. The membrane locking phenomenon is investigated in Chapter 7. The simple curved beams are used to introduce this topic. This is followed up with studies of parasitic shear and incompressible locking. In Chapters 5 to 7 we investigate the effect of constraints on strains due to the physical regimes considered, e.g. vanishing shear strains in thin Timoshenko beams or Mindlin plates, vanishing membrane strains in curved beams and shells etc. In Chapter 8 we proceed to a few problems where consistency between strain and stress functions are needed even where no constraints on the

55

strains appear. These examples show the universal need of the consistency aspect and also the power of the general Hu-Washizu theorem to establish the complete, correct and consistent variational basis for the displacement type finite element method.

5.6 References
5.1 S. P. Timoshenko, On the transverse vibrations of bars of uniform crosssection, Philosophical Magazine 443, 125-131, 1922. 5.2 G. Prathap, The Finite Element Academic Press, Dordrecht, 1993.

Method

in

Structural

Mechanics,

Kluwer

56

Chapter 6 Shear locking

In the preceding chapter, we saw a Ritz approximation of the Timoshenko beam problem and noted that it was necessary to ensure a certain consistent relationship between the trial functions to obtain accurate results. We shall now take up the finite element representation of this problem, which is essentially a piecewise Ritz approximation. Our conclusions from the preceding chapter would therefore apply to this as well.

6.1 The linear Timoshenko beam element

An element based on elementary theory needs two nodes with 2 degrees of freedom at each node, the transverse deflection w and slope dw/dx and uses cubic 1 interpolation functions to meet the C continuity requirements of this theory (Fig. 6.1). A similar two-noded beam element based on the shear flexible 0 Timoshenko beam theory will need only C continuity and can be based on simple linear interpolations. It was therefore very attractive for general purpose applications. However, the element was beset with problems, as we shall presently see.

6.1.1 The conventional formulation of the linear beam element

The strain energy of a Timoshenko beam element of length 2l can be written as the sum of its bending and shear components as,
T T 1 2 EI + 1 2 kGA dx

(6.1) (6.2a) (6.2b)

where

= ,x = w,x

In Equations (6.2a) and (6.2b), w is the transverse displacement and the section rotation. E and G are the Young's and shear moduli and the shear correction factor used in Timoshenko's theory. I and A are the moment of inertia and the area of cross-section, respectively.

(a) 1 w w,x (b) 1 w

2 w w,x 2 w

Fig. 6.1 (a) Classical thin beam and (b) Timoshenko beam elements.

57

In the conventional procedure, displacement field variables as,

linear

interpolations

are

chosen

for

the

N 1 = (1 ) 2

N 2 = (1 + ) 2

(6.3a) (6.3b)

where the dimensionless coordinate =x/l varies from -1 to +1 for an element of length 2l. This ensures that the element is capable of strain free rigid body motion and can recover a constant state of strain (completeness requirement) and that the displacements are continuous within the element and across the element boundaries (continuity requirement). We can compute the bending and shear strains directly from these interpolations using the strain gradient operators given in Equations (6.2a) and (6.2b). These are then introduced into the strain energy computation in Equation (6.1), and the element stiffness matrix is calculated in an analytically or numerically exact (a 2 point Gauss Legendre integration rule) way. For the beam element shown in Fig. 6.1, for a length h the stiffness matrix can be split into two parts, a bending related part and a shear related part, as,

kb

0 0 = EI h 0 0

h2 1 0 h 2 h 2 3 1 0 1 ks = kGt h -1 - h 2 0 0 0 2 1 0 1 h 2 h 6 0 0

h 2 - h 2 h2 6 1 -h 2 2 -h 2 h 3

We shall now model a cantilever beam under a tip load using this element, considering the case of a "thin" beam with E=1000, G=37500000, t=1, L=4, using a fictitiously large value of G to simulate the "thin" beam condition. Table 6.1 shows that the normalized tip displacements are dramatically in error. In fact with a classical beam element model, exact answers would have been obtained with one element for this case. We can carefully examine Table 6.1 to see the trend as the number of elements is increased. The tip deflections obtained, which are several orders of magnitude lower than the correct answer, are directly related to the square of the number of elements used for the idealization. In other words, the discretization process has introduced an error so large that the 2 resulting answer has a stiffness related to the inverse of N . This is clearly unrelated to the physics of the Timoshenko beam and also not the usual sort of discretization errors encountered in the finite element method. It is this very phenomenon that is known as shear locking. Table 6.1 - Normalized tip deflections

No. of elements 1 2 4 8 16

Thin beam -5 0.200 10 -5 0.800 10 -4 0.320 10 -3 0.128 10 -3 0.512 10

58

The error in each element must be related to the element length, and therefore when a beam of overall length L is divided into N elements of equal length h, the additional stiffening introduced in each element due to shear 2 locking is seen to be proportional to h . In fact, numerical experiments showed that the locking stiffness progresses without limit as the element depth t decreases. Thus, we now have to look for a mechanism that can explain how this 2 spurious stiffness of (h/t) can be accounted for by considering the mathematics of the discretization process. The magic formula proposed to overcome this locking is the reduced integration method. The bending component of the strain energy of a Timoshenko beam element of length 2l shown in Equation (6.1) is integrated with a one-point Gaussian rule as this is the minimum order of integration required for exact evaluation of this strain energy. However, a mathematically exact evaluation of the shear strain energy will demand a two-point Gaussian integration rule. It is this rule that resulted in the shear stiffness matrix of rank two that locked. An experiment with a one-point integration of the shear strain energy component causes the shear related stiffness matrix to change as shown below. The performance of this element was extremely good, showing no signs of locking at all (see Table 4.1 for a typical convergence trend with this element).

kb

0 0 = EI h 0 0

0 0 0 1 0 1 0 0 0 1 0 1

h 2 h 2 1 1 h 2 h 2 4 h 2 h 2 4 ks = kGt h 1 h 2 1 h 2 2 2 h 2 h 4 h 2 h 4

6.1.2 The field-consistency paradigm

It is clear from the formulation of the linear Timoshenko beam element using exact integration (we shall call it the field-inconsistent element) that ensuring the completeness and continuity conditions are not enough in some problems. We shall propose a requirement for a consistent interpolation of the constrained strain fields as the necessary paradigm to make our understanding of the phenomena complete. If we start with linear trial functions for w and , as we had done in Equation 6.3 above, we can associate two generalized displacement constants with each of the interpolations in the following manner

w = a0 + a1

(x l ) (x l )

(6.4a) (6.4b)

= b0 + b1

We can relate such constants to the field-variables obtaining in this element in a discretized sense; thus, a1/l=w,x at x=0, b0= and b1/l=,x at x=0. This denotation would become useful when we try to explain how the discretization process can alter the infinitesimal description of the problem if the strain fields are not consistently defined.

59

If the strain-fields are now derived from the displacement fields given in Equation (6.4), we get

= (b1 l ) = (b0 a1 l ) + b1 (x l )

(6.5a) (6.5b)

An exact evaluation of the strain energies for an element of length h=2l will now yield the bending and shear strain energy as

Us

{(b1 l )}2 2 = 1 2 (kGA ) (2l ) {(b0 a1 l )2 + 1 3 b1 }

U B = 1 2 (EI ) (2l )

(6.6a) (6.6b)

It is possible to see from this that in the constraining physical limit of a very thin beam modeled by elements of length 2l and depth t, the shear strain energy in Equation (6.6b) must vanish. An examination of the conditions produced by these requirements shows that the following constraints would emerge in such a limit

b0 a1 l 0 b1 0

(6.7a) (6.7b)

In the new terminology that we had cursorily introduced in Section 5.4, constraint (6.7a) is field-consistent as it contains constants from both the contributing displacement interpolations relevant to the description of the shear strain field. These constraints can then accommodate the true Kirchhoff constraints in a physically meaningful way, i.e. in an infinitesimal sense, this is equal to the condition (-w,x)0 at the element centroid. In direct contrast, constraint (6.7b) contains only a term from the section rotation . A constraint imposed on this will lead to an undesired restriction on . In an infinitesimal sense, this is equal to the condition ,x0 at the element centroid (i.e. no bending is allowed to develop in the element region). This is the `spurious constraint' that leads to shear locking and violent disturbances in the shear force prediction over the element, as we shall see presently.

6.1.3 An error model for the field-consistency paradigm

We must now determine that this field-consistency paradigm leads us to an accurate error prediction. We know that the discretized finite element model will contain an error which can be recognized when digital computations made with these elements are compared with analytical solutions where available. The consistency requirement has been offered as the missing paradigm for the errorfree formulation of the constrained media problems. We must now devise an operational procedure that will trace the errors due to an inconsistent representation of the constrained strain field and obtain precise a priori measures for these. We must then show by actual numerical experiments with the original elements that the errors are as projected by these a priori error models. Only such an exercise will complete the scientific validation of the consistency paradigm. Fortunately, a procedure we shall call the functional reconstitution technique makes it possible to do this verification.

60

6.1.4 Functional re-constitution

We have postulated that the error of shear constraint in Equation (6.7b). case where the inconsistent element depth t. The strain energy for such a shear locking originates from the spurious We must now devise an error model for the is used to model a beam of length L and beam can be set up as,

L {1 2 EI 0

,2 + 1 2 kGA ( w,x )2 dx x

(6.8)

If an element of length 2l is isolated, the discretization process produces energy for the element of the form given in Equation (6.6). In this equation, the constants, which were introduced due to the discretization process, can be replaced by the continuum (i.e. the infinitesimal) description. Thus, we note that in each element, the constants in Equations (6.6a) and (6.6b) can be traced to the constants in Equations (6.4a) and (6.4b) and can be replaced by the values of the field variations , ,x and w,x at the centroid of the element. Thus, the strain energy of deformation in an element is,

e = 1 2 (EI ) (2l ) ( ,x )2 + 1 2 (kGA ) (2l ) ( w,x )2 + 1 6 (kGAl )2 ( ,x )2

(6.9)

Thus the constants in the discretized strain energy functional have been reconstituted into an equivalent continuum or infinitesimal form. From this reconstituted functional, we can argue that an idealization of a beam region of length 2l into a linear displacement type finite element would produce a modified strain energy density within that region of,

e = 1 2 EI + kGAl2 3

) (,

)2

+ 1 2 (kGA ) ( w,x )2

(6.10)

This strain energy density indicates that the original physical system has been altered due to the presence of the inconsistent term in the shear strain field. Thus, we can postulate that a beam of length L modeled by equal elements of length 2l will have a re-constituted functional

2 2 2 {1 2 (EI + kGAl 3 ) (,x ) + 1 2 (kGA ) ( w,x ) }dx 0 L

(6.11)

We now understand that the discretized beam is stiffer in bending (i.e. its flexural rigidity) by the factor kGAl2 3EI . For a thin beam, this can be very large, and produces the additional stiffening effect described as shear locking.

6.1.5 Numerical experiments to verify error prediction

Our functional re-constitution procedure (note that this is an auxiliary procedure, distinct from the direct finite element procedure that yields the stiffness matrix) allows us to critically examine the consistency paradigm. It indicates that an exactly-integrated or field-inconsistent finite element model tends to behave as a shear flexible beam with a much stiffened flexural rigidity I. This can be related to the original rigidity I of the system by comparing the expressions in Equations (6.8) and (6.11) as,

I I = 1 + kGAL2 3EI

(6.12)

61

We must now show through a numerical experiment that this estimate for the error, which has been established entirely a priori, starting from the consistency paradigm and introducing the functional re-constitution technique, anticipates very accurately, the behavior of a field-inconsistent linearly interpolated shear flexible element in an actual digital computation. Exact solutions are available for the static deflection W of a Timoshenko cantilever beam of length L and depth t under a vertical tip load. If W fem is the result from a numerical experiment involving a finite element digital computation using elements of length 2l, the additional stiffening can be described by a parameter as,

efem = W Wfem 1

(6.13)

From Equation (6.12), we already have an a priori prediction for this factor as,

e = I I 1 = kGAl2 3EI

(6.14)

We can now re-interpret the results shown in Table 6.1 for the thin beam case. Using Equations (6.13) and (6.14), we can argue a priori that the inconsistent element will produce normalized tip deflections (W fem W ) = 1 (1 + e ) . Since e>>1, we have

W fem W = N 2 5 10 5

(6.15)

for the thin beam. Table 6.2 shows how the predictions made thus compare with the results obtained from an actual finite element computation using the fieldinconsistent element. This has shown us that the consistency paradigm can be scientifically verified. Traditional procedures such as counting constraint indices, or computing the rank or condition number of the stiffness matrices could offer only a heuristic picture of how and why locking sets in. It will be instructive to note here that conventional error analysis norms in the finite element method are based on the percentage error or equivalent in some computed value as compared to the theoretically predicted value. We have seen now that the error of shear locking can be exaggerated without limit, as the structural parameter that acts as a penalty multiplier becomes indefinitely

Table 6.2 - Normalized tip deflections for the thin beam (Case 2) computed from fem model and predicted from error model (Equation (6.15)).

N 1 2 4 8 16

Computed 0.200 0.800 0.320 0.128 0.512

(fem) -4 10 -4 10 -3 10 -3 10 -3 10

Predicted -4 0.200 10 -4 0.800 10 -3 0.320 10 -3 0.128 10 -3 0.512 10

62

large. The percentage error norms therefore saturate quickly to a value approaching 100% and do not sensibly reflect the relationship between error and the structural parameter even on a logarithmic plot. A new error norm called the additional stiffening parameter, e can be introduced to recognize the manner in which the errors of locking kind can be blown out of proportion by a large variation in the structural parameter. Essentially, this takes into account, the fact that the spurious constraints give rise to a spurious energy term and consequently alters the rigidity of the system being modeled. In many other examples (e.g. Mindlin plates, curved beams etc.) it was seen that the rigidity, I, of the field consistent system and the rigidity, I, of the inconsistent 2 system, were related to the structural parameters in the form, I/I = (l/t) where l is an element dimension and t is the element thickness. Thus, if w is the deflection of a reference point as predicted by an analytical solution to the theoretical description of the problem and wfem is the fem deflection predicted by a field inconsistent finite element model, we would expect the relationship described by Equation 6.14. A logarithmic plot of the new error norm against the parameter (l/t) will show a quadratic relationship that will continue indefinitely as (l/t) is increased. This was found to be true of the many constrained media problems. By way of illustration of the distinction made

Fig. 6.2 Variation of error norms e, E with structural parameter kGL /Et cantilever beam under tip shear force.

for a

63

by this definition, we shall anticipate again, the results above. If we represent the conventional error norm in the form E = (W W fem ) W , and plot both E and the new error norm e from the results for the same problem using 4 FI 2 elements against the penalty multiplier (l/t) on a logarithmic scale, the dependence is as shown in Fig. 6.2. It can be seen that E saturates quickly to a value approaching 100% and cannot show meaningfully how the error propagates as the penalty multiplier increases indefinitely. On the other hand, e captures this relationship, very accurately.

6.1.6 Shear Force Oscillations

A feature of inconsistently modeled constrained media problems is the presence of spurious violently oscillating strains and stresses. It was not understood for a very long time that in many cases, stress oscillations originated from the inconsistent constraints. For a cantilever beam under constant bending moment modeled using linear Timoshenko beam elements, the shear force (stresses) displays a saw-tooth pattern (we shall see later that a plane stress model using 4-node elements will also give an identical pattern on the neutral bending surface). We can arrive at a prediction for these oscillations by applying the functional re-constitution technique. If V is the shear force predicted by a field-consistent shear strain field (we shall see soon how the field-consistent element can be designed) and V the shear force obtained from the original shear strain field, we can write from Equation (6.5b),

V = kGA (b0 a1 l ) V = V + kGA b1 (x l )

(6.16a) (6.16b)

We see that V has a linear term that relates directly to the constant that appeared in the spurious constraint, Equation (6.7b). We shall see below from Equation (6.17) that b1 will not be zero, in fact it is a measure of the bending moment at the centroid of the element. Thus, in a field-inconsistent formulation, this constant will activate a violent linear shear force variation when the shear forces are evaluated directly from the shear strain field given in Equation (6.5b). The oscillation is self-equilibrating and does not contribute to the force equilibrium over the element. However, it contributes a finite energy in Equation (6.9) and in the modeling of very slender beams, this spurious energy is so large as to completely dominate the behavior of the beam and cause a locking effect. Figure 6.3 shows the shear force oscillations in a typical problem - a straight cantilever beam with a concentrated moment at the tip. One to ten equal length field-inconsistent elements were used and shear forces were computed at the nodes of each element. In each case, only the variation within the element at the fixed end is shown, as the pattern repeats itself in a saw-tooth manner over all other elements. At element mid-nodes, the correct shear force i.e. V=0 is reproduced. Over the length of the element, the oscillations are seen to be linear functions corresponding to the kGA b1 (x/l) term. Also indicated by the solid lines, is the prediction made by the functional re-constitution exercise. We shall explore this now.

64

Fig. 6.3 Shear force oscillations in element nearest the root, for N element models of a cantilever of length L = 60.

Consider a straight cantilever beam with a tip shear force Q at the free end. This should produce a linearly varying bending moment M and a constant shear force Q in the beam. An element of length 2l at any station on the beam will now respond in the following manner. Since, a linear element is used, only the average of the linearly varying bending moment is expected in each finite element. If the element is field-consistent, the constant b1 can be associated after accounting for discretization, to relate to the constant bending moment M0 at the element centroid as,

M 0 = EI b1 l b1 = M 0l EI

or (6.17)

In a field-inconsistent problem, due to shear locking, it is necessary to consider the modified flexural rigidity I (see Equation 6.17) that modifies b1 to b1 , that is,

b1 = M 0l EI = {M 0l EI (1 + e )} = b1 (1 + e )
where e = kGAl2 3EI .

(6.18)

Thus, in a field-inconsistent formulation, the constant b1 gets stiffened by the factor e; the constant bending moment M0 is also underestimated by the

65

same factor. Also, for a very thin beam where e>>1, the centroidal moment M0 2 predicted by a field-consistent element diminishes in a t rate for a beam of rectangular cross-section. These observations have been confirmed through digital computation. The field-consistent element will respond with V = V0 = Q over the entire element length 2l. The field-inconsistent shear force V from Equations (6.16) and (6.18) can be written for a very thin beam (e>>1) as,

V = Q + (3M 0 l ) (x l )

(6.19)

These are the violent shear force linear oscillations within each element, which originate directly from the field-inconsistency in the shear strain definition. These oscillations are also seen if field-consistency had been achieved in the element by using reduced integration for the shear strain energy. Unless the shear force is sampled at the element centroid (i.e. Gaussian point, x/l=0), these disturbances will be much more violent than in the exactly integrated version.

6.1.7 The consistent formulation of the linear element

We can see that reduced integration ensures that the inconsistent constraint does not appear and so is effective in producing a consistent element, at least in this instance. We must now satisfy ourselves that such a modification did not violate any variational theorem. The field-consistent element, as we now shall call an element version free of spurious (i.e. inconsistent) constraints, can and has been formulated in various other ways as well. The `trick' is to evaluate the shear strain energy, in this instance, in such a way that only the consistent term will contribute to the shear strain energy. Techniques like addition of bubble modes, hybrid methods etc. can produce the same results, but in all cases, the need for consistency of the constrained strain field must be absolutely met. We explain now why the use of a trick like the reduced integration technique, or the use of assumed strain methods allows the locking problem to be overcome. It is obvious that it is not possible to reconcile this within the ambit of the minimum total potential principle only, which had been the starting point of the conventional formulation. We saw in Chapter 2, an excellent example of a situation where it was necessary to proceed to a more general theorem (one of the so-called mixed theorems) to explain why the finite element method computed strain and stress fields in a `best-fit' sense. We can now see that in the case of constrained media problems, the mixed theorem such as the Hu-Washizu or Hellinger-Reissner theorem can play a crucial role in proving that by modifying the minimum total potential based finite element formulation by using an assumed strain field to replace the kinematically derived constrained field, no energy, or work principle or variational norms have been violated.

66

To eliminate problems such as locking, we look for a consistent constrained strain field to replace the inconsistent kinematically derived strain field in the minimum total potential principle. By closely examining the strain gradient operators, it is possible to identify the order up to which the consistent strain field must be interpolated. In this case, for the linear displacement interpolations, Equations (6.5b), (6.7a) and (6.7b) tell us that the consistent interpolation should be a constant. At this point we shall still not presume what this constant should be, although past experience suggests it is the same constant term seen in Equation (6.7a). Instead, we bring in the Hellinger-Reissner theorem in the following form to see the identity of the consistent strain field clearly. For now, it is sufficient to note that the Hellinger-Reissner theorem is a restricted case of the Hu-Washizu theorem. In this theorem, the functional is stated in the following form,
T T 1 2 EI + EI 1 2 kGA

+ kGA

dx

(6.20)

where

and

are the new strain variables introduced into this multi-field

principle. Since we have difficulty only with the kinematically derived we can have = and recommend the use of a which is of consistent order to replace

. A variation of the functional in Equation (6.20) with respect to the as yet undetermined coefficients in the interpolation for yields

) dx = 0

(6.21)

This orthogonality condition now offers a means to constitute the coefficients of the consistent strain field from the already known coefficients of the kinematically derived strain field. Thus, for given by Equation (6.5b), it is possible to show that = (b0 1 l ) . In this simple instance, the same result is obtained by sampling the shear strain at the centroid, or by the use of onepoint Gaussian integration. What is important is that, deriving the consistent strain-field using this orthogonality relation and then using this to compute the corresponding strain energy will yield a field-consistent element which does not violate any of the variational norms, i.e. an exact equivalence to the mixed element exists without having to go through the additional operations in a mixed or hybrid finite element formulation, at least in this simple instance. We say that the variational correctness of the procedure is assured. The substitute strain interpolations derived thus can therefore be easily coded in the form of strain function subroutines and used directly in the displacement type element stiffness derivations.

6.1.8 Some concluding remarks on the linear beam element

So far we have seen the linear beam element as an example to demonstrate the principles involved in the finite element modeling of a constrained media problem. We have been able to demonstrate that a conceptual framework that includes a condition that specifies that the strain fields which are to be constrained must satisfy a consistency criterion is able to provide a complete scientific basis for the locking problems encountered in conventional displacement type modeling. We have also shown that a correctness criterion (which links the assumed strain variation of the displacement type formulation

67

to the mixed variational theorems) allows us to determine the consistent strain field interpolation in a unique and mathematically satisfying manner. It will be useful now to see how these concepts work if a quadratic beam element is to be designed. This is a valuable exercise as later, the quadratic beam element shall be used to examine problems such as encountered in curved beam and shell elements and in quadrilateral plate elements due to non-uniform mapping.

6.2 The quadratic Timoshenko beam element

We shall now very quickly see how the field-consistency rules explain the behavior of a higher order element. We saw in Chapter 5 that the conventional formulation with lowest order interpolation functions led to spurious constraints and a non-singular assembled stiffness matrix, which result in locking. In a higher order formulation, the matrix was singular but the spurious constraints resulted in a system that had a higher rank than was felt to be desirable. This resulted in sub-optimal performance of the approximation. We can now use the quadratic beam element to demonstrate that this is true in finite element approximations as well.

6.2.1 The conventional formulation

Consider a quadratic beam element designed according to conventional principles, i.e. exact integration of all energy terms arising from a minimum total potential principle. As the beam becomes very thin, the element does not lock; in fact it produces reasonably meaningful results. Fig. 6.4 shows a typical comparison between the linear and quadratic beam elements in its application to a simple problem. A uniform cantilever beam of length 1.0 m, width 0.01 m and 10 depth 0.01 m has a vertical tip load of 100 N applied at the tip. For E=10 2 N/m and =0.3, the engineering theory of beams predicts a tip deflection of w=4.0 m. We shall consider three finite element idealizations of this problem with the linear 2-node field-consistent element considered earlier in this section (2C, on the Figure), the quadratic 3-node field-inconsistent element being discussed now (3I, on the Figure) and the quadratic 3-node fieldconsistent element which we shall derive later (3C). It is seen that for this simple problem, the 3C element produces exact results, as it is able to simulate the constant shear and linear bending moment variation along the beam length. The 3I and 2C elements show identical convergence trends and behave as if they are exactly alike. The curious aspects that call for further investigation are: the quadratic element (3I) behaves in exactly the same way as the fieldconsistent linear element (2C), giving exactly the same accuracy for the same number of elements although the system assembled from the former had nearly twice as many nodes. It also produced moment predictions, which were identical, i.e., the quadratic beam element, instead of being able to produce linearaccurate bending moments could now yield only a constant bending moment within each element, as in the field-consistent linear element. Further, there were now quadratic oscillations in the shear force predictions for such an element. Note now that these curious features cannot be explained from the old arguments, which linked locking to the non-singularity or the large rank or the spectral condition number of the stiffness matrix. We shall now proceed to explain these features using the field-consistency paradigm.

68

Fig. 6.4 A uniform cantilever beam with tip shear force convergence trends of linear and quadratic elements. If quadratic isoparametric functions are used for the field-variables w and in the following manner

w = a0 + a1 (x l ) + a2 (x l )2

= b0 + b1 (x l ) + b2 (x l )2
the shear strain interpolation will be,

= (b0 + b2 3 a1 l ) + (b1 2a2 l ) b2 3 1 3 2

(6.22)

Again, we emphasize the usefulness of expanding the strain field in terms of the Legendre polynomials. When the strain energies are integrated, because of the orthogonal nature of the Legendre polynomials the discretized energy expression becomes the sum of the squares of the coefficients multiplying the Legendre polynomials. Indeed, the strain energy due to transverse shear strain is,

U s = 1 2 (kGA ) (2l )

{(b

2 + b2 3 a1 l )2 + 1 3 (b1 2a2 l )2 + 4b2 45

(6.23)

Therefore, when we introduce the penalty limit condition that for a thin beam the shear strain energies must vanish, we can argue that the coefficients of the strain field expanded in terms of the Legendre polynomials must vanish separately. In this case, three constraints emerge:

(b0 + b2

3 a1 l ) 0

(6.24a)

69

(b1 2a2 l ) 0
b2 0

(6.24b) (6.24c)

Equations (6.24a) and (6.24b) represent constraints having contributions from the field interpolations for both w and . They can therefore reproduce, in a consistent manner, true constraints that reflect a physically meaningful imposition of the thin beam Kirchhoff constraint. This is therefore the fieldconsistent part of the shear strain interpolation. Equation (6.24c) however contains a constant only from the interpolation for . This constraint, when enforced, is an unnecessary restriction on the freedom of the interpolation for , constraining it in fact to behave only as a

linear interpolation as the constraint implies that ,xx0 in a discretized sense over each beam element region. The spurious energy implied by such a constraint does not contribute directly to the discretized bending energy, unlike the linear beam element seen earlier. Therefore, field-inconsistency in this element would not cause the element to lock. However, it will diminish the rate of convergence of the element and would induce disturbances in the form of violent quadratic oscillations in the shear force predictions, as we shall see in the next section.

6.2.2 Functional reconstitution

We can use the functional re-constitution technique to see how the inconsistent terms in the shear strain interpolation alter the description of the physics of the original problem (we shall skip most of the details, as the material is available in greater detail in Ref. 6.1). The b2 term that appears in the bending energy also makes an appearance in the shear strain energy, reflecting its origin through the spurious constraint. We can argue that this accounts for the poor behavior of the field-inconsistent quadratic beam element (the 3I of Fig. 6.4). Ref. 6.1 derives the effect more precisely, demonstrating that the following features can be fully accounted for: i) the displacement predictions of the 3I element are identical to that made by the 2C element on an element by element basis although it has an additional midnode and has been provided with the more accurate quadratic interpolation functions. ii) the 3I element can predict only a constant moment within each element, exactly as the 2C element does. iii) there are quadratic oscillations in the shear force field within each element. We have already discussed earlier that the 3I element (the fieldinconsistent 3-noded quadratic) converges in exactly the same manner as the 2C element (the field-consistent linear). This has been explained by showing using the functional re-constitution technique, that the b2 term, which describes the linear variation in the bending strain and bending moment interpolation, is "locked" to a vanishingly small value. The 3I element then effectively behaves as a 2C element in being able to simulate only a constant bending-moment in each region of a beam, which it replaces.

70

6.2.3 The consistent formulation of the quadratic element

As in the linear element earlier, the field-consistent element (3C) can be formulated in various ways. Reduced integration of the shear strain energy using a 2-point Gauss-Legendre formula was the most popular method of deriving the element so far. Let us now derive this element using the `assumed' strain approach. We use the inverted commas to denote that the strain is not assumed in an arbitrary fashion but is actually uniquely determined by the consistency and the variational correctness requirements. The re-constitution of the field is to be done in a variationally correct way, i.e. we are required to replace in Equation (6.22) which had been derived from the kinematically admissible displacement field interpolations using the strain-displacement operators with an `assumed' strain field which contains terms only upto and including the linear Legendre polynomial in keeping with the consistency requirement. Let us write this in the form

= c0 + c1
The orthogonality condition in Equation (6.21) dictates how

(6.25) should replace

over the length of the element. This determines how c0 and c1 should be constituted from b0, b1 and b2. Fortunately, the orthogonal nature of the Legendre polynomials allows this to be done for this example in a very trivial fashion. The quadratic Legendre polynomial and its coefficient are simply truncated and c0=b0 and c1=b1 represent the variationally correct fieldconsistent `assumed' strain field. The use of such an interpolation subsequently in the integration of the shear strain energy is identical to the use of reduced integration or the use of a hybrid assumed strain approach. In a hybrid assumed strain approach, such a consistent re-constitution is automatically implied in the choice of assumed strain functions and the operations leading to the derivation of the flexibility matrix and its inversion leading to the final stiffness matrix.

6.3 The Mindlin plate elements

A very large part of structural analysis deals with the estimation of stresses and displacements in thin flexible structures under lateral loads using what is called plate theory. Thus, plate elements are the most commonly used elements in general purpose structural analysis. At first, most General Purpose Packages (GPPs) for structural analysis used plate elements based on what are called the 1 C theories. Such theories had difficulties and limitations and a1so attention 0 turned to what are called the C theories. The Mindlin plate theory [6.2] is now the most commonly used basis for the development of plate elements, especially as they can cover applications to moderately thick and laminated plate and shell constructions. It has been estimated that in large scale production runs using finite element packages, the simple four-node quadrilateral plate element (the QUAD4 element) may account for as much as 80% of all usage. It is therefore important to understand that the evolution of the current generation of QUAD4 elements from those of yester-year, over a span of nearly three decades was made difficult by the presence of shear locking. We shall now see how this takes place.

71

The history behind the discovery of shear locking in plate elements is quite interesting. It was first recognized when an attempt was made to represent the behavior of shells using what is called the degenerate shell approach [6.3]. In this the shell behavior is modeled directly after a slight modification of the 3D equations and shell geometry and domain are represented by a 3D brick element but its degrees of freedom are condensed to three displacements and two section rotations at each node. Unlike classical plate or shell theory, the transverse shear strain and its energy is therefore accounted for in this formulation. Such an approach was therefore equivalent to a Mindlin theory formulation. These elements behaved very poorly in representing even the trivial example of a plate in bending and the errors progressed without limit, as the plates became thinner. The difficulty was attributed to shear locking. This is in fact the two-dimensional manifestation of the same problem that we encountered for the Timoshenko beam element; ironically it was noticed first in the degenerate shell element and was only later related to the problems in designing Timoshenko beam and Mindlin plate elements [6.4]. The remedy proposed at once was the reduced integration of the shear strain energy [6.5,6.6]. This was only partially successful and many issues remained unresolved. Some of these were, i) the 22 rule failed to remove shear locking in the 8-node serendipity plate element, ii) the 22 rule in the 9-node Lagrangian element removed locking but introduced zero energy modes, iii) the selective 23 and 32 rule for the transverse shear strain energies from xz and yz recommended for a 8-node element also failed to remove shear locking, iv) the same selective 23 and 32 rule when applied to a 9-noded element is optimal for a rectangular form of the element but not when the element was distorted into a general quadrilateral form, v) even after reduced integration of the shear energy terms, the degenerate shell elements performed poorly when trying to represent the bending of curved shells, due to an additional factor, identified as membrane locking [6.7], originating now from the need for consistency of the membrane strain interpolations. We shall consider the membrane-locking phenomenon in another section. We shall confine our study now to plate elements without going into the complexities of the curved shell elements. In Kirchhoff-Love thin plate theory, the deformation is completely described by the transverse displacement w of the mid-surface. In such a description, the transverse shear deformation is ignored. To account for transverse shear effects, it is necessary to introduce additional degrees of freedom. We shall now consider Mindlin's approximations, which have permitted such an improved description of plate behavior. The degenerate shell elements that we discussed briefly at the beginning of this section can be considered to correspond to a Mindlin type representation of the transverse shear effects. In Mindlin's theory [6.2], deformation is described by three quantities, the section rotations x and y (i.e. rotations of lines normal to the midsurface of the undeformed plate) and the mid-surface deflection w. The bending strains are now derived from the section rotations and do not cause any

72

difficulty when a finite element model is made. The shear strains are now computed as the difference between the section rotations and the slopes of the neutral surfaces, thus,

xz = x w,x yz = y w,y

(6.26)

The stiffness matrix of a Mindlin plate element will now have terms from the bending strain energy and the shear strain energy. It is the inconsistent representation of the latter that causes shear locking.

6.3.1 The 4-node plate element

The 4-node bi-linear element is the simplest element based on Mindlin theory that could be devised. We shall first investigate the rectangular form of the element [6.4] as it is in this configuration that the consistency requirements can be easily understood and enforced. In fact, an optimum integration rule can be found which ensures consistency if the element is rectangular. It was established in Ref. 6.4 that an exactly integrated Mindlin plate element would lock even in its rectangular form. Locking was seen to vanish for the rectangular element if the bending energy was computed with a 22 Gaussian integration rule while a reduced 1-point rule was used for the shear strain energy. This rectangular element behaved very well if the plate was thin but the results deteriorated as the plate became thicker. Also, after distortion to a quadrilateral form, locking re-appeared. A spectral analysis of the element stiffness matrix revealed a rank deficiency - there were two zero energy mechanisms in addition to the usual three rigid body modes required for such an element. It was the formation of these mechanisms that led to the deterioration of element performance if the plate was too thick or if it was very loosely constrained. It was not clear why the quadrilateral form locked even after reduced integration. We can now demonstrate from our consistency view-point why the 1-point integration of the shear strain energy is inadequate to retain all the true Kirchhoff constraints in a rectangular thin plate element. However, we shall postpone the discussion on why such a strategy cannot preserve consistency if the element was distorted to a later section. Following Ref. [6.4], the strain energy for an isotropic, linear elastic plate element according to Mindlin theory can be constituted from its bending and shear energies as,

U = UB + US

24 1

Et2
2

) { [

2 x ,x

2 + y ,y +2 x ,x y , y

+ (1 ) 2 x,y + y,x 6k (1 ) t
2

)2 ] dx dy

(6.27)

2 2 ( x w,x ) + ( y w,y ) dx dy

73

Fig. 6.5 Cartesian and natural coordinate system for a four-node rectangular plate element. where x, y are Cartesian co-ordinates (see Fig. 6.5), w is the transverse displacement, x and y are the section rotations, E is the Young's modulus, is the Poisson's ratio, k is the shear correction factor and t is the plate thickness. The factor k is introduced to compensate for the error in approximating the shear strain as a constant over the thickness direction of a Mindlin plate. Let us now examine the field-consistency requirements for one of the shear strains, xz, in the Cartesian system. The admissible displacement field interpolations required for a 4-node element can be written in terms of the Cartesian co-ordinates itself as,

w = a0 + a1x + a2 y + a3 xy = b0 + b1x + b2 y + b3 xy
The shear strain functions is, field derived from these kinematically

(6.28a) (6.28b) admissible shape

xz = (b0 a1 ) + (b2 a3 )y + b1x + b3 xy

(6.29)

As the plate thickness is reduced to zero, the shear strains must vanish. The discretized constraints that are seen, to be enforced as xz 0 in Equation (6.29) are, b0 a1 0 (6.30a)

b2 a3 0 b1 0 b3 0

(6.30b) (6.30c) (6.30d)

74

The constraints shown in Equations (6.30a) and (6.30b) are physically meaningful and represent the Kirchhoff condition in a discretized form. Constraints (6.30c) and (6.30d) are the cause for concern here - these are the spurious or `inconsistent' constraints which lead to shear locking. Thus, in a rectangular element, the requirement for consistency of the interpolations for the shear strains in the Cartesian co-ordinate system is easily recognized as the polynomials use only Cartesian co-ordinates. Let us now try to derive the optimal element and also understand why the simple 1-point strategy of Ref. 6.4 led to zero energy mechanisms.

It is clear from Equations (6.29) and (6.30) that the terms b1x and b3xy are the inconsistent terms which will contribute to locking in the form of spurious constraints. Let us now look for optimal integration strategies for removing shear locking without introducing any zero energy mechanisms. We shall consider first, the part of the shear strain energy contributed by xz. We must integrate exactly, terms such as (b0-a1), (b2-a3)y, b1x, and b3xy. We now identify terms such as (b0-a1), (b2-a3), b1, and b3 as being equivalent to the quantities (x-w,x)0, (x-w,x),y0, (x,x)0, and (x,xy)0 where the subscript 0 denotes the values at the centroid of the element (for simplicity, we let the centroid of the element lie at the origin of the Cartesian co-ordinate system). An exact integration, that is a 22 Gaussian integration of the shear strain energy leads to

xzdx dy = 4lh [( x w,x )0 + h

2 2

2 3 ( x w,x )2 , y0 + l 2 3 ( x ,x )0 + h 2l 2 9 x ,xy

2 )0 ]

(6.31)

In the penalty limit of a thin plate, these four quantities act as constraints. The first two reproduce the true Kirchhoff constraints and the remaining two act as spurious constraints that cause shear locking by enforcing x,x0 and x,xy0 in the element. If a 12 Gaussian integration is used, we have,

xzdx dy = 4lh [( x w,x )0 + h

2 2

3 ( x w,x )2 , y0

(6.32)

Thus, only the true constraints are retained and all spurious constraints are removed. This strategy can also be seen to be variationally correct in this case; we shall see later that in a quadrilateral case, it is not possible to ensure variational correctness exactly. By a very similar argument, we can show that the part of the shear strain energy from yz will require a 21 Gaussian integration rule. This element would be the optimal rectangular bi-linear Mindlin plate element. Let us now look at the 1-point integration strategy used in Ref. 6.4. This will give shear energy terms such as,

xz dx dy = 4lh [( x w,x )0 ]
2 2

(6.33)

75

We have now only one true constraint each for the shear energy from xz and yz respectively while the other Kirchhoff constraints ( x w,x ), y0 0 and are lost. This introduces two zero energy modes and accounts for the consequent deterioration in performance of the element when the plates are thick or are very loosely constrained, as shown in Ref. 6.4. We have seen now that it is a very simple procedure to re-constitute field-consistent assumed strain fields from the kinematically derived fields such as shown in Equation (6.29) so that they are also variationally correct. This is not so simple in a general quadrilateral where the complication arising from the isoparametric mapping from a natural co-ordinate system to a Cartesian system makes it very difficult to see the consistent form clearly. We shall see the difficulties associated with this form in a later section.

( y w,y ),x0 0

6.3.2 The quadratic 8-node and 9-node plate elements

The 4-node plate element described above is based on bi-linear functions. It would seem that an higher order element based on quadratic functions would be far more accurate. There are now two possibilities, an 8-node element based on what are called the serendipity functions and a 9-node element based on the Lagrangian bi-quadratic functions. There has been a protracted debate on which version is more useful, both versions having fiercely committed protagonists. By now, it is well known that the 9-node element in its rectangular form is free of shear locking even with exact integration of shear energy terms and that its performance is vastly improved when its shear strain energies are integrated in a selective sense (23 and 32 rules for xz and yz terms respectively). It is in fact analogous to the quadratic Timoshenko beam element, the fieldinconsistencies not being severe enough to cause locking. This is however not true for the 8-node element which was derived from the Ahmad shell element [6.3] and which actually pre-dates the 4-node Mindlin element. An exact integration of bending and shear strain energies resulted in an element that locked for most practical boundary suppressions even in its rectangular form. Many ad-hoc techniques e.g. the reduced and selective integration techniques, hybrid and mixed methods, etc. failed or succeeded only partially. It was therefore regarded for some time as an unreliable element as no quadrature rule seemed to be able to eliminate locking entirely without introducing other deficiencies. It seems possible to attribute this noticeable difference in the performance of the 8- and 9-node elements to the missing central node in the former. This makes it more difficult to restore consistency in a simple manner.

6.3.3 Stress recovery from Mindlin plate elements

The most important information a structural analyst looks for in a typical finite element static analysis is the state of stress in the structure. It is therefore very important for one to know points of optimal stresses in the Mindlin plate elements. It is known that the stress recovery at nodes from displacement elements is unreliable, as the nodes are usually the points where the strains and stresses are least accurate. It is possible however to determine points of optimal stress recovery using an interpretation of the displacement method as a procedure that obtains strains over the finite element domain in a least-squares accurate sense. In Chapter 2, we saw a basis for this interpretation. We can apply this rule to determine points of accurate stress

76

recovery in the Mindlin plate elements. For a field-consistent rectangular 4node element, the points are very easy to determine [6.8] (note that in a fieldinconsistent 4-node element, there will be violent linear oscillations in the shear stress resultants corresponding to the inconsistent terms). Thus, Ref. 6.8 shows that bending moments and shear stress resultants Qxz and Qyz are accurate at the centroid and at the 12 and 21 Gauss points in a rectangular element for isotropic or orthotropic material. It is coincidental, and therefore fortuitous, that the shear stress resultants are most accurate at the same points at which they must be sampled in a selective integration strategy to remove the fieldinconsistencies! For anisotropic cases, it is safest to sample all stress resultants (bending and shear) at the centroid. Such rules can be extended directly to the 9-node rectangular element. The bending moments are now accurate at the 22 Gauss points and the shear stress resultants in an isotropic or orthotropic problem are optimal at the same 23 and 32 Gauss points which were used to remove the inconsistencies from the strain definitions. However, accurate recovery of stresses from the 8-node element is still a very challenging task because of the difficulty in formulating a robust element. The most efficient elements known today are variationally incorrect even after being made field-consistent and need special filtering techniques before the shear stress resultants can be reliably sampled. So far, discussion on stress sampling has been confined to rectangular elements. When the elements are distorted, it is no simple matter to determine the optimal points for stress recovery - the stress analyst must then exercise care in applying these rules to seek reliable points for recovering stresses.

6.4 Concluding remarks

We can conclude this section on shear locking by noting that the available understanding was unable to resolve the phenomena convincingly. The proposed improvement, which was the consistency paradigm, together with the functional re-constitution procedure, allowed us to derive an error estimate for a case under locking and we could show through numerical (digital) experiments that these estimates were accurate. In this way we are convinced that a theory with the consistency paradigm is more successful from the falsifiability point of view than one without.

6.5 References
6.1 G. Prathap and C. R. Babu, Field-consistent strain interpolations for the quadratic shear flexible beam element, Int. J. Num. Meth. Engng. 23, 19731984, 1986. 6.2 R. D. Mindlin, Influence of rotary inertia and shear on flexural motion of elastic plates, J. Appl. Mech. 18, 31-38, 1951. 6.3 S. Ahmad, B. M. Irons and O. C. Zienkiewicz, Analysis of thick and thin shell structures by curved finite elements, Int. J. Num. Meth. Engng. 2, 419-451, 1970. 6.4 T. J. R. Hughes, R. L. Taylor and W. Kanoknukulchal, A simple and efficient finite element for plate bending, Int. J. Num. Meth. Engng. 411, 15291543, 1977. 6.5 S. F. Pawsey and R. W. Clough, Improved numerical integration of thick shell finite elements, Int. J. Num. Meth. Engng. 43, 575-586, 1971.

77

6.6 O. C. Zienkiewicz, R. L. Taylor and J. M. Too, Reduced integration technique in general analysis of plates and shells, Int. J. Num. Meth. Engng. 43, 275290, 1971. 6.7 H. Stolarski and T. Belytschko, Membrane locking and reduced ingression for curved elements, J. Appl. Mech. 49, 172-178, 1982. 6.8 G. Prathap and C. R. Babu, Accurate force evaluation with a simple bi-linear plate bending element, Comp. Struct. 25, 259-270, 1987.

78

Chapter 7 Membrane locking, parasitic shear and incompressible locking

7.1 Introduction
The shear locking phenomenon was the first of the over-stiffening behavior that was classified as a locking problem. Other such pathologies were noticed in the behavior of curved beams and shells, in plane stress modeling and in modeling of three dimensional elastic behavior at the incompressible limit (0.5). The field-consistency paradigm now allows all these phenomena to be traced to inconsistent representations of the constrained strain fields. A unifying pattern is therefore introduced to the understanding of the locking problems in constrained media elasticity - whether it is shear locking in a straight beam or flat plate element, membrane locking in a curved beam or shell element, parasitic shear in 2D plane stress and 3D elasticity or incompressible locking in 2D plane strain and 3D elasticity as the Poisson's ratio 0.5. This chapter will now briefly summarize this interpretation. Unlike the chapter on shear locking, a detailed analysis is omitted as it would be beyond the scope of the present book and readers should refer to the primary published literature.

7.2 Membrane locking

Earlier, we argued that the most part of structural analysis deals with the behavior of thin flexible structures. One popular and efficient form of construction of such thin flexible structures is the shell - enclosing space using one or more curved surfaces. A shell is therefore the curved form of a plate and its structural action is a combination of stretching and bending. It is possible to perform a finite element analysis of a shell by using what is called a facet representation - i.e. the shell surface is replaced with flat triangular and/or quadrilateral plate elements in which a membrane stiffness (membrane element) is superposed on a bending stiffness (plate bending element). Such a model is understandably inaccurate in that with very coarse meshes, they do not capture the bending-stretching coupling of thin shell behavior. Hence, the motivation for designing elements with mid-surface curvature taken into account and which can capture the membrane-bending coupling correctly. There are two ways in which this could be done. One is to use elements based on specific shell theories (e.g. the Donnell, Flugge, Sanders, Vlasov theories, etc.). There are considerable controversies regarding the relative merits of these theories as each has been obtained by carrying out approximations to different degrees when the 3-dimensional field equations are reduced to the particular class of shell equations. The second approach is called the degenerate shell approach three dimensional solid elements can be reduced (degenerated) into shell elements having only mid-surface nodal variables - these are no longer dependent on the various forms of shell theories proposed and should be simple to use. They are in fact equivalent to a Mindlin type curved shell element. However, accurate and robust difficult to design - a problem decades. One key issue behind this locking - the very poor behavior of curved shell elements have been extremely challenging researchers for nearly three difficulty is the phenomenon of membrane curved elements, as they become thin. This

79

Fig. 7.1 The classical thin curved beam element. was neither recognized nor understood for a very long time. It was believed at first that this was because curved elements could not satisfy the strain-free rigid body motion condition because of their inherent curved geometry. In this section, we shall apply the field-consistency paradigm to the problem to recognize that a consistent representation of membrane strains is desired to avoid the phenomenon called membrane locking and that it is the need for consistency rather than the requirement for strain-free rigid body motion which is the key factor. This is most easily identified by working with the onedimensional curved beam (or arch) elements. In Chapter 6, we studied structural problems where both bending and shear deformation was present. We saw that simple finite element models for the shear deformable bending action of a beam or plate behaved very poorly in thin beam or plate regimes where the shear strains become vanishingly small when compared with the bending strains. This phenomenon is known as shear locking. We also saw how the field-consistency paradigm was needed to explain this behavior in a scientifically satisfying way. We shall now see that a similar problem appears in curved finite elements in which both flexural and membrane deformation take place. Such elements perform very poorly in cases where inextensional bending of the curved structure was predominant, i.e. the membrane strains become vanishingly small when compared to the bending strains. This phenomenon has now come to be called membrane locking. In this section, we shall examine in turn, the classical thin curved beam element and the linear and quadratic shear flexible curved beam elements in a curvilinear co-ordinate system.

7.2.1 The classical thin curved beam element

Figure 7.1 describes the simplest curved beam element of length 2l and radius of curvature R based on classical thin beam theory. The displacement degrees of freedom required are the circumferential displacement u and the radial displacement w. The co-ordinate s follows the middle line of the curved beam. The membrane strain and the bending strain are described by the straindisplacement relations

= u,s + w R = u,s R w,ss

A C description for u and a C1 description w is admissible displacement interpolations for u and w are
0

(7.1a) (7.1b) required. Kinematically

80

Fig. 7.2 Geometry of a circular arch

u = a0 + a1 w = b0 + b1 + b2
2

(7.2a)

+ b3

(7.2b)

where = s l and a0 to b3 are the generalized degrees of freedom which can be related to the nodal degrees of freedom u, w and w,s at the two nodes. The strain field interpolations can be derived as

= a1 Rl 2b2 l 2 6b3 l 2 = (a1 l + b0 R + b2 3R ) + (b1 R + 3b3 5R ) b 2

2

) ( ) 3R (1 3 ) b

(7.3a)

5R 3 5 3

(7.3b)

If a thin and deep arch (i.e. having a large span to rise ratio so that L/t>>1 and R/H is small, see Fig. 7.2) is modeled by the curved beam element, the physical response is one known as inextensional bending such that the membrane strain tends to vanish. From Equation (7.3b) we see that the inextensibility condition leads to the following constraints:

a1 l + b0 R + b2 3R 0 b1 + 3b3 5 0 b2 0 b3 0

(7.4a) (7.4b) (7.4c) (7.4d)

Following the classification system we introduced earlier, we can observe that constraint (7.4a) has terms participating from both the u and w fields. It can therefore represent the condition = u,s + w R 0 in a physically meaningful way,

81

i.e. it is a true constraint. However, the three remaining constraint (7.4b) to (7.4d) have no participation from the u field. Let us now examine in turn what these three constraints imply for the physical problem. From the three constraints we have the conditions b1 0, b2 0 and b3 0. Each in turn implies the conditions w,s 0, w,ss 0 and w,sss 0 . These are the spurious constraints. This element has been studied extensively and the numerical evidence given there confirms that membrane locking is mainly determined by the first of these three constraints. Next, we shall briefly describe the functional re-constitution analysis that will yield semi-quantitative error estimates for membrane locking in this element. We can now use the functional re-constitution technique, (which in Chapter 6 revealed how the discretization process associated with the linear Timoshenko beam element led to large errors known as shear locking) to see how the inconsistent constraints in the membrane strain field interpolation (Equations (7.4b) to (7.4d)) modify the physics of the curved beam problem in the discretization process. The analysis is a little more difficult than for the shear locking case and would not be described in full here. For a curved beam of length 2l, moment of inertia I and area of cross-section A, the strain energy can be written as,

e = 1 2 EI T + 1 2 EA T ds

(7.5)

where and are as defined in Equation (7.1). For simplicity assume that the cross-section of the beam is rectangular so that if the depth of the beam is t, then I = At2 12 . After discretisation, and reconstitution of the functional, we can show that the inconsistent membrane energy terms disturb the bending energy terms - the principal stiffening term is physically equivalent to a spurious inplane or membrane stiffening action due to an axial load N = EAl2 3R 2 . This action is quite complex and it is not as simple to derive useful universally valid error estimates as was possible with the linear Timoshenko beam element. We therefore choose a problem of simple configuration which can reveal the manner in which the membrane locking action takes place. We consider a shallow circular arch of length L and radius of curvature R, simply supported at the ends. The loading system acting on it is assumed to induce nearly inextensional bending. We also assume that the transverse deflection w can be approximated in the form w=csin(s/L). If the entire arch is discretized by curved beam elements of length 2l each, the strain energy of the discretized arch (again assuming that =-w,ss and neglecting the consistently modeled part of the membrane energy) one can show that due to membrane locking, the bending rigidity of the arch has now been altered in the following form:

(EI ) (EI ) = { + (4 1

2 (Ll Rt )2

(7.6)

We can therefore expect membrane locking to depend on the structural parameter 2 of the type (Ll/Rt) .

82

Fig. 7.3 The Mindlin curved beam element. Numerical experiments with various order integration rules applied to evaluate the various components of the stiffness matrix indicated that the inconsistent elements locked exactly as predicted by the pattern describe above. The consistent element was entirely free of locking. This classical curved beam element (the cubic-linear element) was discredited for a very long time because it was not possible to understand why the element should perform so poorly. We have now found that the explanation for its very poor behavior could be attributed entirely to the inconsistency in the membrane strain field interpolations. The belief that the lack of the strain free rigid body motion in the conventional formulation was the cause of these errors cannot be substantiated at all. We have also seen how we could restore the element into a very accurate one (having the accuracy of the straight beam models) by designing a new membrane strain field which met the consistency requirements.

7.2.2 The Mindlin curved beam elements

The exercise above showed that the locking behavior of curved beams could be traced to the inconsistent definition of the membrane strain field. The curved beam element used to demonstrate this was based on classical thin beam theory and used a curvilinear co-ordinate description. Of interest to us now is to understand how the membrane locking phenomenon can take place in a general curved shell element. Most general shell elements which have now been accepted into the libraries of general purpose finite element packages are based on what is called the degenerate shell theory which is equivalent to a shear deformable theory and use quadratic interpolations (i.e. three nodes on each edge to capture the curved geometry in a Cartesian system). To place the context of membrane locking in general shell elements correctly, it will be useful to investigate the behavior of the Mindlin curved beam elements. These are now described by the circumferential (tangential) displacement u, the radial (normal) displacement w and the rotation of a normal to the midline . As in the Timoshenko beam element, this rotation differs from the rotation of the mid-line w,s by an amount = w,s which describes the shear strain at the section. Fig. 7.3 shows how the Mindlin curved beam element is specified. The strain energy in a Mindlin curved beam is U=UM+UB+US where the respective contributions to U arise from the membrane strain , the bending

83

strain and the transverse shear strain . For an element of length 2l and radius R,

U M = (1 2 ) (EA ) 2 ds where = u,s + w R U B = (1 2 ) (EI ) 2 ds where = u,s R ,s U S = (1 2 ) (kGA ) 2 ds where = w ,s

where EA, EI and kGA represent the respective rigidities.

(7.7a) (7.7b) (7.7c)

There are no difficulties with the bending strain. However, depending on the regime, both shear strains and membrane strains can be constrained to vanish leading to shear and membrane locking. The arguments concerning shear locking are identical to those expressed in Chapter 6, i.e. while the inconsistency in the shear strain field is severe enough to cause locking in the former (a significant stiffening which caused errors that propagated indefinitely as the structural parameter kGl2 Et2 became very large), it is much milder in the latter and is responsible only for a poorer rate of convergence (the performance of the element having reduced to that of the linear field-consistent element). However with membrane locking, the inconsistencies at the quadratic level can cause locking - therefore both the linear and quadratic elements must be examined in turn to get the complete picture. For the linear Mindlin curved beam element, as a C displacement type formulation is sufficient, we can choose linear interpolations of the following form,
0

u = a0 + a1 w = b0 + b1 = c0 + c1
where = s l is the natural co-ordinate along the arch mid-line. interpolations lead to the following membrane strain field interpolation:

(7.8a) (7.8b) (7.8c) These

= (a1 l + b0 R ) + (b1 R )

(7.9)

There is now one single spurious constraint in the limit of inextensional bending of an arch, i.e. b10, which implies a constraint of the form w,s0 at the centroid of the element. The membrane locking action is similar to but simpler than that seen for the classical thin curved beam element above. In a thin beam, in the Kirchhoff limit, w,s, which implies that membrane locking induces the section rotations to lock. It can therefore be seen that a constraint of this form leads to an in-plane stiffening action of the form corresponding to the introduction of a spurious energy

(1 2 ) (EAl2

3R 2

) w ,

2 s

ds

(7.10)

84

The membrane locking action is therefore identical to that predicted for the classical thin curved beam. A variationally correct field-consistent element can easily be derived using the Hu-Washizu theorem. This is a very accurate element. The quadratic Mindlin curved beam element can also be interpreted in the light of our studies of the classical curved beam earlier. Spurious constraints at both the linear (i.e. w,s0) and quadratic (i.e w,ss0) levels lead to rate respectively. We errors that propagate at a (Ll Rt )2 rate and l 2 Rt expect now that in a quadratic Mindlin curved beam element, the latter constraint is the one that is present and that these errors are large enough to make the membrane locking effect noticeable. This also explains why the degenerate shell elements, which are based on quadratic interpolations, suffer from membrane locking.
2

A quadratic (three-node) element having nodes at s=-l, 0 and l (see Fig. 7.3) would have a displacement field specified as,

u = a0 + a1 + a2 2 w = b0 + b1 + b2
2

(7.11a) (7.11b) (7.11c)

= c0 + c1 + c2 2
The membrane strain field can be expressed as:

= (a1 l + b0 R + b2 3R ) + (2a2 l + b1 R ) b2 3R 1 3 2

(7.12)

The spurious constraint is now b20 a discretized equivalent of imposing the constraint w,ss0 at the centroid of each element. To understand that this can cause locking, we note that in a thin curved beam, ,s w,ss and therefore a spurious constraint of this form acts to induce a spurious bending energy. A functional re-constitution approach will show that the bending energy for a thin curved beam modifies to,
2 2 U B = (1 2 ) (EI ) 1 + (4 15 ) l 2 Rt w,ss ds

(7.13)

This can be interpreted as spurious additional stiffening of the actual bending rigidity of the structure by the (4 15 ) l 2 Rt 2 term. It is a factor of this form that caused the very poor performance of the exactly integrated quadratic shell elements. A variationally correct field consistent element is obtained by using the Hu-Washizu theorem or very simply by using a two-point integration rule for the shear and membrane strain energy energies. This is an extremely accurate element.

i.e.

(l

It is interesting to compare the error norm for locking for this element,
2

Rt2

with the corresponding error norm for locking in the inconsistent

85

linear Mindlin curved beam element, which is L l Rt 2 . It is clear that since L>>l when an arch is discretized with several elements, the locking in the quadratic element is much less than in the linear element. One way to look at this is to examine the rate at which locking is relieved as more elements are added to the mesh, i.e. as l is reduced. In the case of the inconsistent linear 2 element, locking is relieved only at an l rate whereas in the case of the 4 quadratic element, locking is relieved at the much faster l rate. However, locking in the quadratic element is significant enough to degrade the performance of the quadratic curved beam and shell elements to a level where it is impractical to use in its inconsistent form. This is clearly seen from the fact that to reduce membrane-locking effects to acceptable limits, one must use element lengths l << Rt . Thus for an arch with R=100 and t=1, one must use about a 100 elements to achieve an accuracy that can be obtained with two to four field-consistent elements.

7.2.3 Axial force oscillations

Our experience with the linear and quadratic straight beam elements showed that inconsistency in the constrained strain-field is accompanied by spurious stress oscillations corresponding to the inconsistent constraints. We can now understand that the inconsistent constraint b10 will induce linear axial force oscillations in the linear curved beam element. Similarly, the inconsistent quadratic constraint will trigger off quadratic oscillations in the quadratic element.

7.2.4 Concluding remarks

So far, we have examined the phenomenon of membrane locking in what were called the simple curved beam elements. These were based on curvilinear geometry and the membrane locking effect was easy to identify and quantify. However, curved beam and shell elements which are routinely used in general purpose finite element packages are based on quadratic formulations in a Cartesian system. We shall call these the general curved beam and shell elements. Membrane locking is not so easy to anticipate in such a formulation although the predictions made in this section have been found to be valid for such elements.

7.3 Parasitic shear

Shear locking and membrane locking are phenomena seen in finite element modeling using structural elements based on specific theories, e.g. beam, plate and shell theories. These are theories obtained by simplifying a more general continuum description, e.g. the three dimensional theory of elasticity or the two dimensional theories of plane stress, plane strain or axisymmetric elasticity. The simplifications are made by introducing restrictive assumptions on strain and/or stress conditions like vanishing of shear strains or vanishing of transverse normal strains, etc. In the beam, plate and shell theories, this allows structural behavior to be described by mid-surface quantities. However, these restrictive assumptions impose constraints which make finite element modeling quite difficult unless the consistency aspects are taken care of. We shall now examine the behavior of continuum elements designed for modeling continuum elasticity directly, i.e. 2D elements for plane strain or stress or axisymmetric elasticity and solid elements for 3D elasticity. It appears that

86

Fig. 7.4 (a) Constant strain triangle, (b) Rectangular bilinear element, (c) Eight-node rectangle and (d) Nine-node rectangle. under constraining physical regimes, locking behavior is possible. One such problem is parasitic shear. We shall see the genesis of this phenomenon and show that it can be explained using the consistency concepts. The problems of shear locking, membrane locking and parasitic shear obviously are very similar in nature and the unifying factor is the concept of consistency. The finite element modeling of plane stress, plane strain, axisymmetric and 3D elasticity problems can be made with two-dimensional elements and threedimensional elements (sometimes called continuum elements to distinguish them from structural elements), of which the linear 3-noded triangle and bi-linear 4noded rectangle (Fig. 7.4) are the simplest known in the 2D case. In most applications, these elements are reliable and accurate for determining general two-dimensional stress distributions and improvements can be obtained by using the higher order elements, the quadratic 6-noded triangle and 8- or 9-noded quadrilateral elements (Fig. 7.4). However, in some applications, e.g. the plane stress modeling of beam flexure or in plane strain and axisymmetric applications where the Poisson's ratio approaches 0.5 (for nearly incompressible materials or for materials undergoing plastic deformation), considerable difficulties are seen. Parasitic shear occurs in the plane stress and 3D finite element modeling of beam flexure. Under pure bending loads, for which the shear stress should be zero, the elements yield large values of shear stress except at certain points (e.g. centroids in linear elements). In the linear elements, the errors progress in the same fashion as for shear locking in the linear Timoshenko beam element, i.e. the errors progress indefinitely as the element aspect ratio increases. Reduced integration of the shear strain energy or addition of bubble modes for the 4-noded

87

Fig. 7.5 Plane stress model of beam flexure. rectangular element completely eliminates this problem. The 8- and 9-noded elements do not lock but improve in performance with reduced integration. These elements, in their original (i.e. field-inconsistent) form show severe linear and quadratic shear stress oscillations. Here, we shall extend the consistency paradigm to examine this phenomenon in detail, confining attention to the rectangular elements.

7.3.1 Plane stress model of beam flexure

Figure 7.5 shows a two-dimensional plane stress description of the problem of beam flexure. The beam is of length L, depth T and thickness (in the normal to the plane direction) b. Two independent field variables are required to describe the problem, the displacements u and v (Fig. 7.5). The strain energy functional for this problem consists of energies from the normal strains, UE and shear strains, UG. For an isotropic case, where E is the Young's modulus, G is the shear modulus and is Poisson's ratio,

U = U E + UG
Eb L T u 2 ,x + v2 ,y +2 u,x v , y + Gb 2 u, y +v ,x 2 dx dy = 0 0 2 2 1-

)(

(7.14)

A one-dimensional beam theory can be obtained by simplifying this problem. For a very slender beam, the shear strains must vanish and this would then yield the classical Euler-Bernouilli beam theory. Note that although shear strains vanish, shear stresses and shear forces would remain finite. We shall see now that it is this constraint condition that leads to difficulties in plane stress finite element modeling of such problems.

7.3.2 The 4-node rectangular element

An element of length 2l and depth 2t is chosen. The rectangular form is chosen so that the issues involved can be identified clearly. There are now 8 degrees of freedom, u1-u4 and v1-v4 at the four nodes. The field-variables are now interpolated in the following fashion:

88

u = a0 + a1x + a2 y + a3 xy v = b0 + b1x + b2 y + b3 xy

(7.15a) (7.15b)

We need to examine only the shear strain field that will become constrained (i.e. become vanishingly small) while modeling the flexural action of a thin beam. This will be interpolated as,

= u,y +v,x = (a2 + b1 ) + a3 x + b3 y

and the shear strain energy within the element will be,

(7.16)

U e (G ) = -l Gb 2 a2 + b1 + a3 x + b3 y
l -t

)2 dx dy

= (Gb 2 ) (4tl )

[(a

+ b1 )2 + (a3l )2 3 + (b3t )2 3

(7.17)

Equation (7.17) suggests that in the thin beam limit, the following constraints will emerge:

a2 + b1 0 a3 0 b3 0

(7.18a) (7.18b) (7.18c)

We can see that Equation (7.18a) is a constraint that comprises constants from both u and v functions. This can enforce the true constraint of vanishing shear strains in a realistic manner. In contrast, Equations (7.18b) and (7.18c) impose constraints on terms from one field function alone in each case. These are undesirable constraints. We shall now see how these lead to the stiffening effect that is called parasitic shear. In the plane stress modeling of a slender beam, we would be using elements which are very long in the x direction, i.e. elements with l>t would be used. 2 2 One can then expect from Equation (7.17) that as l >>t , the constraint a30 will be enforced more rapidly than the constraint b30 as the beam becomes thinner. The spurious energies generated from these two terms will also be in a similar proportion. In Equation (7.17), one can interpret the shear strain field to comprise a constant value, which is `field-consistent' and two linear variations which are `field-inconsistent'. We see now that the variation along the beam length is the critical term. Therefore, let us consider the use of the shear strain along y=0 as a measure of the averaged shear strain across a vertical section. This gives

= (a2 + b1 ) + a3 x

(7.19)

Thus, along the length of the beam element, have a constant term that reflects the averaged shear strain at the centroid, and a linear oscillating term along the length, which is related to the spurious constraint in Equation (7.18b). As in the linear beam element in Section 4.1.6, this oscillation is selfequilibrating and does not contribute to the force equilibrium over the element. However, it contributes a finite energy in Equation (7.17) and in the modeling

89

of very slender beams, this spurious energy is so large as to completely dominate the model of the beam behavior and cause a locking effect. We shall now use the functional re-constitution technique to determine the magnitude of locking and the extent of the shear stress oscillations triggered off by the inconsistency in the shear strain field. We shall simplify the present example and derive the error estimates from the functional re-constitution exercise. We shall denote by a3 , the coefficient determined in a solution of the beam problem by using a field-consistent representation of the total energy in the beam, and by a3 , the value from a solution using energies based on inconsistent terms. The field-consistent and field-inconsistent terms will differ by,

a3 a3 = 1 + Gl 2 Et 2

)
stiffening

(7.20) parameter that

This factor, e L = 1 + Gl2 Et2 is the additional determines the extent of parasitic shear.

We can relate a3 to the physical parameters relevant to the problem of a beam in flexure to show that if e L >> 1 , as in a very slender beam, there are shear force oscillations,

V = V0 + (3M 0 l ) (x l )

(7.21)

The oscillations are proportional to the centroidal moment, and for e L >> 1 , the oscillations are inversely proportional to the element length, i.e. they become more severe for smaller lengths of the plane stress element. It is important to note that a field-inconsistent two node beam element based on linear interpolations for w and (see Equation (7.19)) produces an identical stress oscillation.

7.3.3 The field-consistent elements

It is well known that the plane stress elements behave reliably in general 2D applications where no constraints are imposed on the normal or shear strain fields. The need for field-consistency becomes important only when a strain field is constrained. Earlier, we saw that in the plane stress modeling of thin beam flexure, the shear strain field is constrained and that these constraints are enforced by a penalty multiplier Gl2 Et2 . The larger this term is, the more severely constrained the shear strain becomes. Of the many strategies that are available for removing the spurious constraints on the shear strain field, the use of a variationally correct re-distribution is recommended. This requires the field-consistent to be determined from using the orthogonality condition

) dx dy = 0

(7.22)

Thus instead of Equation (7.16), we can show that the shear strain-field for the 4-node rectangular plane stress element will be,

90

 = (a2 + b1 )

(7.23)

It is seen that the same result is achieved if reduced integration using a 1 Point Gaussian integration is applied to shear strain energy evaluation. In a very similar fashion, the consistent shear strain field can be derived for the quadratic elements.

7.3.4 Concluding remarks on the rectangular plane stress elements

The field-consistency principle and the functional re-constitution technique can be applied to make accurate error analyses of conventional 4-node and 8-node plane stress elements. These a priori estimates have been confirmed through numerical experiments. A field-consistent re-distribution strategy allows these elements to be free of locking and spurious stress-oscillations in modeling flexure. These options can be built into a shape function sub-routine package in a simple manner so that where field-consistency requirements are paramount, the appropriate smoothed shape functions will be used.

7.3.5 Solid element model of beam/plate bending

So far, we have dealt with the finite element modeling of problems which were simplified to one-dimensional or two-dimensional descriptions. There remains a large class of problems which need to be addressed directly as three dimensional states of stress. Solid or three dimensional elements are needed to carry out the finite element modeling of such cases. A variety of 3D solid elements exist, e.g. tetrahedral, triangular prism or hexahedral elements. In general cases of 3D stress analysis, no problems are encountered with the use of any of these elements as long as a sufficiently large number of elements are used and no constrained media limits are approached. However, under such limits the tetrahedral and triangular prism elements cannot be easily modified to avoid these difficulties. It is possible however to improve the hexahedral 8-node and 27-node elements so that they are free of locking and of other ill-effects like stress oscillations. One problem encountered is that of parasitic shear or shear locking when solid elements are used to model regions where bending action is predominant. In fact, in such models, parasitic shear and shear locking merge indistinguishably into each other. We briefly examine this below. The 8-noded brick element is based on the standard tri-linear interpolation functions and is the three dimensional equivalent of the bi-linear plane stress element. It would appear now that the stiffness matrix for the element can be derived very simply by carrying out the usual finite element operations, i.e. with strain-displacement matrix, stress-strain matrix and numerically exact integration (i.e. 222 Gaussian integration). Although this element performs very well in general 3D stress analysis, it failed to represent cases of pure bending (parasitic shear) and cases of near incompressibility (incompressible locking). The element is unable to represent the shear strains in a physically meaningful form when the shear strains are to be constrained, as in a case of pure bending. The phenomenon is the 3D analogue of the problem we noticed for

91

the bi-linear plane stress element above. The problem of designing a useful node brick element has therefore received much attention and typical attempts alleviate this have included techniques such as reduced integration, addition bubble functions and assumed strain hybrid formulations. We shall now look this from the point of view of consistency of shear strain interpolation.

8to of at

We can observe here that the shear strain energy is computed from Cartesian shear strains xy , xz and yz . In the isoparametric formulation, these are to be computed from terms involving the derivatives of u, v, w with respect to the natural co-ordinates , , and the terms from the inverse of the three dimensional Jacobian matrix. It is obvious that it will be a very difficult, if not impossible task, to assure the consistency of the interpolations for the Cartesian shear strains in a distorted hexahedral. Therefore attention is confined to a regular hexahedron so that consistency requirements can be easily examined. To facilitate understanding, we shall restrict our analysis to a rectangular prismatic element so that the (x, y, z) and (, , ) systems can be used interchangeably. We can consider the tri-linear interpolations to be expanded in the following form, i.e.

u = a1 + a2 x + a3 y + a4z + a5 xy + a6 yz + a7 xz + a8 xyz v = b1 + b2 x + b3 y + b4z + b5 xy + b6 yz + b7 xz + b8 xyz w = c1 + c2 x + c3 y + c4z + c5 xy + c6 yz + c7 xz + c8 xyz

(7.24a) (7.24b) (7.24c)

where a1 to a8 and b1 to b8 are related to the nodal degrees of freedom u1 to u8 and v1 to v8 respectively. We shall now consider a case of a brick element undergoing a pure bending response requiring a constrained strain field xy = u,y +v ,x 0 . From Equations (7.24a) and (7.24b) we have,

xy = (a3 + b2 ) + (a6 + b7 )z + a5 x + b5 y + a8 xz + b8 yz
When Equation (7.25) is constrained, the conditions appear: a50; b50; a80 and b80. following spurious

(7.25) constraint

Our arguments projecting these as the cause of parasitic shear would be proved beyond doubt if we can derive a suitable error model. These error models are now much more difficult to construct, but have been achieved using a judicious mix of a priori and a posteriori knowledge of zero strain and stress conditions to simplify the problem. These error models have been convincingly verified by numerical computations [7.1]. A field-consistent representation will be one that ensures that the interpolations for the shear fields will contain only the consistent terms. A simple way to achieve this is to retain only the consistent terms from the original interpolations. Thus, for xy we have,

xy = (a3 + b2 ) + (a6 + b7 )z

(7.26)

92

Similarly, the other shear strain fields can be modified. These strain fields should now be able to respond to bending modes in all three planes without locking.

7.4 Incompressible locking

Conventional displacement formulations of 2D plane strain and 3D elasticity fail when Poisson's ratio 0.5, i.e. as the material becomes incompressible. Such situations can arise in modeling of material such as solid rocket propellants, saturated cohesive soils, plastics, elastomers and rubber like materials and in materials that flow, e.g. incompressible fluids or in plasticity. Displacement fields lock and highly oscillatory stresses are seen - drawing an analogy with shear and membrane locking, we can describe the phenomenon as incompressible locking. In fact, a penalty term (bulk modulus as 0.5) is easily identifiable as enforcing incompressibility. Since this multiplies the volumetric strain, it is easy to argue that a lack of consistency in the definition of the volumetric strain will cause locking and will appear as oscillatory mean stresses or pressures.

7.4.1 A simple one-dimensional case - an incompressible hollow sphere

It is instructive to demonstrate how field-consistency operates in the finite element approximation of a problem in incompressible (or nearly-incompressible) elasticity by taking the simplest one-dimensional case possible. A simple problem of a pressurized elastic hollow sphere serves to illustrate this. Consider an elastic hollow sphere of outer radius b and inner radius a with an internal pressure P. If the shear modulus G is 1.0, the radial displacement field for this case is given by,

1 (1 2 ) r u= + 2 2 (1 + ) 4r
where = Pba3 G b 3 a3 ratio.

(7.27)

[(

)]

and r=R/b non-dimensional radius and is the Poisson's

A finite element approximation approaches the problem from the minimum total potential principle. Thus, a radial displacement u leads to the following strains:

r = u,r
and the volumetric strain is

= u r

= u r

= r + + = u,r + 2u r
The elastic energy stored in the sphere (when G=1) is
1

U =

2 r

2 2 + + + 1 2 2 r 2 dr

(7.28)

93

and = 2 (1 2 ) when G=1 and a = a b is the dimensionless inner radius. The displacement type finite element representation of this very simple problem runs into a lot of difficulties when 0.5, i.e., when in the functional represented by Equation (7.28). These arguments now point to the volumetric strain field in Equation (7.28) as the troublesome term. When , this strain must vanish. In a physical situation, and in an infinitesimal interpretation of Equation (7.28), this means that 0. However, when a finite element discretization is introduced, this does not happen so simply. To see this, let us examine what will happen when a linear element is used to model the problem.

7.4.2 Formulation of a linear element

The displacement field u and the radial distance r are interpolated by linear isoparametric functions as

r = 0.5 (r1 + r2 ) + 0.5 (r2 r1 ) u = 0.5 (u1 + u2 ) + 0.5 (u2 u1 )

where the nodes of the element are at r1 and r2.

(7.29a) (7.29b)

The troublesome term in the potential energy functional is the volumetric strain and this contribution can be written as,

2 2

r dr =

(ru,r +2u )

dr

(7.30)

It is the discretized representation of this term

(ru,r +2u )

that will

determine whether there will be locking when . Thus, unlike in previous sections, the constraint is not on the volumetric strain but on r because of the use of spherical coordinates to integrate the total strain energy. It is therefore meaningful to define this as a pseudo-strain quantity =r. To see how this is inconsistently represented, let us now write the interpolation for this derived kinematically from the displacement fields as,

(ru,r +2u ) =

(r1 + r2 ) (u2 u1 ) + (u1 + u2 ) + [3 2 (u2 u1 )] 2 (r2 r1 )

(7.31)

Consider now the constraints that are enforced when this discretized field is used to develop the strain energy arising from the volumetric strain as shown in Equation (7.31) when .

(r1 + r2 ) ( u2 u1 ) + (u1 + u2 ) 0 2 (r2 r1 ) (u2 u1 ) 0

(7.32a)

(7.32b)

94

Condition (7.32a) represents a meaningful discretized version of the condition of vanishing volumetric strain and is called a true constraint in our field-consistency terminology. However, condition (7.32b) leads to an undesired restriction on the u field as it implies u,r 0 . This is therefore the spurious constraint that leads to locking and spurious pressure oscillations.

7.4.3 Incompressible 3D elasticity

We now go directly to modeling with 3D elasticity. We re-write the strain energy for 3D elasticity in the following form:

U =1 2G

d D d dV + 1 2 n n dV
T T

(7.33)

In Equation (7.33), d is the distortional strain and n is the volumetric strain. In this form, it is convenient to describe the elasticity matrix in terms of the shear modulus, G and the bulk modulus K, where G = E [2 (1 + )] ,

K = E [3 (1 - 2 )] and = (K 2G 3 ) = E

[(1 + ) (1 - 2 )] .

Consider the fields given in Equations (7.24) for the 8-noded brick element. An elastic response for incompressible or nearly-incompressible materials will require a constrained strain field,

n = u,x +v , y +w,z 0

(7.34)

where n is the volumetric strain. From Equations (7.24) and (7.34) we have,

n = (a2 + b3 + c4 ) + (b5 + c7 )x + (a5 + c6 )y + (a7 + b6 )z + (a8 yz + bz xz + c8 xy )

(7.35)

Equation (7.35) can be constrained to zero only when the coefficients of each of its terms vanish, giving rise to the constraint conditions,

a2 + b3 + c4 0 b5 + c7 0 a5 + c6 0 a7 + b6 0 a8 0 b8 0 c8 0

(7.36a) (7.36b) (7.36c) (7.36d) (7.36e) (7.36f) (7.36g)

Equation (7.36b) to (7.36) are the inconsistent terms that cause locking. To ensure a field-consistent representation, only the consistent terms should be retained, giving

n = a2 + b3 + c4
This field will now be able to respond to the incompressible or incompressible strain states and the element should be free of locking.

(7.37) nearly

95

7.4.4 Concluding remarks

It is clear from this section that incompressible locking is analogous to shear locking, etc. A robust formulation must therefore ensure consistent representations of volumetric strain in cases where incompressibility or near incompressibility is expected.

7.5 References
7.1 G. Prathap, The Finite Element Academic Press, Dordrecht, 1993.

Method

in

Structural

Mechanics,

Kluwer

96

Chapter 8 Stress consistency

8.1 Introduction
In Chapters 6 to 7 we examined the difficulties experienced by the displacement approach to the finite element formulation of problems in which some strain fields are constrained. To obtain accurate solutions at reasonable levels of discretization, it was necessary to modify these strain fields and use these in computing the stiffness matrix and also in stress recovery. The criterion governing the relationship between the various terms in the modified strain field interpolation was described as consistency - i.e. the strain field interpolations must maintain a consistent balance internally of its contributing terms. This allows the constraints that emerge after discretization to remain physically faithful to the continuum problem. This was the rule that guided the locking-free design of all elements discussed so far. In this chapter, we take a look at a class of problems where no constraints are imposed on the strains but there is a need to relax the satisfaction of the constitutive relationship linking discretised stress to discretised strain so that again a certain degree of consistency is maintained. Such situations develop where there are structural regions in which the rigidity varies spatially due to varying elastic moduli or cross-sectional area and also in initial strain problems, of which the thermal strain problem is the most commonly encountered. In structural regions with varying rigidity the spatial variation of strain-fields and stress or stress-resultant fields will not match. In the discretization of such cases, it is necessary to consider a form of external consistency requirement between the discretized strain fields and the discretized stress or stress-resultant fields. This is necessary so that a correct interpretation of computed stresses and stress resultants is possible; otherwise, oscillations will be seen. In initial strain problems, the initial strain variation and the total strain variations may be of different order. This is a familiar problem in finite element thermal stress analysis. Here, it is necessary to obtain kinematically equivalent thermal nodal forces from temperature fields and also ensure that the discretised thermal strains corresponding to this are consistent with the discretised description of total strains. Again, one must carefully identify the conflicting requirements on the order of discretised functions to be used.

8.2 Variable moduli problems 8.2.1 A tapered bar element

We shall now conduct a simple numerical experiment with a tapered bar element to show that the force field computed directly from strains in a structural element of varying sectional rigidities has extraneous oscillations. A linear element would have sufficed to demonstrate the basic principles involved. However, we use a quadratic tapered bar element so that the extraneous oscillations, which

97

for a general case can be of cubic form, are not only vividly seen but also need special care to be reduced to its consistent form. In a linear element, this exercise becomes very trivial, as sampling at the centroid of the element gives the correct stress resultant. We shall consider an isoparametric formulation for a quadratic bar element of length 2l with mid-node exactly at the mid-point of the element. Then the interpolations for the axial displacement u and the cross sectional area A in terms of their respective nodal values are,

u = u2 + (u3 u1 ) 2 + (u1 2u2 + u3 ) 2 2 A = A 2 + (A 3 A1 ) 2 +

(A1 2A 2 + A3 )

2 2

We shall first examine how the element stiffness is formulated when the minimum total potential principle is used. We start with a functional written as,

= 1 2 N T dx W
where,

= du dx N = E A (x )
W = pu dx p E

the axial strain, the kinematically constituted axial force, the potential of external forces, distributed axial load, Young's modulus of elasticity

After discretization, will be a linear function of but N will be a cubic function of . The strain energy of deformation is then expressed as,

U =1 2

dx

From this product the terms of the stiffness matrix emerge. Due to the orthogonal nature of the Legendre polynomials, terms from N which are linked with the quadratic and cubic Legendre polynomials, N3 and N4 respectively, will not contribute to the energy and therefore will not provide terms to the stiffness matrix! It is clear that, the displacements recovered from such a formulation cannot recognize the presence of the quadratic and cubic terms N3 and N4 in the stress field N as these have not been accounted for when the stiffness matrix was computed. Hence, in a displacement type finite element formulation, the stresses recovered from the displacement vector will have extraneous oscillations if N3 and N4 are not eliminated from the stress field during stress recovery. We shall designate by BAR3.0, the conventional element using N for stiffness matrix evaluation and recovery of force resultant. Next, we must see how the consistent representation of the force field denoted by N must be made. N should comprise only the terms that will contribute to the stiffness and strain energy and the simplest way to do this is to expand the kinematically determined N in terms of Legendre polynomials, and retain only terms that will meaningfully contribute to the energy in

(N ) .
T

98

Thus, N terms:

must be consistent with , i.e. in this case, retain only up to linear

N = N 1 + N 2

(8.1)

Such an element is denoted by BAR3.1. It uses N for stiffness matrix evaluation and recovery of forces resultant. To see the variational basis for the procedure adopted so far, the problem is re-formulated according to the Hu-Washizu principle which allows independent fields for assumed strain and assumed stress functions.

8.2.2 Numerical experiments

We shall perform the computational exercises with the two versions of the element; note that in both cases, the stiffness matrices and computed displacements are identical. Figure 8.1 shows a tapered bar clamped at node-1 and subjected to an axial force P at node-3. The taper is defined by the parameters, = (A 3 A1 ) 2A 2 and = (A1 2A 2 + A 3 ) 2A 2 Finite element results from a computational exercise using the two versions described above for a bar with cross section tapering linearly from the root to the tip for which =0 are obtained. Thus is given by,

= (A 3 A1 ) (A1 + A3 )
Thus can vary from 0 to 1.0. Fig. 8.2 shows the axial force patterns obtained from the finite element digital computation for a case with =-0.9802 (with A3=0.01 A1) It can be noted here that the results are accurate at

= 1

3.

Fig. 8.1 A cantilever bar modeled with a single element.

99

Fig. 8.2 Axial force pattern for linearly tapered bar (=0.9802 and =0.0) with A3=0.01 A1.

A general case of taper is examined next where A1=1.0, A2=0.36 and A3=0.04, and the area ratios are =-4/3 and =4/9. Fig. 8.3 shows the results from the finite element computations. Due to the presence of both quadratic and cubic oscillations, there are no points which can be easily identified for accurate force recovery! Therefore it is necessary to perform a re-constitution of the force resultant fields on a consistency basis using the orthogonality principle as done here before reliable force recovery can be made.

8.2.3 Reconstitution of the stress-resultant field using the Hu-Washizu principle

In forming the Hu-Washizu functional for the total potential, an as yet undetermined assumed force function N is introduced but the assumed strain field can be safely retained as (note that in a constrained media problem it will be required to introduce a field-consistent that will be different from the kinematically derived and therefore usually field-inconsistent ) - the functional now becomes

{1 2 (EA )

+ N T ( ) dx W

100

Fig. 8.3 Axial force pattern for bar with combined linear and quadratic taper (=-4/3 and =4/9).

A variation of the Hu-Washizu energy functional with respect to the kinematically admissible degree of freedom u, gives the equilibrium equation,
T u {dN dx p} dx = 0

Variation with respect constitutive relation

to

the

assumed

strain

field

gives

rise

to

{- N + EA } dx = 0
N

(8.2) gives rise to the (8.3)

and variation with respect to the assumed force field condition

{ - } dx = 0

Now Equations (8.2) and (8.3) are the orthogonality conditions required for reconstituting the assumed fields for the stress resultant and the strain. The consistency paradigm suggests that the assumed stress resultant field N

101

should be of the same order as the assumed strain field . Then Equation (8.3) gives the orthogonality condition for strain field-redistribution. On the other hand, orthogonality condition (8.2) can be used to from the kinematically reconstitute the assumed stress-resultant field N derived field N. Now, Equation (8.2) can be written as,

{N - N } dx = 0

(8.4)

Thus, if N is expanded in terms of Legendre polynomials, it can be proved that N which is consistent and orthogonally satisfies Equation (8.4) is obtained very simply by retaining all the Legendre polynomial terms that are consistent with , i.e. as shown in Equation (8.1). Thus the procedure adopted in the previous section has variational legitimacy.

8.3 Initial strain/stress problems

Finite element thermal stress analysis requires the formulation of what is called an initial strain problem. The temperature fields which are imposed must be converted to discretised thermal (initial) strains and from this kinematically equivalent thermal nodal forces must be computed. The usual practice in general purpose codes is to use the same shape functions to interpolate the temperature fields and the displacement fields. Thermal stresses computed directly from stress-strain and strain-displacement matrices after the finite element analysis is performed thus show large oscillating errors. This can be traced to the fact that the total strains (which are derived from displacement fields) are one order higher than the thermal strains (derived from temperature fields). Some useful rules that are adopted to overcome this difficulty are that the temperature field used for thermal stress analysis should have the same consistency as the element strain fields and that if element stresses are based on Gauss points, the thermal stresses should also be based on these Gauss point values. This strategy emerged from the understanding that the unreliable stress predictions originate from the mismatch between the element strain and the initial strain due to temperature 0. We shall now show that this is due to the lack of consistency of their respective interpolations within the element. Earlier in this chapter, we saw that stress resultant fields computed from strain fields in a displacement type finite element description of a domain with varying sectional rigidities showed extraneous oscillations. This was traced to the fact that these stress resultant fields were of higher interpolation order than the strain fields and that the higher degree stress resultant terms did not participate in the stiffness matrix computations. In this section, we show that much the same behavior carries over to the problem of thermal stress computations.

8.3.1 Description of the problem

With the introduction of initial strains 0 due to thermal loading, the stress to strain relationship has to be written as

102

 = D ( - 0 ) = D m

(8.5)

The strain terms now need to be carefully identified. {} is the total strain and {m} is the mechanical or elastic strain. The free expansion of material produces initial strains

0 = T

(8.6)

where T is the temperature relative to a reference value at which the body is free of stress and is the coefficient of thermal expansion. The total strains (i.e. the kinematically derived strains) are defined by the strain-displacement matrix,

{}=[B]{d}

(8.7)

where {d} is the vector of nodal displacements. In a finite element description, the displacements and temperatures are interpolated within the domain of the element using the same interpolation functions. The calculation of the total strains {} involves the differentiation of the displacement fields and the strain field functions will therefore be of lower order than the shape functions. The initial strain fields (see Equation (8.6)) involve the temperature fields directly and this is seen to result in an interpolation field based on the full shape functions. The initial strain matrix is of higher degree of approximation than the kinematically derived strain fields if the temperature fields can vary significantly over the domain and are interpolated by the same isoparametric functions as the displacement fields. It is this lack of consistency that leads to the difficulties seen in thermal stress prediction. This originates from the fact that the thermal load vector is derived from a part of the functional of the form,

{ }

D {0 } dV

(8.8)

Again, the problem is that the `higher order' components of the thermal (or initial) stress vector are not sensed by the total strain interpolations in the integral shown above. In other words, the total strain terms "do work" only on the consistent part of the thermal stress terms in the energy or virtual work integral. Thus, a thermal load vector is created which corresponds to a initial strain (and stress) vector that is `consistent' with the total strain vector. The finite element displacement and total strain fields which are obtained in the finite element computation then reflect only this consistent part of the thermal loading. Therefore, only the consistent part of the thermal stress should be computed when stress recovery is made from the nodal displacements; the inclusion of the inconsistent part, as was done earlier, results in thermal stress oscillations. We demonstrate these concepts using a simple bar element.

8.3.2 The linear bar thermal load vector

element

Derivation

of

stiffness

matrix

and

Consider a linear bar element of length 2l. The axial displacement u, the total strain , the initial strain 0 and stress are interpolated as follows:

u = (u1 + u2 ) 2 + (u2 u1 ) 2

(8.9a)

103

Fig. 8.4 Linear bar element.

= u,x = (u2 u1 ) 2l 0 = {(T1 + T2 ) 2 + (T2 T1 ) 2} = E (u2 u1 ) 2l E {(T1 + T2 ) 2 + (T2 T1 ) 2}

where

(8.9b) (8.9c) (8.9d)

=x l

(see Fig.8.4), E is the modulus of elasticity, A the area of

cross-section of the bar and the coefficient of expansion. To form the element stiffness and thermal load vector, an integral of the form

dx

(8.10)

is to be evaluated. This leads to a matrix equation for each element of the form,

EA 2l

1 1

1 1

u1 (T1 + T2 ) 1 F1 = EA u2 2 1 F2

(8.11)

where F1 and F2 are the consistently distributed nodal loads arising from the distributed external loading. By observing the components of the interpolation fields in Equations (8.9b) to (8.9d) carefully (i.e. constant and linear terms) and tracing the way they participate in the `work' integral in Equation (8.10), it is clear that the (T2 T1 ) 2 term associated with the linear (i.e. ) term in

(originating from 0) cannot do work on the constant term in T and therefore

vanishes from the thermal load vector; see Equation (8.11). Thus the equilibrium equations that result from the assembly of the element equilibrium equations represented by Equation (8.11) will only respond to the consistent part of the initial strain and will give displacements corresponding to this part only. If these computed displacements are to be used to recover the initial strains or thermal stresses, only the consistent part of these fields should be used. The use of the original initial strain or stress fields will result in oscillations corresponding to the inconsistent part. We shall work these out by

104

Fig. 8.5 Clamped bar subject to varying temperature hand using a simple example below and compare it with the analytical solution.

8.3.3 Example problem

Figure 8.5 shows a bar of length L=4l clamped at both ends and subjected to a varying temperature field. We shall consider a case where two conventionally derived linear elements are used, so that we require the nodal temperatures T1, T2 and T3 as input. The nodal reactions are F1 and F3 (corresponding to the clamped conditions u1=u3=0). We have the assembled equations as,

1 1 EA 1 2 2l 1

1 1

0 F1 EA u2 = 0 + 2 0 F 3

(T1 + T2 ) (T1 T3 ) (T + T ) 2 3

(8.12)

From this, we can compute the displacements and nodal reactions as,

u2 = l (T1 T3 ) 2 and F1 = F3 = EA (T1 + 2T2 + T3 ) 4

(8.13)

and these are the correct answers one can expect with such an idealization. If the stress in the bar is computed from the nodal reactions, one gets a constant stress = E (T1 + 2T2 + T3 ) 4 , in both elements, which is again the correct answer one can expect for this idealization. It is a very trivial exercise to show analytically that a problem in which the temperature varies linearly from 0 to T at both ends will give a constant stress field = ET 2 , which the above model recovers exactly. Problems however appear when the nodal displacement computed from Equations (8.12) is used in Equations (8.9b) to (8.9d) to compute the initial strains and stresses in each element. We would obtain now, the stresses as (subscripts 12 and 23 denote the two elements)

12 = E (T1 + 2T2 + T3 ) 4 E (T2 T1 ) 2 23 = E (T1 + 2T2 + T3 ) 4 E (T3 T2 ) 2

(8.14a) (8.14b)

It is now very clear that a linear oscillation is introduced into each element and this corresponds to that inconsistent part of the initial strain or stress interpolation which was not sensed by the total strain term. This offers us a straightforward definition of what consistency is in this problem - retain only that part of the stress field that will do work on the strain term in the

105

functional. To see how this part can be derived in a variationally correct manner, we must proceed to the Hu-Washizu theorem.

8.3.4 A note on the quadratic bar element

We may note now that if a quadratic bar element had been the basis for the finite element idealization, the total strains would have been interpolated to a linear order; the initial strain and thermal stress field will now have a quadratic representation (provided the temperature field has a quadratic variation) and the inconsistency will now be of the 1 3 2 type; thus thermal stresses derived using a formal theoretical basis will show these quadratic oscillations which will vanish to give correct answers at the points corresponding to = 1 integration rule.

3 ; i.e. the points corresponding to the 2-point Gauss

8.3.5 Re-constitution of the thermal strain/stress field using the HuWashizu principle
We now seek to find a variational basis for the use of the re-constituted consistent initial strain and stress interpolations in the Hu-Washizu principle. The minimum total potential principle states the present problem as, find the minimum of the functional,

MTP =

m 2 + P dV

(8.15)

where and m are as defined in (8.5) to (8.7) and the displacement and strain fields are interpolated from the nodal displacements using the element shape functions and their derivatives and P is the potential energy of the prescribed loads. We now know that the discretized stress field thus derived, , is inconsistent to the extent that the initial strain field 0 is not of the same order as the total strain field . It is therefore necessary to reconstitute the discretized stress field into a consistent stress field without violating any variational norm. In the examples above, we had seen a simple way in which this was effected. To see how we progress from the inconsistent discretized domain (i.e. involving 0 and ) to the consistent discretized domain (i.e. introducing 0 and , it is again convenient to develop the theory from the generalized HuWashizu mixed theorem. We shall present the Hu-Washizu theorem from the point of view of the need to re-constitute the inconsistent 0 to a consistent 0 without violating any variational norms. We proceed thus: Let the continuum linear elastic problem have a discretized solution based on the minimum total potential principle described by the displacement field u, strain field and stress field (we project that the strain field is derived from the displacement field through the strain-displacement gradient operators

106

are computed using the constitutive relationships, i.e. = D ( 0 ) . It is clear that 0 is an approximation of the strain field 0. Note that we also argue that we can use = as there is no need to introduce such a distinction here (in a constrained media elasticity problem it is paramount that be derived as the consistent substitute for .) fields What the Hu-Washizu theorem does, following the interpretation given by de Veubeke, is to introduce a "dislocation potential" to augment the usual total potential. This dislocation potential is based on a third independent stress field which can be considered to be the Lagrange multiplier removing the lack of compatibility appearing between the kinematically derived strain field 0 and the independent strain field 0 . The three-field Hu-Washizu theorem can be stated as,

of the theory of elasticity and that the stress field is derived from the strain field through the constitutive laws as shown in (8.2). Let us now replace the discretized domain by another discretized domain corresponding to the application of the Hu-Washizu principle and describe the computed state to be defined by the quantities , 0 and , where again, we take that the stress

HW = 0
where

(8.16)

HW = T m 2 +

( m

m ) + P dV

(8.17)

where m = ( 0 ) . At this stage we do not know what they are to be of consistent order with .

or 0 are except that

In the simpler minimum total potential principle, which is the basis for the derivation of the displacement type finite element formulation in most textbooks, only one field (i.e. the displacement field u), is subject to variation. However, in this more general three field approach, all three fields are subjected to variation and leads to three sets of equations which can be grouped and classified as follows: Variation on Nature Equilibrium Orthogonality (Compatibility) Equation,

+ terms from P=0

(8.18a) (8.18b)

(0 0 ) dV

=0

Orthogonality (Equilibrium)

dV = 0

(8.18c)

107

Fig. 8.6 (a) Clamped bar under linear temperature variation and its (b) Bar element model, (c) Plane stress model. Let us first examine the orthogonality condition in (8.18c). We can interpret this as a variational condition to restore the equilibrium imbalance between and . In this instance this condition reduces to = . Note that in a problem where the rigidity modulus D can vary significantly over the element volume, this condition allows consistent way.

to be reconstituted from

in a

The orthogonality condition in (8.18b) is now very easily interpreted. Since we have = consistent with , this condition shows us how to smooth 0 to 0 to maintain the same level of consistency as . This equation therefore provides the variationally correct rule or procedure to determine consistent thermal strains and stresses from the inconsistent definitions.

Therefore, we now see how the consistent initial strain field 0 can be derived from the inconsistent initial strain field 0 without violating any variational norm. We thus have a variationally correct procedure for reconstituting the consistent initial strain field - for the bar element above, this can be very trivially done by using an expansion of the strain fields in terms of the orthogonal Legendre polynomials, as done in the previous section for the consistent stress-resultants.

108

8.3.6 Numerical examples

In this section we take up the same example of a bar with fixed ends subjected to linearly varying temperature along its axis and model it with the bar element and a plane stress element for demonstrating the extraneous stress oscillations resulting from the lack of consistency in initial strain definition. A bar of length 8 units, depth 2 units subjected to a linear temperature distribution as shown in Fig 8.6a is modeled first with two bar elements (see Fig 4. 18.6b). Let BAR.0 and BAR.1 represent the bar element versions with inconsistent and consistent thermal stress fields. As already predicted, both give accurate displacements (-0.002 units at the mid point), see Equation (9.13)); BAR.1 gives exact stress throughout while BAR.0 shows linear oscillations (see Equations (8.14a) and (8.14b)) as shown in Fig 8.7.

Fig.8.7 Consistent and inconsistent thermal stress recovery for a clamped bar problem.

109

8.3.7 Concluding remarks

Before we close this section, it will be worthwhile to discuss the "correctness" of the various approaches from the variational point of view. Note that Equation (9.2) and (9.3) would have been fulfilled if an assumed stress-resultant field of the same order as N had been used in place of the N of consistent order. This would result in the Hu-Washizu formulation yielding the same stress predictions as the minimum potential principle. Thus, both N and N are equally "correct" in the Hu-Washizu sense, but only N is "correct" with the minimum potential energy formulation. Where they differ is in the consistency aspect i.e. N is consistent with whereas N isn't. It follows that from the theoretical point of view of the variational or virtual work methods, the potential energy formulation and the Hu-Washizu formulation describe the continuum physics in exactly the same way. Both the potential energy and HuWashizu formulations are equally valid as far as variational correctness is concerned. However, when approximations are introduced in the finite element method, the Hu-Washizu formulation gives decidedly better results than the potential energy formulation because of the flexibility it provides to modify stress fields and strain fields to satisfy the consistency requirements in situations where they play a prominent role. Note that the consistency paradigm, in this case that used to justify the truncation of N (obtained directly from the displacement approach) to a consistent N on the argument that the truncated terms are not sensed in the discretized strain energy computations, lies clearly outside the variational correctness paradigm. Similarly, in the earlier chapters of this book, we saw another variation of the consistency paradigm - the need to ensure a proper balance in the discretized representation of constrained strain fields, which again lies outside the scope of the variational correctness paradigm. It is very important to understand therefore that both consistency and correctness paradigms are mutually exclusive but are needed together to ensure that the finite element formulations are correctly and consistently done. The stress field-consistency paradigm introduced here also applies to finite element applications where certain terms in either the stress or strain fields do not participate in the displacement recovery due to their inconsistent representation e.g. initial stress/strain problems, problems with varying material properties within an element etc. In the next section, we shall extend this paradigm to extract consistent thermal stress and/or strains from a displacement type formulation. In this section, we have demonstrated another variation of the consistency paradigm - in applications to evaluation of stresses and stress resultants in finite element thermal stress computations. It is shown that such stresses must be computed in a consistent way, recognizing that the displacement type finite element procedure can sense only stresses which have the same consistency as the total strain field interpolations used within each element. Thus, in an element where the temperature field interpolations do not have the same consistency as the strain fields, it is necessary to derive a consistent initial strain-field for purposes of recovery of stresses from the computed displacement field. A simple and variationally correct way to do this has been established from the Hu-Washizu theorem.

110

Chapter 9 Conclusion
9.1 Introduction
Technology and Science go hand in hand now, each fertilizing the other and providing an explosive synergy. This was not so in the beginning; as there was Technology much before there was Science. In fact, there was technology before there were humans and it was technology which enabled the human race to evolve; science came later. The FEM is also an excellent example of a body of knowledge that began as Art and Technology; Science was identified more slowly. This chapter will therefore sum up the ideas that constitute the science of the finite element method and also point to the future course of the technology.

9.2 The C-concepts: An epistemological summing-up

One would ideally like to find a set of cardinal or first principles that govern the entire finite element discretisation procedure, ensuring the quality of accuracy and efficiency of the method. At one time, it was believed that two rules, namely continuity and completeness, provided a necessary and sufficient basis for the choice of trial functions to initiate the discretisation process. Chapter 3 of this book reviewed the understanding on these two aspects. However, finite element practitioners were to be confronted with a class of problems where the strict implementation of the continuity condition led to conflicts where multiple strain fields were required. This was the locking problem which we examined carefully in chapters 5 to 7. Let us recapitulate our fresh insight in an epistemological framework below.

9.2.1 Continuity conflict and the consistency paradigm

We saw in chapters 5 to 7 that in a class of problems where some strain fields need to be constrained, the implementation of the continuity rules on the displacement trial functions blindly led to a conflict whereby excessive constraining of some strain terms caused a pathological condition called locking. Let us again take up the simple problem of a Timoshenko beam element. Two strain energy components are important; the bending energy which is constituted from the bending strain as = ,x and the shear strain energy which is based on the shear strain = w,x . The simplistic understanding that prevailed for a

,x and w,x terms in the energy expression demanded that the trial functions selected for the and w fields 0 must be C continuous; i.e. continuity of and w at nodes (or across edges in a 2D problem, etc.) was necessary and sufficient.
long time was that the existence of the However, this led to locking situations. To resolve this conflict, it became necessary to treat the two energy components on separate terms. The shear strain energy was constrained in the regime of thin beam bending to become vanishingly small. This required that the shear strains must vanish without

111

producing spurious constraints. The idea of field consistent formulation was to anticipate these conditions on the constrained strain field. Thus for w,x to vanish without producing gratuitous constraints, the trial functions used for and w in the shear strain field alone must be complete to exactly the same order - this is the consistency paradigm. We can also interpret this as a multiple continuity requirement. Note that although the bending strain requires to be 0 0 C continuous, in the shear strain component, we require not to have C -1 continuity, i.e. say a C continuity, is. Thus, such a conflict is resolved by operating with consistently represented constrained strain fields.

9.2.2 The correctness requirement

The consistency requirement now demands that when trial functions are chosen for displacement fields, these fields (such as in = w,x ) which need to be constrained must have a reconstituted consistent definition. The minimum total potential principle, being a single field variational principle, does not provide a clue as to how this reconstitution of the strain fields can be performed without violating the variational rules. In the course of this book, we have come to understand that a multi-field variational principle like the generalized Hu-Washizu theorem provided exactly that flexibility to resolve the conflicting demands between consistency and continuity. This, we called the correctness requirement - that the reconstitution of the strain field which needed a reduced-continuity satisfying displacement field from an inconsistent strain field kinematically derived from the original continuous displacement fields must be done strictly according to the orthogonality condition arising from the HW theorem.

9.2.3 The correspondence principle

The study so far with consistency and correctness showed the primacy of the HuWashizu theorem in explaining how the finite element method worked. It also became very apparent that if the internal working of the discretisation procedure is examined more carefully, it turns out that it is strain and stress fields which are being approximated in a "best-fit" sense and not the displacement fields as was thought earlier. This brings us to the stress correspondence paradigm, and the discovery that this paradigm can be axiomatised from the Hu-Washizu theorem, as shown in chapter 2.

9.3 How do we learn? How must we teach? In a chronological sense or in a logical sense?
In an epistemological sense, often, it is phenomena which first confronts us and is systematically recognized and classified. The underlying first principles that explain the phenomena are identified very much later, if not at all. The clear understanding of first principles can also lead to the unraveling of new phenomena that was often overlooked. Thus, sequence in and derive pedagogical if and when the first principles are available, then the logical which one teaches a subject will be to start with first principles the various facets of the phenomena from these principles. In practice, we can therefore have a choice from two courses of action:

112

teach in a chronological sense, explaining how the ideas unfolded themselves to us, or teach in the logical sense, from first principles to observable facts. It is not clear that one is superior to the other. In my first book [9.1], I choose the chronological order. The conflicts created by the locking phenomena were resolved by inventing the consistency paradigm. This led to the correctness principle and the recognition that the HuWashizu principle formed a coherent basis for understanding the fem procedure. From these, one could axiomatise the correspondence rule - that the fem approach manipulates stresses and strains directly in a best-approximation sense. Thus, chapter 12 of Reference 9.1 sums up the first principles. In the present book, I have preferred a logical sequence - the HW and correspondence rules are taken up first and then the various phenomena like locking, etc.

9.4 Finite element technology: Motivations and future needs

Advances in computer science and technology have had a profound influence on structural engineering, leading to the emergence of this new discipline we call computational structural mechanics (CSM). Along with it a huge software industry has grown. CSM has brought together ideas and practices from several disciplines - solid and structural mechanics, functional analysis, numerical analysis, computer science, and approximation theory. CSM has virtually grown out of the finite element method (FEM). Algorithms for the use of the finite element method in a wide variety of structural and thermomechanical applications are now incorporated in powerful general purpose software packages. The use of these packages in the Computer-AidedDesign/Computer-Aided-Manufacturing cycle forms a key element in new manufacturing technologies such as the Flexible Manufacturing Systems. These allow for unprecedented opportunities for increase in productivity and quality of engineering by automating the use of structural analysis techniques to check designs quickly for safety, integrity, reliability and economy. Very large structural calculations can be performed to account for complex geometry, loading history and material behavior. Such calculations are now routinely performed in aerospace, automotive, civil engineering, mechanical engineering, oil and nuclear industries. Modern software packages, called general purpose programmes, couple FEM software with powerful graphics software and the complete cycle of operations involving pre-processing (automatic description of geometry and subsequent subdivision of the structure) and post-processing (projecting derived information from FEM analysis on to the geometry for color coded displays to simplify interpretation and make decision making that much easier). Already artificial intelligence in the form of knowledge based expert systems and expert advisers and optimization procedures are being coupled to FEM packages to reduce human intervention in structural design to a bare minimum. It is not difficult to delineate many compelling reasons for the vigorous development of CSM [9.2]. These are:

113

1. There are a large number of unsolved practical problems of current interest which still await experimental and/or numerical solutions. Some of these demand large computational power. Some of the examples described in reference 9.2 are: simulation of response of transportation vehicles to multidirectional crash impact forces, dynamics of large flexible structures taking into account joint nonlinearities and nonproportional damping, study of thermoviscoelatic response of structural components used in advanced propulsion systems etc. In many structural problems, the fundamental mechanics concepts are still being studied (e.g. in metal forming, adequate characterization of finite strain inelasticity is still needed). 2. Computer simulation is often required to reduce the dependence on extensive and expensive testing; in certain mission critical areas in space, computer modeling may have to replace tests. Thus, for large space structures (e.g. large antennas, large solar arrays, the space station), it may not be possible for ground-test technology in 1-g environment to permit confident testing in view of the large size of the structures, their low natural frequencies, light weight and the presence of many joints. 3. Emerging and future computer systems are expected to provide enormous power and potential to solve very large scale structural problems. To realize this potential fully, it is necessary to develop new formulations, computational strategies, and numerical algorithms that exploit the capabilities of these new machines (e.g. parallelism, vectorization, artificial intelligence). Noor and Atluri [9.2] also expect that high-performance structures will demand the following technical breakthroughs: 1. Expansion of the scope of engineering problems modeled, such as: a) examination of more complex phenomena (e.g. damage tolerance of structural components made of new material systems); b) study of the mechanics of high-performance modern materials, such as metalmatrix composites and high-temperature ceramic composites; c) study of structure/media interaction phenomena (e.g. hydrodynamic/structural coupling in deep sea mining, thermal/control/structural coupling in space exploration, material/aerodynamic/structural coupling in composite wing design, electromagnetic/thermal/structural coupling in microelectronic devices); d) use of stochastic models to account for associated with loads, environment, and material variability e) development of efficient high-frequency nonlinear dynamic modeling capabilities (with applications to impulsive loading, high-energy impact, structural penetration, and vehicle crash-worthiness); f) improved representation of structural details such as damping and flexible hysteritic joints: g) development of reliable life-prediction methodology for structures made of new materials, such as stochastic mechanisms of fatigue, etc.

114

h) analysis and design of intelligent structures with active and/or passive adaptive control of dynamic deformations, e.g. in flight vehicles, large space structures, earthquake-resistant structures i) Computer simulation of manufacturing processes interface mechanics, superplastic forming. such as solidification,

2.

Development of practical measures for assessing the reliability of the computational models and estimating the errors in the predictions of the major response quantities. Continued reduction of cost and/or engineering design/analysis problems. time for obtaining solutions to

3.

Special hardware and software requirements must become available to meet the needs described above. These include: 1. Distributed computing environment having high-performance computers for large scale calculations, a wide range of intelligent engineering workstations for interactive user interface/control and moderate scale calculations.

2. User-friendly engineering workstations with high-resolution and high speed graphics, high speed long distance communication etc. 3. Artificial intelligence-based expert systems, incorporating the experience and expertise of practitioners, to aid in the modeling of the structure, the adaptive refinement of the model and the selection of the appropriate algorithm and procedure used in the solution. 4. Computerized symbolic manipulation capability calculations and increase their reliability. to automate analytic

5. Special and general purpose application software systems that have advanced modeling and analysis capabilities and are easy to learn and use.

9.5 Future directions

It is expected that CSM will continue to grow in importance. Three areas which are expected to receive increasing attention are: 1) modeling of complex structures; 2) predata and postdata processing and 3) integration of analysis programs into CAD/CAM systems. The accurate analysis of a complex structure requires the proper selection of mathematical and computational models. There is therefore a need to develop automatic model generation facilities. Complex structures will require an enormous amount of data to be prepared. These can be easily performed by using predata and postdata processing packages using high resolution, high throughput graphic devices. The generation of three dimensional color movies can help to visualize the dynamic behavior of complex structural systems.

115

9.6 Concluding remarks

In this concluding chapter, we have summarized the ideas that should inform the rational development of robust finite elements and we have also described the current trends in the use of structural modeling and analysis software in engineering design. It has been seen that CSM has greatly improved our capabilities for accurate structural modeling of complex systems. It has the potential not only to be a tool for research into the basic behavior of materials and structural systems but also as a means for design of engineering structures.

9.7 References
9.1 G. Prathap, The Finite Element Method in Structural Mechanics, Kluwer Academic Press, Dordrecht, 1993. 9.2 A. K. Noor and S. N. Atluri, Advances and trends in computational structural mechanics, AIAA J, 25, 977-995 (1987).

116