- Email: [email protected]
Critical features of the finite element method
CriticalFeaturesof the Finite Element Method by R.
J. MELOSH
MARC
Analysis
Research
Corporation,
Palo Alto,
California
This paper examines the essential features of finite element analysis and identifies the inherent limitations of the method. It defines those aspects of finite element analysis which are common to all applications of the method. It discusses and compares the requirements and form of element models for computer simulation of structural networks and continua. It generalizes the terminology to clarify use of the method for a variety of physical systems. The perspective of the paper is that the finite element method is a discipline and problem invariant technique for synthesizing models of any engineering system; that it encompasses a wide variety of analysis objectives including developing exact, approximate, and bounding response predictions; that it can be fleshed into a variational, weighted residual, finite differenence, or a semi-empirical analysis method; and that its intrinsic limitations on fidelity rest on the relevance of physical measurements and constitutive laws hung on its framework.
ABSTRACT:
I.
Introduction
Templeton (1)observes that the finite element method is the most significant advance in engineering analysis since the 1920’s. This is at once a general recognition of the increased power and scope of engineering analysis credited to the method and a particular testimony to its practical value in application to the complex structural and heat transfer problems of the automobile industry. This dramatic tribute belies the few and simple ideas which are the skeleton of the method. Chief among these is the concept that a mathematical model of any engineering system can be synthesized from models of disjoint parts. The fact that this concept inbues an analysis modularity well-suited to digital computer implementation is of major importance to widespread use of the method. The consequence that the framework of finite element analysis is independent of engineering discipline, geometric and topological complexity, and response intricacy is the source of its power. The engineering, physics, and mathematical literature project the finite element method as a variational or weighted residual method. It expresses it in scalar, matrix, or tensor form involving differential, integral, or difference equations. It relates it to networks, continua and compound systems. It describes finite element analysis in the context of structures, elasticity, heat transfer, electricity, hydrodynamics, and/or magnetism. All of which may lead to some confusion about what the method is and what its limitations are. This paper seeks to identify the skeleton of finite element analysis and establish its inherent range of application. Most of the paper presents the method as a method of analysis of structural systems which behave as linear 489
R. J. Melosh
systems. It provides an analysis of the method from a phenomenological viewpoint: i.e. it evokes limitations of the method as reflected in the numerics which represent it. The presentation moves from a description of the finite element skeleton to study of its bones. Section II reviews the basic steps of analysis. Section III and IV consider analysis of networks and continua and the relation between these analysis classes. Section V replaces structural concepts with more abstract ones in identifying the interdisciplinary scope of the method. The last section briefly summarizes the key ideas of the paper. Though sparks of the finite element method can be found in the works of Courant (2), Prager (3), Synge (4) and others, the flames were kindled when Turner et al. (5) put the method in practical perspective. Since then, the ideas have undergone many transformations, but it is to the Turner work that this paper owes its largest debt. II. The Finite Element
Analysis Method
Regarding the finite element method as a means for simulating structures, this section reviews and discusses implications of its synthesis process. The building block of the analysis is the finite element model. The model pertains to a disjoint region of the system. It relates forces at nodes on the boundary of the element to corresponding nodal displacements, kii=U
(I)
where k is element stiffness matrix, a square matrix, ii is the vector of displacements of the nodes and U is the vector of nodal forces. The vectors ii and U are cogradient. Thus, if the axes to which they are referenced are rotated and translated, they are expressed in the new coordinate system by ij=Rii
Q=RU
where Q and Q are the vectors of joint displacements new axes. In the new basis, Eq. (1) is given by RkR-‘ij
(2) and forces related to the
= Q.
(3)
The potential energy, V, for the finite element model is defined as v=~liku-iiU=~~RkR-‘q-~Q
(4)
where the terms involving the 4 factor represent the strain energy and those involving joint forces represent external work. Thus, Eq. (1) defines a configuration state which is a stationary value of the potential. Equation (4) defines a scalar functional in which the element equations are embedded. Equations of two or more elements include equations of each element and constraint equations requiring that nodes common to the two elements have common values for corresponding components of displacement. Thus, if the 490
Journal
of The
Franklin
Institute
Critical Features of the Finite Element Method equations
for two elements
are represented
by
(5) the constraint
equations
can be expressed
by
0;! =
[mj,)
(6)
where C is the constraint matrix, iz is the vector of distinct displacement components as represented by & and &. Substituting Eqs. (5) and (6) in (4), and requiring that the function be stationary with respect to displacement changes yields the equations for the pair of elements in the form, C’KCiiz
= Qz
with
(Qz>= [C](!f).
Here, the QX are energy equivalent joint forces compliant with the continuity constraint. Introduction of constraints involves post- and pre-multiplication by the constraint matrix and its transpose as long as a potential exists. Some distinction between these transformations and coordinate transformations is lost if orthogonal coordinate systems are involved for then, R-r = R’.
(8)
Boundary conditions can be parsed in the form of Eq. (6). Some of the 41 will be constants if displacements are prescribed. If the values of some element joint forces are prescribed, coefficients of the constraint matrix can be abstracted from an element model equations, of Eqs. (l), thus defining a constraint. The equations that result from combining equations from all elements are the system load-deflection relations. These are solved for the unknown deflections. With all deflections known, the element equations (1) furnish the nodal forces on each element. The finite element synthesis process accommodates structural systems: systems which involve only equilibrium and displacement continuity conditions across disjoint elements. The behavior characteristics of the structure reside only in the element models and interelement nodal displacement constraints. The stiffness matrix for the complete system reflects the character of element matrices. With an appropriate choice of coordinate axes (e.g. orthogonal), both types of transformations are of contragredient form. Thus, if all element matrices are symmetric, the system matrix is. If all element matrices are positive semi-definite, the system stiffness matrix must also be positive semidefinite or positive-definite as a direct consequence of the definition of definiteness. The sparseness of the system stiffness matrix is well-known. In the usual case, the non-zero-coefficients in any column of the stiffness matrix can occur
Vol. 302, Nos. 5 & 6, November/December
1916
491
R. J. Melosh
only between a degree of freedom and its topologically closest neighbors. This sparseness is a direct consequence of the zero coefficients in the interelement constraint matrix of Eqs. (6). The requirement that finite elements represent disjoint substructures limits the variables included in Eqs. (6) to displacements for shared nodes and thus enforces sparseness of the system stiffness matrix. It is noted that Eqs. (1) through (8) are the basis for non-linear as well as linear analysis, despite the linear form of Eq. (1). These analyses involve tracking behavior through a sequence of loading steps each of which is assumed to involve a conservative system. In most formulations, each step involves an incremental loading and linearized element models (6)-(9). In some, nonlinearity is admitted (lo), (11). III. Analysis
of Networks
A network element is a finite element which represents uniquely and without mathematical approximation all admissible modes of behavior of the disjoint region it models. A network is a collection of network elements. A structural network is any physical system satisfactorily modeled as a network. Based on these definitions, the next paragraphs review the mathematrical requirements of network elements and the characteristics of network modeling. Figure 1 depicts some structural networks, networks, and network elements. These include rigid link-spring, torque tube, truss, and frame structures. These networks exhibit features common to all networks. Network elements involve one or two nodes. Any two nodes in the network can be connected by only one element. To distinquish networks from continua, network elements are called “branches” and nodes, “joints”. Since a branch must represent all structural characteristics without approximation, its mathematical model must be (a) based on displacement variables which can represent all possible boundary deformations, (b) include all elastic deformation states within the element, and (c) contain an explicit statement of branch equilibrium. Violation of any of these requirements disqualifies a finite element as a branch. For example, if only lateral displacements are used as nodal variables for a beam element, it cannot be a branch because the assumption of the beam equations that plane sections remain plane can not be enforced across element boundaries with this beam model. Satisfaction of requirements (a) and (b) results in load-deflection relations for the branch, corresponding to Eq. (l), which involve a symmetric positivedefinite stiffness matrix. Symmetry arises because the model is linear and conservative. It expresses the Maxwell-Betti theorem and is reflected in the self-adjoint differential equations of the beam element. Positive definiteness is a consequence of the fact that loads do work upon the structure in inducing deformation. Branch load-deflection relations, of the form of Eqs. (l), which involve only elastic deformations, can be modified to require satisfaction of macroscopic
492
Journal of The Franklin
Institute
Critical Features of the Finite Element Method
Mechanical
Drive
system
shafl
Tuss
Cd)
H Frame Structural
network
Network
element
Networks
FIG. 1. Structural networks, network, and network elements. equilibrium.
The equilibrium
requirements U,=CU,
(9)
where U, are undefined joint forces (reactions), augment the potential energy function. Minimization of this function with respect to the unknown displacements and the LaGrange multipliers then produces the equations:
(10) where the LaGrange multipliers have been interpreted as displacements in the reaction degrees of freedom. It is noted the linear dependence of the last r of Eqs. (10) signal a degeneracy of the element stiffness matrix reducing the rank by as many equations as the number of independent equilibrium requirements. The primary consequences of the branch characteristics are that the system stiffness matrix must be symmetric, positive semi-definite, and imply satisfac- . tion of equilibrium for each branch.. The symmetry is exploited to reduce the data processing task. With constraints introduced to eliminate system motions Vol. 302, Nos. 5 & 6, November/December
1976
493
R. J. Melosh
which involve no branch deformations, a positive definite stiffness matrix remains to insure a unique and minimum potential solution. Joint and branch equilibrium are often used as the basis for checks on combining of branch equations and solving of the system equations, despite unsatisfactory numerical characteristics of this test. Branch equations can be developed analytically or empirically. Analytical development involves defining a mathematical model and establishing its parameters by testing. For example, a differential equation can be used to model the beam branch. This model is particularized by measurements of the Young’s and Shear moduli for the material. To avoid mathematical approximation, exact solutions of the beam equation can be used in compliance with the definition of stiffness coefficients as joint forces associated with unit and zero joint displacements. Then the accuracy of the branch model depends on the accuracy of measurements of elastic constants. Network branch parameters can be uniquely particularized from physical tests on the structural system. Tests may involve measurements of displacements alone. This follows directly from the fact that since any two joints can be connected by a single branch, the coefficients of every non-diagonal partition of Eq. (10) can be uniquely defined by tests of the complete system. The remaining coefficients of Eqs. (10) can then be developed from the equilibrium requirements for each branch. For the frame of Fig. l(d), for example, displacements and rotations of each of the joints can be taken for a unit force at each joint and for a unit moment at each joint. Inversion of this flexibility matrix produces the system stiffness matrix. The off-diagonal joint partition for joints i and i (a 2 X 2 matrix) provides data for evaluating the three independent elastic coefficients needed for the load-deflection coefficients of Eqs. (1) for branch i-j. Similarly, equations for every other branch can be constructed without concern for the details of geometry and material composition. The scope of network analysis is limited by the validity of the form of the branch equations, Eqs. (1) and the accuracy of the measured displacements from which equation coefficients are deduced. Because branch equations cannot involve mathematical approximations, application is usually restricted to structures, like those of Fig. 1, whose topology articulates slender branches. IV. Analysis
of Continua
and Compound Systems
A continuum element is a finite element which approximates behavior of a disjoint part of a system. A continua is a collection of continuum elements. A structural continuum is a structural system which is satisfactorily modeled as a continuum. A compound system consists of both network and continuum elements. Using these definitions, the next paragraphs review mathematical requirements for continua and the characteristics of structural continua and compound systems. Figure 2 shows some structural continua, continuum elements, and a continua. These include membrane, plate, shell, and three-dimensional solid 494
Journal of The Franklin Institute
CriticalFeatures of the Finite Element Method
(0)
Membrane El---8
(b)
Plate
(cl
Shell
(d)
Solid
(a) Structural
continua
(b)Contlnuum
element
(c)Continua
FIG. 2. Structural continua, continuum elements, and continua. systems. Continua engage finite elements which usually have three or more nodes. A particular pair of nodes may be connected by one or more elements. As a minimum, the continuum element must satisfy the requirements for a branch as the element becomes vanishingly small. Thus, it must (a) involve boundary displacements which imply satisfaction of displacement continuity across element boundaries, (b) represent the stress-strain relations, and (c) imply satisfaction of microscopic equilibrium conditions everywhere in the element interior and across boundaries. Since the differential equations of elasticity can be cast in self-adjoint form, stiffness matrices must approach symmetry as the element size approaches zero. Since stress-strain relations of most structural material represent Green’s materials, stiffness matrices are usually positive semi-definite or positive definite. Conversely, elements of finite size displacements need not define all possible modes of boundary behavior, be continuous across boundaries, satisfy macroscopic equilibrium or imply satisfaction of compatibility within an element. Stiffness matrices need not be symmetric nor positive definite. An interesting aspect of finite element analysis of continua is the evolution of ad hoc criteria to insure the adequacy of element models. These criteria are phenomenological-they qualify the element matrix without regard for the means used in establishing the stiffness coefficients. Of prime importance is Iron’s patch test (12). This test is based on the proven hypothesis (13) that if deformations associated with a possible constant strain state are imposed on the boundaries of a set of the elements and the Vol. 302, Nos. 5 & 6, November/December
1976
495
R. J. Melosh
analysis yields the exact solution, then the element guarantees to lead to the solution of the equations of elasticity for any structure, when enough finite elements are used. In principle, this test must be used for all admissible constant states (e.g. x-strain, y-strain, and shear strain for a membrane), all element shapes, and all material models. This criteria brings respectability to many of the finite element models developed early in the evolution finite element technology. Irons (14) shows that with this criteria in mind, one can generate efficient element models (high accuracy/element) without concern for more sophisticated criteria. Taylor et al. (15)illustrate addition of behavior states to an element with regard for this criteria, to improve its efficiency. The patch test is closely related to an observation by Key (16) that adequacy of an element is related to its ability to reflect the modes of behavior associated with its matrix of elastic constants. For example, the material matrix for a membrane involves dilatation, stretching along one axis and contracting along another, and shear. A linear combination of those modes transforms them to constant strain states. The data of Table 1 illustrates that a stronger statement can be made-the validity and efficiency of an element model depends on the fidelity with which it represents the modes of the material matrix when the mesh is sufficiently fine. The uniformly thick membrane shown in Table 1 is articulated by 35 rectangular elements. Analyses for displacements were made for the four elements defined-each analysis involving 64 equilibrium equations. The mean error in displacements was determined by comparing solution results with an analysis with 1050 equations. The ratio of the maximum to minimum eigenvalue (the spectral norm) was evaluated as a scalar measure of the element matrices involved. The elements are listed in accordance with the closeness of the norm to the material matrix norm. The fact that this ordering corresponds exactly with the solution error suggests that when many elements are used, the element model whose natural behavior (as reflected by its eigenvalues and vectors) most closely corresponds with microscopic material natural modes will be most efficient. Recently a numerical sufficiency test has been proposed to determine whether an element model, given in the form of its stiffness matrix, must involve monotonic convergence (20). This test also checks that element models satisfy macroscopic equilibrium requirements. Thus, in combination with the patch test, there exists at least one set of numerical tests for qualifying a continuum element stiffness matrix. Can non-destructive tests of a structural continuum be used to evaluate coefficients of an acceptable continuum element stiffness matrix? A simple “no” answer does not do justice to the question. The no answer can be based on the fact that for the general mesh, there are no stiffness coefficients in the global matrix which can be uniquely identified with an element which shares all its boundaries with other elements. A more interesting statement can be made: the disjoint element stiffness 496
Journal of The Franklin Institute
Critical Features of the Finite Element Method
Relation
TABLEI of Norms to Accuracy
Undeformed panel
Deformed panel
Basic element
4 Element ortlculatlon
Material
N.A.
Solution error
S pectro I norm
NA
4.n
37re55
rectangle (5) Hyperboltc rectangle (18) Biased triangles
(5)
Laminated triangles (19)
DO59
4.5
+ El
0,068
3.2
lzlk1
0. I30
3.0
0.136
2.0
matrix for the continuum element has no physical counterpart. This is so because no stiffness matrix with a finite number of degrees of freedom can include the infinite number of boundary displacement modes of the continuum. As a consequence, forces are usually required at all joints of a continuum when a unit displacement is imposed at one joint to prevent motions in all other degrees of freedom. Since this condition is the physical interpretation for coefficients in the global stiffness matrix, it follows that a disjoint continuum element has no physical realization. We conclude that some approximations must be introduced to derive a continuum finite element model from test data. The general approach is to use coupon tests to establish elastic constants and apply established methods of mathematics to develop coefficients of the element modeling equations. In perspective, the approach includes the following steps: (1) Assume displacement, and/or stress distribution, over the element domain as a function of nodal displacement parameters. Vol. 302, Nos. 5 & 6, November/December
1976
497
R. J. Melosh
(2) Introduce these functions into differential or integral equilibrium equations for the domain. (3) Choose a function by which error in satisfying equilibrium is to be weighed. (4) Make the error stationary with respect to the nodal displacement variables, thereby attaining a “best” estimate of the element Eqs. (1). Within this general approach, at least 768 analysis “methods” can be uniquely defined (17). Some of these are defined in detail in the other papers of this issue. The method selected for element development is important because it provides the basis for interpreting results. It establishes the types of discontinuities implied, the basis for determining analysis accuracy, and whether the solution is a general or bounding result. A question often asked by the analyst is which element model should be used for a specific structural continuum. The question must relate to efficiency since whatever approximations are introduced in steps 1-4, they must vanish as the element size approaches zero. With respect to elements of finite size, the question can only be answered by numerical experiments. As an incidental diversion, we observe that there is no best element for a given problem with a given mesh. Since disjoint element models are physically unrealizable, the most efficient analysis will use different element models in different regions of the system because the best assumed boundary displacement modes vary throughout the system. As in the case of networks, the accuracy of finite element analysis must depend only on the accuracy of measurements of material characteristics as element size gets small. Thereby, the fact that there is no physical counterpart to the finite element model is not a block to obtaining relevant response predictions by the approach. V. Finite Elements Analysis of Engineering Systems This section elaborates on the idea that the finite element method provides a process for analysis of engineering systems. The synthesis process accounts for systems which are bound together by continuity and balance conditions. Such systems include a wide variety of engineering systems. These conditions are expressed mathematically in terms of two types of variables corresponding to the deflection and stress variables of structural systems. These are “extrinsic” and “intrinsic” variables, respectively. Extrinsic variables, like deflection, are characterized by the fact that they are usually measurable by instrumentation external to the system. They are measured relative to some arbitrary datum. They are the natural variables in which to specify continuity requirements. Intrinsic variables, like stress, are usually not directly measureable in physical tests. They may be deductible from measurements of extrinsic variables. Their absolute rather than their relative value is of importance. They are the natural variables for equilibrium (balance) equations. Table II defines the extrinsic and intrinsic variables for a variety of engineering systems. The set of variables for a particular system is not unique as 498
Journal of The
Franklin
Institute
Critical Features of the Finite Element Method TABLE II Natural Variablesfor Engineering Systems
System
Intrinsic variable
Extrinsic variable
Fluids
Deflection strain Head
Stress internal force Flow rate
Heat transfer
Temperature
Flow rate
Electrical
Voltage
Current
Magnetic
Magnetomotive force
Flux
Structures
Material characteristic Stiffness damping Viscisity compressibility Conductivity capacitance Resistance inductance capacitance Reluctance
suggested by entries for the structural system. Other special models could be added to this list. When the system has a potential (plane stress for structures; irrotational and invisicid flow for fluids), the potential is the extrinsic variable and its derivatives, the intrinsic. References (21)-(26) illustrate the application of the finite element method to some of these technologies. Thus, the full disciplinary scope of finite element analysis, as represented in Table II, is being realized. Most finite element analysis is based on use of the extrinsic variables as the independent variable. (Even when stress functions are used for finite element analysis (27), this is the case.) It has been shown that selection of the intrinsic variables as the independent variables can result in four times as many calculations as by using the extrinsic (28). VI. Conclusion The features of the finite element method common to all papers on the method are: (1) The synthesis process, which includes treatment of both continuity and balance conditions. (2) Finite element models which identify with disjoint regions of the system and relate the extrinsic and intrinsic variables to provide the sole characterization of the behavior of the material of the system. Condition (1) allows the approach to span many engineering systems and may be taken as a definition of same. Condition (2) ultimately restricts accuracy of the method to the accuracy of measurements of material behavior. Study of’ the element models in the context of networks and continua indicate that it is the disjoint region constraint which evokes diversity of element models-a diversity which must disappear as the elements become vanishingly small. In particular: (a) For networks, a disjoint region representation is natural when extrinsic variables are used as independent variables. Vol.302,Nos.5 &6, November/December
1976
499
R. J. Melosh
Element stiffness coefficients are directly deductible by measuring extrinsics of the physical systems. (b) For continua, disjoint models are unnatural. Approximations, at least for behavior along element boundaries, are necessary to construct a disjoint model. The mathematical requirement of analysis is that whatever approximations are used, the element behavior must approach that of the material as the element size approaches zero. The principal limitation associated with components 1 and 2 above, is that finite element analysis cannot be applied to synergistic systems.
References (1) R. E. Templeton, “Introduction, Int. Conf. on Vehicle Structural Mech.“, SAE, Detroit, Michigan, March, 1974. (2) R. Courant, “Variational methods for the solution of problems of equilibrium and vibrations”, Bull. Am. Math. Sot., Vol. 49, pp. l-23, 1943. (3) W. Prager and J. L. Synge, “Approximations in elasticity based on the concept of function spaces”, Qt. Appl. Math., Vol. 5, pp. 211-269, 1947. (4) J. L. Synge, “The Hypercircle Method in Mathematical Physics”, Cambridge Press, Cambridge, pp. 168-213, 1957. (5) M. J. Turner, R. W. Clough, H. C. Martin, and L. P. Topp, “Stiffness and deflection analysis of complex structures”, J. Aeron. Sci., Vol. 23, pp. 805-823, 1956. (6) MARC-CDC, “Nonlinear Finite Element Analysis Program”, User’s Manual, Vol. 1, Control Data Corp., 1975. (7) H. Armen, H. Levine, A. Pifko, and A. Levy, “Nonlinear Analysis of Structures”, Grumman Research Dept. Rep., RE-454, May, 1973. (8) H. A. Balmer, J. St. Doltsinis, and M. Koenig, “Elasto-plastic and creep analysis with the ASKA program system”, Comput. Methods Appl. Mech. & Engng, Vol. 3, pp. 178-194, 1974. (9) K. J. Bathe and E. L. Wilson, “NONSAP-A general finite element program for nonlinear dynamic analysis of complex structures”, Proc. 2nd Int. Conf. on Structural Mechanics in Reactor Technology, Papo M3/1. Berlin, Sept. 1973. (10) R. H. Mallett and L. A. Schmit, “Nonlinear structural analysis by energy search”, .I. Struct. Div., ASCE, Vol. 93, No. ST3, pp. 221-234, 1967. (11)M. P. Kamat, R. J. Melosh, D. E. Killian, and G. W. Swift, “Theoretical Contractor’s Report, Palo Alto, Calif., 1976. (12) G. P. Bazeley, Y. K. Chueng, B. M. Irons, and 0. C. Zienkiewicz, “Triangular elements in plate bending+onforming and non-conforming solutions”, Proc. Conf. on Matrix Methods in Strut. Mech., pp. 547-576, AFFDL-TR-66-80, Wright-Patterson A.F.B., Ohio, 1966. (13) G. Strang and G. J. Fix, “An Analysis of the Finite Element Method”, Prentice Hall, Englewood Cliffs, New Jersey, 1973. (14) B. M. Irons, “The Semiloof Shell Element”, University of Calgary, Rep. TZNlN4, Alberta, Canada, 1975. (15) R. L. Taylor, P. J. Beresford, and E. L. Wilson, “A non-conforming element for stress analysis”, Int. J. Num. Methods in Engng., pp. 236-244, 1976. (16) S. N. Key, Private communication on the relation between finite element and constitutive models, 1964.
500
Journal of The Franklin
Institute
Critical Features of the Finite Element Method
(17)R. J. Melosh, “The optimum Struct., pp. 241-263,
approach Vol. 1, 1971.
to analysis
of elastic continua”,
Comput.
&
(18) R. J. Melosh, “Basis for deviation of matrices for the direct stiffness matrix”, AIAA J. pp. 1681-1637, 1963. (19) R. H. MacNeal, Ed., “The NASTRAN Theoretical Manual NASA SP-221(01), National Aeronautics and Space Administration, April, 1972. (20) R. J. Melosh and D. W. Lobitz, “On a numerical sufficiency test for monotonic convergence of finite element models”, AIAA J., Vol. 93, pp. 675-678, May, 1975. (21) V. C. Zienkiewicz, A. K. Bahiani and P. L. Arlett, “Three-dimensional fluid problems by finite elements”, Engineer Vol. 224, pp. 547-560, 1967. (22) J. T. Oden, “Finite element applications in mathematical physics”, Conf. on Finite Elements and Applications, Uxbridge, England, 1972. (23) M. A. Biot, “Thermoelasticity and irreversible thermodynamics”, J. AppZ. Phys., Vol. 27, pp. 240-253, 1956. (24) A. J. Baker, “Computational fluid mechanics using the finite element method”, Bell Aerospace Rep. 8500-601001, July, 1973. of finite elements (25) 0. C. Zienkiewicz, P. L. Arlett, and A. K. Bahrani, “Application J. IEEE, Vol. 115, pp. 1762-1766, 1968. to solution of Helmholz’ equations”, “Finite element formulation of general (26) J. T. Oden and B. E. Kelley, Int. J. Num. Methods in Engng., pp. 161electrothermo-elasticity problems”, 180, 1971. (27) B. M. Fraeijs De Veubeke, “Displacement and equilibrium models in the finite element method”, Stress Analysis, Chap. 9, pp. 184-216, Wiley, New York, 1965. (28) R. M. Bamford, “Force method vs the displacement method”, JPL Memo, 12, April, 1970.
Vol. 302, Nos. 5 & 6, November/December
1976
501
J. MELOSH
MARC
Analysis
Research
Corporation,
Palo Alto,
California
This paper examines the essential features of finite element analysis and identifies the inherent limitations of the method. It defines those aspects of finite element analysis which are common to all applications of the method. It discusses and compares the requirements and form of element models for computer simulation of structural networks and continua. It generalizes the terminology to clarify use of the method for a variety of physical systems. The perspective of the paper is that the finite element method is a discipline and problem invariant technique for synthesizing models of any engineering system; that it encompasses a wide variety of analysis objectives including developing exact, approximate, and bounding response predictions; that it can be fleshed into a variational, weighted residual, finite differenence, or a semi-empirical analysis method; and that its intrinsic limitations on fidelity rest on the relevance of physical measurements and constitutive laws hung on its framework.
ABSTRACT:
I.
Introduction
Templeton (1)observes that the finite element method is the most significant advance in engineering analysis since the 1920’s. This is at once a general recognition of the increased power and scope of engineering analysis credited to the method and a particular testimony to its practical value in application to the complex structural and heat transfer problems of the automobile industry. This dramatic tribute belies the few and simple ideas which are the skeleton of the method. Chief among these is the concept that a mathematical model of any engineering system can be synthesized from models of disjoint parts. The fact that this concept inbues an analysis modularity well-suited to digital computer implementation is of major importance to widespread use of the method. The consequence that the framework of finite element analysis is independent of engineering discipline, geometric and topological complexity, and response intricacy is the source of its power. The engineering, physics, and mathematical literature project the finite element method as a variational or weighted residual method. It expresses it in scalar, matrix, or tensor form involving differential, integral, or difference equations. It relates it to networks, continua and compound systems. It describes finite element analysis in the context of structures, elasticity, heat transfer, electricity, hydrodynamics, and/or magnetism. All of which may lead to some confusion about what the method is and what its limitations are. This paper seeks to identify the skeleton of finite element analysis and establish its inherent range of application. Most of the paper presents the method as a method of analysis of structural systems which behave as linear 489
R. J. Melosh
systems. It provides an analysis of the method from a phenomenological viewpoint: i.e. it evokes limitations of the method as reflected in the numerics which represent it. The presentation moves from a description of the finite element skeleton to study of its bones. Section II reviews the basic steps of analysis. Section III and IV consider analysis of networks and continua and the relation between these analysis classes. Section V replaces structural concepts with more abstract ones in identifying the interdisciplinary scope of the method. The last section briefly summarizes the key ideas of the paper. Though sparks of the finite element method can be found in the works of Courant (2), Prager (3), Synge (4) and others, the flames were kindled when Turner et al. (5) put the method in practical perspective. Since then, the ideas have undergone many transformations, but it is to the Turner work that this paper owes its largest debt. II. The Finite Element
Analysis Method
Regarding the finite element method as a means for simulating structures, this section reviews and discusses implications of its synthesis process. The building block of the analysis is the finite element model. The model pertains to a disjoint region of the system. It relates forces at nodes on the boundary of the element to corresponding nodal displacements, kii=U
(I)
where k is element stiffness matrix, a square matrix, ii is the vector of displacements of the nodes and U is the vector of nodal forces. The vectors ii and U are cogradient. Thus, if the axes to which they are referenced are rotated and translated, they are expressed in the new coordinate system by ij=Rii
Q=RU
where Q and Q are the vectors of joint displacements new axes. In the new basis, Eq. (1) is given by RkR-‘ij
(2) and forces related to the
= Q.
(3)
The potential energy, V, for the finite element model is defined as v=~liku-iiU=~~RkR-‘q-~Q
(4)
where the terms involving the 4 factor represent the strain energy and those involving joint forces represent external work. Thus, Eq. (1) defines a configuration state which is a stationary value of the potential. Equation (4) defines a scalar functional in which the element equations are embedded. Equations of two or more elements include equations of each element and constraint equations requiring that nodes common to the two elements have common values for corresponding components of displacement. Thus, if the 490
Journal
of The
Franklin
Institute
Critical Features of the Finite Element Method equations
for two elements
are represented
by
(5) the constraint
equations
can be expressed
by
0;! =
[mj,)
(6)
where C is the constraint matrix, iz is the vector of distinct displacement components as represented by & and &. Substituting Eqs. (5) and (6) in (4), and requiring that the function be stationary with respect to displacement changes yields the equations for the pair of elements in the form, C’KCiiz
= Qz
with
(Qz>= [C](!f).
Here, the QX are energy equivalent joint forces compliant with the continuity constraint. Introduction of constraints involves post- and pre-multiplication by the constraint matrix and its transpose as long as a potential exists. Some distinction between these transformations and coordinate transformations is lost if orthogonal coordinate systems are involved for then, R-r = R’.
(8)
Boundary conditions can be parsed in the form of Eq. (6). Some of the 41 will be constants if displacements are prescribed. If the values of some element joint forces are prescribed, coefficients of the constraint matrix can be abstracted from an element model equations, of Eqs. (l), thus defining a constraint. The equations that result from combining equations from all elements are the system load-deflection relations. These are solved for the unknown deflections. With all deflections known, the element equations (1) furnish the nodal forces on each element. The finite element synthesis process accommodates structural systems: systems which involve only equilibrium and displacement continuity conditions across disjoint elements. The behavior characteristics of the structure reside only in the element models and interelement nodal displacement constraints. The stiffness matrix for the complete system reflects the character of element matrices. With an appropriate choice of coordinate axes (e.g. orthogonal), both types of transformations are of contragredient form. Thus, if all element matrices are symmetric, the system matrix is. If all element matrices are positive semi-definite, the system stiffness matrix must also be positive semidefinite or positive-definite as a direct consequence of the definition of definiteness. The sparseness of the system stiffness matrix is well-known. In the usual case, the non-zero-coefficients in any column of the stiffness matrix can occur
Vol. 302, Nos. 5 & 6, November/December
1916
491
R. J. Melosh
only between a degree of freedom and its topologically closest neighbors. This sparseness is a direct consequence of the zero coefficients in the interelement constraint matrix of Eqs. (6). The requirement that finite elements represent disjoint substructures limits the variables included in Eqs. (6) to displacements for shared nodes and thus enforces sparseness of the system stiffness matrix. It is noted that Eqs. (1) through (8) are the basis for non-linear as well as linear analysis, despite the linear form of Eq. (1). These analyses involve tracking behavior through a sequence of loading steps each of which is assumed to involve a conservative system. In most formulations, each step involves an incremental loading and linearized element models (6)-(9). In some, nonlinearity is admitted (lo), (11). III. Analysis
of Networks
A network element is a finite element which represents uniquely and without mathematical approximation all admissible modes of behavior of the disjoint region it models. A network is a collection of network elements. A structural network is any physical system satisfactorily modeled as a network. Based on these definitions, the next paragraphs review the mathematrical requirements of network elements and the characteristics of network modeling. Figure 1 depicts some structural networks, networks, and network elements. These include rigid link-spring, torque tube, truss, and frame structures. These networks exhibit features common to all networks. Network elements involve one or two nodes. Any two nodes in the network can be connected by only one element. To distinquish networks from continua, network elements are called “branches” and nodes, “joints”. Since a branch must represent all structural characteristics without approximation, its mathematical model must be (a) based on displacement variables which can represent all possible boundary deformations, (b) include all elastic deformation states within the element, and (c) contain an explicit statement of branch equilibrium. Violation of any of these requirements disqualifies a finite element as a branch. For example, if only lateral displacements are used as nodal variables for a beam element, it cannot be a branch because the assumption of the beam equations that plane sections remain plane can not be enforced across element boundaries with this beam model. Satisfaction of requirements (a) and (b) results in load-deflection relations for the branch, corresponding to Eq. (l), which involve a symmetric positivedefinite stiffness matrix. Symmetry arises because the model is linear and conservative. It expresses the Maxwell-Betti theorem and is reflected in the self-adjoint differential equations of the beam element. Positive definiteness is a consequence of the fact that loads do work upon the structure in inducing deformation. Branch load-deflection relations, of the form of Eqs. (l), which involve only elastic deformations, can be modified to require satisfaction of macroscopic
492
Journal of The Franklin
Institute
Critical Features of the Finite Element Method
Mechanical
Drive
system
shafl
Tuss
Cd)
H Frame Structural
network
Network
element
Networks
FIG. 1. Structural networks, network, and network elements. equilibrium.
The equilibrium
requirements U,=CU,
(9)
where U, are undefined joint forces (reactions), augment the potential energy function. Minimization of this function with respect to the unknown displacements and the LaGrange multipliers then produces the equations:
(10) where the LaGrange multipliers have been interpreted as displacements in the reaction degrees of freedom. It is noted the linear dependence of the last r of Eqs. (10) signal a degeneracy of the element stiffness matrix reducing the rank by as many equations as the number of independent equilibrium requirements. The primary consequences of the branch characteristics are that the system stiffness matrix must be symmetric, positive semi-definite, and imply satisfac- . tion of equilibrium for each branch.. The symmetry is exploited to reduce the data processing task. With constraints introduced to eliminate system motions Vol. 302, Nos. 5 & 6, November/December
1976
493
R. J. Melosh
which involve no branch deformations, a positive definite stiffness matrix remains to insure a unique and minimum potential solution. Joint and branch equilibrium are often used as the basis for checks on combining of branch equations and solving of the system equations, despite unsatisfactory numerical characteristics of this test. Branch equations can be developed analytically or empirically. Analytical development involves defining a mathematical model and establishing its parameters by testing. For example, a differential equation can be used to model the beam branch. This model is particularized by measurements of the Young’s and Shear moduli for the material. To avoid mathematical approximation, exact solutions of the beam equation can be used in compliance with the definition of stiffness coefficients as joint forces associated with unit and zero joint displacements. Then the accuracy of the branch model depends on the accuracy of measurements of elastic constants. Network branch parameters can be uniquely particularized from physical tests on the structural system. Tests may involve measurements of displacements alone. This follows directly from the fact that since any two joints can be connected by a single branch, the coefficients of every non-diagonal partition of Eq. (10) can be uniquely defined by tests of the complete system. The remaining coefficients of Eqs. (10) can then be developed from the equilibrium requirements for each branch. For the frame of Fig. l(d), for example, displacements and rotations of each of the joints can be taken for a unit force at each joint and for a unit moment at each joint. Inversion of this flexibility matrix produces the system stiffness matrix. The off-diagonal joint partition for joints i and i (a 2 X 2 matrix) provides data for evaluating the three independent elastic coefficients needed for the load-deflection coefficients of Eqs. (1) for branch i-j. Similarly, equations for every other branch can be constructed without concern for the details of geometry and material composition. The scope of network analysis is limited by the validity of the form of the branch equations, Eqs. (1) and the accuracy of the measured displacements from which equation coefficients are deduced. Because branch equations cannot involve mathematical approximations, application is usually restricted to structures, like those of Fig. 1, whose topology articulates slender branches. IV. Analysis
of Continua
and Compound Systems
A continuum element is a finite element which approximates behavior of a disjoint part of a system. A continua is a collection of continuum elements. A structural continuum is a structural system which is satisfactorily modeled as a continuum. A compound system consists of both network and continuum elements. Using these definitions, the next paragraphs review mathematical requirements for continua and the characteristics of structural continua and compound systems. Figure 2 shows some structural continua, continuum elements, and a continua. These include membrane, plate, shell, and three-dimensional solid 494
Journal of The Franklin Institute
CriticalFeatures of the Finite Element Method
(0)
Membrane El---8
(b)
Plate
(cl
Shell
(d)
Solid
(a) Structural
continua
(b)Contlnuum
element
(c)Continua
FIG. 2. Structural continua, continuum elements, and continua. systems. Continua engage finite elements which usually have three or more nodes. A particular pair of nodes may be connected by one or more elements. As a minimum, the continuum element must satisfy the requirements for a branch as the element becomes vanishingly small. Thus, it must (a) involve boundary displacements which imply satisfaction of displacement continuity across element boundaries, (b) represent the stress-strain relations, and (c) imply satisfaction of microscopic equilibrium conditions everywhere in the element interior and across boundaries. Since the differential equations of elasticity can be cast in self-adjoint form, stiffness matrices must approach symmetry as the element size approaches zero. Since stress-strain relations of most structural material represent Green’s materials, stiffness matrices are usually positive semi-definite or positive definite. Conversely, elements of finite size displacements need not define all possible modes of boundary behavior, be continuous across boundaries, satisfy macroscopic equilibrium or imply satisfaction of compatibility within an element. Stiffness matrices need not be symmetric nor positive definite. An interesting aspect of finite element analysis of continua is the evolution of ad hoc criteria to insure the adequacy of element models. These criteria are phenomenological-they qualify the element matrix without regard for the means used in establishing the stiffness coefficients. Of prime importance is Iron’s patch test (12). This test is based on the proven hypothesis (13) that if deformations associated with a possible constant strain state are imposed on the boundaries of a set of the elements and the Vol. 302, Nos. 5 & 6, November/December
1976
495
R. J. Melosh
analysis yields the exact solution, then the element guarantees to lead to the solution of the equations of elasticity for any structure, when enough finite elements are used. In principle, this test must be used for all admissible constant states (e.g. x-strain, y-strain, and shear strain for a membrane), all element shapes, and all material models. This criteria brings respectability to many of the finite element models developed early in the evolution finite element technology. Irons (14) shows that with this criteria in mind, one can generate efficient element models (high accuracy/element) without concern for more sophisticated criteria. Taylor et al. (15)illustrate addition of behavior states to an element with regard for this criteria, to improve its efficiency. The patch test is closely related to an observation by Key (16) that adequacy of an element is related to its ability to reflect the modes of behavior associated with its matrix of elastic constants. For example, the material matrix for a membrane involves dilatation, stretching along one axis and contracting along another, and shear. A linear combination of those modes transforms them to constant strain states. The data of Table 1 illustrates that a stronger statement can be made-the validity and efficiency of an element model depends on the fidelity with which it represents the modes of the material matrix when the mesh is sufficiently fine. The uniformly thick membrane shown in Table 1 is articulated by 35 rectangular elements. Analyses for displacements were made for the four elements defined-each analysis involving 64 equilibrium equations. The mean error in displacements was determined by comparing solution results with an analysis with 1050 equations. The ratio of the maximum to minimum eigenvalue (the spectral norm) was evaluated as a scalar measure of the element matrices involved. The elements are listed in accordance with the closeness of the norm to the material matrix norm. The fact that this ordering corresponds exactly with the solution error suggests that when many elements are used, the element model whose natural behavior (as reflected by its eigenvalues and vectors) most closely corresponds with microscopic material natural modes will be most efficient. Recently a numerical sufficiency test has been proposed to determine whether an element model, given in the form of its stiffness matrix, must involve monotonic convergence (20). This test also checks that element models satisfy macroscopic equilibrium requirements. Thus, in combination with the patch test, there exists at least one set of numerical tests for qualifying a continuum element stiffness matrix. Can non-destructive tests of a structural continuum be used to evaluate coefficients of an acceptable continuum element stiffness matrix? A simple “no” answer does not do justice to the question. The no answer can be based on the fact that for the general mesh, there are no stiffness coefficients in the global matrix which can be uniquely identified with an element which shares all its boundaries with other elements. A more interesting statement can be made: the disjoint element stiffness 496
Journal of The Franklin Institute
Critical Features of the Finite Element Method
Relation
TABLEI of Norms to Accuracy
Undeformed panel
Deformed panel
Basic element
4 Element ortlculatlon
Material
N.A.
Solution error
S pectro I norm
NA
4.n
37re55
rectangle (5) Hyperboltc rectangle (18) Biased triangles
(5)
Laminated triangles (19)
DO59
4.5
+ El
0,068
3.2
lzlk1
0. I30
3.0
0.136
2.0
matrix for the continuum element has no physical counterpart. This is so because no stiffness matrix with a finite number of degrees of freedom can include the infinite number of boundary displacement modes of the continuum. As a consequence, forces are usually required at all joints of a continuum when a unit displacement is imposed at one joint to prevent motions in all other degrees of freedom. Since this condition is the physical interpretation for coefficients in the global stiffness matrix, it follows that a disjoint continuum element has no physical realization. We conclude that some approximations must be introduced to derive a continuum finite element model from test data. The general approach is to use coupon tests to establish elastic constants and apply established methods of mathematics to develop coefficients of the element modeling equations. In perspective, the approach includes the following steps: (1) Assume displacement, and/or stress distribution, over the element domain as a function of nodal displacement parameters. Vol. 302, Nos. 5 & 6, November/December
1976
497
R. J. Melosh
(2) Introduce these functions into differential or integral equilibrium equations for the domain. (3) Choose a function by which error in satisfying equilibrium is to be weighed. (4) Make the error stationary with respect to the nodal displacement variables, thereby attaining a “best” estimate of the element Eqs. (1). Within this general approach, at least 768 analysis “methods” can be uniquely defined (17). Some of these are defined in detail in the other papers of this issue. The method selected for element development is important because it provides the basis for interpreting results. It establishes the types of discontinuities implied, the basis for determining analysis accuracy, and whether the solution is a general or bounding result. A question often asked by the analyst is which element model should be used for a specific structural continuum. The question must relate to efficiency since whatever approximations are introduced in steps 1-4, they must vanish as the element size approaches zero. With respect to elements of finite size, the question can only be answered by numerical experiments. As an incidental diversion, we observe that there is no best element for a given problem with a given mesh. Since disjoint element models are physically unrealizable, the most efficient analysis will use different element models in different regions of the system because the best assumed boundary displacement modes vary throughout the system. As in the case of networks, the accuracy of finite element analysis must depend only on the accuracy of measurements of material characteristics as element size gets small. Thereby, the fact that there is no physical counterpart to the finite element model is not a block to obtaining relevant response predictions by the approach. V. Finite Elements Analysis of Engineering Systems This section elaborates on the idea that the finite element method provides a process for analysis of engineering systems. The synthesis process accounts for systems which are bound together by continuity and balance conditions. Such systems include a wide variety of engineering systems. These conditions are expressed mathematically in terms of two types of variables corresponding to the deflection and stress variables of structural systems. These are “extrinsic” and “intrinsic” variables, respectively. Extrinsic variables, like deflection, are characterized by the fact that they are usually measurable by instrumentation external to the system. They are measured relative to some arbitrary datum. They are the natural variables in which to specify continuity requirements. Intrinsic variables, like stress, are usually not directly measureable in physical tests. They may be deductible from measurements of extrinsic variables. Their absolute rather than their relative value is of importance. They are the natural variables for equilibrium (balance) equations. Table II defines the extrinsic and intrinsic variables for a variety of engineering systems. The set of variables for a particular system is not unique as 498
Journal of The
Franklin
Institute
Critical Features of the Finite Element Method TABLE II Natural Variablesfor Engineering Systems
System
Intrinsic variable
Extrinsic variable
Fluids
Deflection strain Head
Stress internal force Flow rate
Heat transfer
Temperature
Flow rate
Electrical
Voltage
Current
Magnetic
Magnetomotive force
Flux
Structures
Material characteristic Stiffness damping Viscisity compressibility Conductivity capacitance Resistance inductance capacitance Reluctance
suggested by entries for the structural system. Other special models could be added to this list. When the system has a potential (plane stress for structures; irrotational and invisicid flow for fluids), the potential is the extrinsic variable and its derivatives, the intrinsic. References (21)-(26) illustrate the application of the finite element method to some of these technologies. Thus, the full disciplinary scope of finite element analysis, as represented in Table II, is being realized. Most finite element analysis is based on use of the extrinsic variables as the independent variable. (Even when stress functions are used for finite element analysis (27), this is the case.) It has been shown that selection of the intrinsic variables as the independent variables can result in four times as many calculations as by using the extrinsic (28). VI. Conclusion The features of the finite element method common to all papers on the method are: (1) The synthesis process, which includes treatment of both continuity and balance conditions. (2) Finite element models which identify with disjoint regions of the system and relate the extrinsic and intrinsic variables to provide the sole characterization of the behavior of the material of the system. Condition (1) allows the approach to span many engineering systems and may be taken as a definition of same. Condition (2) ultimately restricts accuracy of the method to the accuracy of measurements of material behavior. Study of’ the element models in the context of networks and continua indicate that it is the disjoint region constraint which evokes diversity of element models-a diversity which must disappear as the elements become vanishingly small. In particular: (a) For networks, a disjoint region representation is natural when extrinsic variables are used as independent variables. Vol.302,Nos.5 &6, November/December
1976
499
R. J. Melosh
Element stiffness coefficients are directly deductible by measuring extrinsics of the physical systems. (b) For continua, disjoint models are unnatural. Approximations, at least for behavior along element boundaries, are necessary to construct a disjoint model. The mathematical requirement of analysis is that whatever approximations are used, the element behavior must approach that of the material as the element size approaches zero. The principal limitation associated with components 1 and 2 above, is that finite element analysis cannot be applied to synergistic systems.
References (1) R. E. Templeton, “Introduction, Int. Conf. on Vehicle Structural Mech.“, SAE, Detroit, Michigan, March, 1974. (2) R. Courant, “Variational methods for the solution of problems of equilibrium and vibrations”, Bull. Am. Math. Sot., Vol. 49, pp. l-23, 1943. (3) W. Prager and J. L. Synge, “Approximations in elasticity based on the concept of function spaces”, Qt. Appl. Math., Vol. 5, pp. 211-269, 1947. (4) J. L. Synge, “The Hypercircle Method in Mathematical Physics”, Cambridge Press, Cambridge, pp. 168-213, 1957. (5) M. J. Turner, R. W. Clough, H. C. Martin, and L. P. Topp, “Stiffness and deflection analysis of complex structures”, J. Aeron. Sci., Vol. 23, pp. 805-823, 1956. (6) MARC-CDC, “Nonlinear Finite Element Analysis Program”, User’s Manual, Vol. 1, Control Data Corp., 1975. (7) H. Armen, H. Levine, A. Pifko, and A. Levy, “Nonlinear Analysis of Structures”, Grumman Research Dept. Rep., RE-454, May, 1973. (8) H. A. Balmer, J. St. Doltsinis, and M. Koenig, “Elasto-plastic and creep analysis with the ASKA program system”, Comput. Methods Appl. Mech. & Engng, Vol. 3, pp. 178-194, 1974. (9) K. J. Bathe and E. L. Wilson, “NONSAP-A general finite element program for nonlinear dynamic analysis of complex structures”, Proc. 2nd Int. Conf. on Structural Mechanics in Reactor Technology, Papo M3/1. Berlin, Sept. 1973. (10) R. H. Mallett and L. A. Schmit, “Nonlinear structural analysis by energy search”, .I. Struct. Div., ASCE, Vol. 93, No. ST3, pp. 221-234, 1967. (11)M. P. Kamat, R. J. Melosh, D. E. Killian, and G. W. Swift, “Theoretical Contractor’s Report, Palo Alto, Calif., 1976. (12) G. P. Bazeley, Y. K. Chueng, B. M. Irons, and 0. C. Zienkiewicz, “Triangular elements in plate bending+onforming and non-conforming solutions”, Proc. Conf. on Matrix Methods in Strut. Mech., pp. 547-576, AFFDL-TR-66-80, Wright-Patterson A.F.B., Ohio, 1966. (13) G. Strang and G. J. Fix, “An Analysis of the Finite Element Method”, Prentice Hall, Englewood Cliffs, New Jersey, 1973. (14) B. M. Irons, “The Semiloof Shell Element”, University of Calgary, Rep. TZNlN4, Alberta, Canada, 1975. (15) R. L. Taylor, P. J. Beresford, and E. L. Wilson, “A non-conforming element for stress analysis”, Int. J. Num. Methods in Engng., pp. 236-244, 1976. (16) S. N. Key, Private communication on the relation between finite element and constitutive models, 1964.
500
Journal of The Franklin
Institute
Critical Features of the Finite Element Method
(17)R. J. Melosh, “The optimum Struct., pp. 241-263,
approach Vol. 1, 1971.
to analysis
of elastic continua”,
Comput.
&
(18) R. J. Melosh, “Basis for deviation of matrices for the direct stiffness matrix”, AIAA J. pp. 1681-1637, 1963. (19) R. H. MacNeal, Ed., “The NASTRAN Theoretical Manual NASA SP-221(01), National Aeronautics and Space Administration, April, 1972. (20) R. J. Melosh and D. W. Lobitz, “On a numerical sufficiency test for monotonic convergence of finite element models”, AIAA J., Vol. 93, pp. 675-678, May, 1975. (21) V. C. Zienkiewicz, A. K. Bahiani and P. L. Arlett, “Three-dimensional fluid problems by finite elements”, Engineer Vol. 224, pp. 547-560, 1967. (22) J. T. Oden, “Finite element applications in mathematical physics”, Conf. on Finite Elements and Applications, Uxbridge, England, 1972. (23) M. A. Biot, “Thermoelasticity and irreversible thermodynamics”, J. AppZ. Phys., Vol. 27, pp. 240-253, 1956. (24) A. J. Baker, “Computational fluid mechanics using the finite element method”, Bell Aerospace Rep. 8500-601001, July, 1973. of finite elements (25) 0. C. Zienkiewicz, P. L. Arlett, and A. K. Bahrani, “Application J. IEEE, Vol. 115, pp. 1762-1766, 1968. to solution of Helmholz’ equations”, “Finite element formulation of general (26) J. T. Oden and B. E. Kelley, Int. J. Num. Methods in Engng., pp. 161electrothermo-elasticity problems”, 180, 1971. (27) B. M. Fraeijs De Veubeke, “Displacement and equilibrium models in the finite element method”, Stress Analysis, Chap. 9, pp. 184-216, Wiley, New York, 1965. (28) R. M. Bamford, “Force method vs the displacement method”, JPL Memo, 12, April, 1970.
Vol. 302, Nos. 5 & 6, November/December
1976
501