- Email: [email protected]
Inherited error in finite element analyses of structures
Compvrrrs & Slructures, Vo1.3, pp. 1205 -12 17. Pergamon Pms
1973. Printed in Great Britain
INHERITED ERROR IN FINITE ELEMENT ANALYSES OF STRUCTURES-f ROBERT J. MELOSH 3939 Fabian Way, Palo Alto, California, 94303, U.S.A. A~~~-Mathematicians are quick to point out that round-off and truncation errors induced by the digital computer are only half of the manipulation errors in numerical analyses. The other half am the errors in quantizing the mathematical problem for computer solution: errors inherited for the equation solving process. This paper examines the relevance of inherited errors in structural analyses using the finite element concept and the digital computer. It illustrates error magnitudes using numerical experiments of simple structures. It constructs a theory explaining the errors for these systems. Having identified the most significant controllable computer error source, it describes a process for minimizing its contribution to the inherited error. The paper concludes that orders of magnitude between errors reported by various investigators can be explained by differences in inherited error. The most significant effect of these errors can be identified with inconsistency in problem formuiation. These inconsistencies can be eliminated by exploiting the existence of rigid body states in the finite eIement models. Thereby, solution errors introduced by inherited errors can be reduced to intrinsic errors in parameters defining the geometry, material characteristics, and boundary conditions.
1. INTRODUCTION How CAN a digital computer develop inaccurate data when relatively few calculations are involved? This can occur in the well known GIG0 (garbage in, garbage out) and infamous PICA (precision-in, critical arithmetic) modes. In structural analysis, it is convenient to divide the analysis at the point where the coefficients of the equations of motion are available. The input errors in the GIG0 mode are ‘inherited’ errors to the equation solving calculations. PICA errors, though they may arise in coefficient generation, more often occur in the equation solving calculations. In critical arithmetic, large relative errors are induced in few calculations because the information of importance is embedded in the lower bits of the floating point mantissa. In structural analyses, this can occur in solving the equations of motion when the equation sequence is bad, [1] or, in displacement method analysis, when one or more structural elements is nearly rigid. Though these conditions can induce large errors, these sources of error are now known. Modern structural analysis codes provide safeguards to sense these conditions and take appropriate action. The importance of inherited errors has long been recognized. Rosanoff et af. [2] in 1968 posed a simple problem, involving but two simultaneous equations, for which a 27 bit mantissa is insufficient to obtain two digit solution accuracy. t Presented at the National Symposium on Computerized Structural Analysis and Design at the School of Engineering and Applied Science, George Washington University, Washington, D.C., 27-29 March 1972. 1205
ROBERTJ. MEJBSH
1206
Figure 1 shows the magnitude of inherited errors in a set of structural problems. It illustrates inherited and manipulation errors in computer analysis of a cantilevered beam. Inherited error estimates are predicated on change in the last binary place of the mantissa [3]. These curves show that even though intrusion occurs only into the last binary place, solution error is comparable to the error induced by all the subsequent arithmetic truncation errors in the solving calculations.
I
0001
0
100
Numbers
FIG. 1. Inherited
I
1
200
300
400
of equations,
and manipulation
500
600
N=2J
error growth for beams.
This paper uses the case method to study inherited error rather than a bounding or probabilistic approach. The bounding approach leads to formulas defining the absolute worst case error. These formulas are often costly to evaluate and very conservative compared with expected errors [4]. The probabilistic approach, on the other hand, is not satisfactory because it only indirectly reflects how bad errors can actually become. The case method uses a worst case, but realistic problem, as the basis for evaluating error characteristics and magnitudes. The problem found worst for PICA mode errors is the cantilever beam [5, 61. Because this structure permits systematic accumulation of error in successive steps of the solution process, it will be assumed to be a worst case system for this study. Choice of the cantilevered beam problem as the case also has secondary advantages. It facilitates an exact analysis of error. It permits comparison of the effects of inherited errors with other effects. The problem of predicting response by the displacement method under a static loading will be considered. In this case, equation solving involves developing the solution of a set of linear, simultaneous equations. Since this is the central mathematical problem in structural analysis; be it linear, nonlinear, time dependent or not; the conclusions should have broad application.
Inherited Error in Finite Element Analyses of Structures
1207
For structural analysis, this paper reviews GIG0 mode errors and their effect. Of primary concern are errors that can be automatically reduced or compensated. The next section describes the relative importance of various sources of inherited error. The third section develops simple error estimation formulas. The fourth discusses ~ni~zing the effects of one source of inherited error. The last section summarizes report conclusions. 2. INHERITED
ERRORS
Of the sources of equation solving inherited error, discretization has justifiably received the most attention. Errors in generating equation coefficients, however, also have a significant effect on the accuracy of response predictions, and should be easier to minimize. This section leads to these conclusions by examining the sources and relative ma~itudes of inherited errors. Error sources
Idealization, numerical modeling, and coefficient generation are the principle sources of errors. Idealization errors are incurred in selecting the mathematical model for the system. Numerical modeling errors are implied in particularizing parameters of the mathematical model so it represents the system of interest. Generation errors arise in calculating the c~~cients in the load deflection equations from numerical input data. In finite element analysis of structures, idealization errors are evoked by selecting the distribution of ‘joints’ and ‘cuts’, electing finite element models to be used, choosing material stress-strain relations and failure criteria and selecting joint force and displacement boundary condition equations. Since implementation of the analysis is almost always directed by a computer code, many of these decisions are implicit in the selection of the computer code itself. Idealization errors are attributable to insufficient detail in the equations used to represent the structure. Lack of detail may be due to insu&ient discretization of the structure or omission of terms in the equations. Insufficient discretization occurs, for example, when the uniform load on a non-prismatic beam represented by a cubic displacement function is portrayed by energy equivalent joint forces and moments. An illustration of errordue to omission of terms is the error due to neglecting shear deformations in analyzing a beam. In a broad sense, idealization errors are errors involved in defining the abstract mathematical equations into which the physical system will be mapped. In this sense then, numerical modelIing errors are the errors in defining those parametersin the equations which are to be the knowns. Thus, in finite element analysis these errors are associated with the definition of geometric parameters (joint locations, element thicknesses, joint eccentricities), material parameters (elastic constants, coefficients of thermal expansion, density, strength limits) and boundary condition parameters (magnitudes of prescribed generalized displacements, load distribution coefficients, joint partial restraint values). Numerical modeling errors are introduced by humans. They may involve errors in judgment on the part of the analyst (assignment errors) or errors in communicating the numbers to the computer (tmnscription or keypunching errors). Generation errors may be caused by manipulation or programming errors. Manipulation errors include round-off and truncation errors evoked by limited precision arithmetic. Programming errors include errors in program logic and formula evaluation. In computer impIemented analysis of structures, it is advantageous to define generation errors as all errors induced from the operation of reading input data until the coefficients inherited
RCXERT J.
I 208
MELOSH
of the load-deflection equation are expressed in a floating point representation. Thus in the displacement method these errors include those induced in card reading, input conversion, scaling, nondimensionalizing, coordinate transformations, and evaluation of formulas for stiffness, loading, and stress coefficients. Of the errors discussed, only idealization errors of the discretization type and generation errors of the manipulative type are candidates for completely automatic error reduction or compensation. Though some automatic correction is possible for the other types of errors, those types will not be considered further here. Error magnitudes
The relative importance of discretization and generation manipulation errors is illustrated by the curves of Fig. 2. This figure shows discretization is the more important source of error for coarse meshes. Inherited error is potentially more important with finer meshes. Whereas discretization errors are bounded and will decrease with increasing mesh refinement, generation manipulation errors can induce unbounded errors in the solution and can monotonically increase with mesh refinement.
;, L z 0) z z G n
005 -
/
(diagonals
00055
600
200
Number
FIG. 2.
only)
of equations,
800
2J
Inherited and discretization error growth for beams.
The discretization error leading to the deflection errors of Fig. 2 are due to approximating the loading with a lumped parameter rather than an energy equivalent loading. The inherited error includes only the truncation errors induced in summing element stiffness matrices. A sound basis for automatic control of discretization error has been laid down by Carrel and Barker [7]. Their important work provides the keystone for increasing the understanding of how to distribute mesh points and what refinement to use. Carroll, Barker. McNeice and Marcal [8] have already illustrated some applications in structural analyses. Thus this paper can concentrate on the generation errors and their effects.
Inherited Error in Finite Element Analyses of Structures
1209
The solution errors of the lower curve of Fig. 2 are by no means the maximum generation induced errors that can occur in cantilevered beam analysis. They reflect only the effects of errors on the diagonal of the stiffness matrix. Figure 3 illustrates the error in deflection for a more representative generation error. This figure shows the error in deflection predictions for the beam as a function of the scale factor used for the analysis for a beam represented by 500 finite elements. The element stiffness matrices were divided by the abscissa values to simulate differences induced by using different dimensional units in defining the problem. The smallest errors provide a measure of the manipulation errors in the solving calculations. The largest errors, which are more than 103 x the smallest, reflect the effect of this inherited error on solution error. The dashed curve of Fig. 2 shows how the error for the factor of 11 varies with the number of equations representing the beam. The form of the curves for other scale factors is the same. Comparing this curve with the lower curve of Fig. 2 leads to the conclusion that this generation error source involves solution errors that are about 102 greater than those due to summation errors. Comparing the errors of the dashed curve of Fig. 2 with the manipulation errors in Fig. 1 shows that these inherited errors induce solution errors that are three times greater than those of manipulation in solving the equations. These data expose the importance of exercising care in performing numerical experiments. The variance in error may provide a basis for explaining the divergent conclusions published about error magnitudes. They illustrate that inherited errors may induce solution errors that are much greater than those due to manipulation errors introduced in solving the equations. In seeking an explanation for the variation of these solution errors with scale factor for beams, a basis for minimizing this type of error for more general finite elements will become evident.
2J = 1000 p = 27 bits R
0
.o + 0 Q
.N ;j E
b E
Element stiffness matrix divisor FIG. 3. Error variation with scale factor.
1210
ROBERT J. MELOSH
3. ERROR ESTIMATION To explain the results in Fig. 3, this section examines the difference equations for the beam and energy relations for the finite element. With this basis, the variation in solution errors is attributed to the inaccuracy in representing rigid body modes introduced by inherited errors. Dlyerence equation solution
The load-deflection relation for the beam finite element neglecting shear deformation is expressed by, 2aEf 3-a
Pl
3aG -6C 3aG
where
2a2F
_
-3aG a2F
ml
66 -3aG
.
(1)
p2
2a21:
m2
is the scalar multiplier is Young’s Modulus, is the bending moment of inertia, r is the element axial length, a wt,Pr are the displacement and force acting normal to the neutral axis at the joint at end 1 and @2sm2 are the angular rotation and moment associated with bending of the joint at end 2 of the element.
;
and
F=l+e,. G=l+e,
(2)
with e, and eG the relative errors introduced by scaling. There are only two distinct relative errors because factors of two among coefficients introduce no change in relative error using a binary floating point number representation. Forming the load-deflection relations for the complete prismatic and uniformly loaded system results in equilibrium equations for the beam, --6w,_ rG-3@,_ raG+12w,G--6w,,
tG+3&
laG=P
(3)
3~,_~G+30,_raF+4B,aF--3w,+tG+t3,,taF=M where P is the generalized normal force at each joint, P=pa3/aBl, and M is the corresponding generalized moment (M=ma2/aEI). With the equations in this form it is seen that the effect of a is merely to scale the loading coefficients. Hence errors in this solution, due to a, are inversely proposal to CL For the problems of interest, the boundary conditions are: w(r=O)=O,
e(r=o)=o,
(4)
1211
Inherited Error in Finite Element Analysesof Structures 3w,_ 1G-3wNG+8N_ -6w,_
1G+6~NG-38N_
laF+20+F=A4 laG-3BNaG=P.
If F= G= a = l., the difference equations are exactly represented. equations (3 and 4) for tip deflection is given by
w,,,=--&3N4+12N3-1N2-2N)
The sol$ion
of
(5)
when m =p = a = 1.0 and N is the number of finite elements. Suppose the solution of equations (3 and 4) are expanded in a power series in /3 where /l=6(F-G)/F
(6)
Then, subtracting equation (5) from the fist terms of this expansion and dividing by equation (5) yields a measure of the relative error of the solution, i.e. (7) where the convergence of the expansion is guaranteed when N2/?< 1. Equation (7) suggests that if NzB>O, the error is bounded. If NzpcO, the error IS unbounded. (The complete expansion will confirm this since eT1p can then be expressed, for a given N, as a constant and the cos NJ/? for 8~0 and a constant and cash N,/p for B>O.) For /INz% 1, and /? ~0; equation (7) shows that the tip deflection error increases approximately as the square of the number of finite elements. For /?N2< 1, and b
The solution of the difference equations furnishes the basis for explaining the erratic results depicted in Fig. 3. For this series of problems, the changing of scale factors results in significant changes in the relative importance of manipulation and inherited error effects. The cases considered can be divided into two classes: those which involve negligible equilibrium imbalance as measured by /3, and those where this imbalance is significant. Each case can be classified by calculating a-l and /? for the relevant divisor. These calculations yield the data given in Table 1. For this set of divisors, the equilibrium imbalance is zero. These data were developed on the same computer as used to generate the data of Fig. 3. Table 1 also lists the solution error for these cases.
1212
ROE~T 3. MEL~SH TABLE1. Errorst due to misrepresenting elasticity: Divisor l/3 115 l/7 1,2,4,8,16 7,14 15
a Error
Solution error:
0 0 0 0 7.5-7 3.7-7
622 -6 629 -6 6.27 -6 &Sg-6 6.17 -6 6.28 -6
t Relative errors from experiments with a 27 bit mantissa and truncation arithmetic. $ All errors involve numbers which are less than the exact values.
These data show that for these cases, solution errors are the result of only manipulation errors in the solution process. Solution errors are essentially unaffected by changing c1 errors from zero to nonzero. Since the relative error in the solution is three orders of magnitude greater than the CLerror, the error in misrepresenting stiffness for these cases is due to truncation errors in decomposition, forward, and back substitution. Thus the solution errors for all the small error cases, the circled points on Fig. 3, are due to manipulation errors. Data for the cases in which the equilibrium imbalance is not zero are listed in Table 2. Because the CIerrors for these cases make a negli~ble contribution these c1errors have been omitted. For convenience, the divisors incurring negative and positive /I errors have been grouped separately. Figure 4 provides a log-log plot of the solution errors as a function of the number of equations. Curves are included for four of the sets of divisors listed in Table 2. The
Number
of
equations,
N
FIQ. 4. Solution error due to equilibrium imbalance.
Inherited Error in Finite Element Analyses of Structures
1213
divisors with negative errors are shown by continuous curves; those for 3, 6, and 12 by a dashed curve. The relative manipulation error has been subtracted from the solution error.
Thus results in Fig. 4 pertain to only equilibrium errors. The continuous curves demonstrate that the negative #I cases have solution errors that increase as N2. The errors increase monotonically with 8, as prescribed by the analytical results. The positive /I case shows that as N increases, the error increases monotonically. However, the errors of sign opposite to the manipulation do not necessarily compensate for truncation, probably because they arise in a different way. This is illustrated by examining the data in Table 2. TABLE2. Errorsg due to equilibrium imbalance Divisor
/I Error
17 9 5.10 11 3,6.12 13
-158-S -1.67-8 -3.71-8 -6.14-s 2.23 -7 4.84 -8
Solution Error 1.48 -5 152-s 7.21-3 1.19 -2 150-s 6.18-6
5 Relative errors from experiments with a 27 bit mantissa and truncation arithmetic. 11All errors involve numbers which are less than the exact values.
Thus the large errors in Fig. 3 are caused by ‘negative’ equilibrium imbalances. The small ones involve manipulation errors or ‘negative’ interaction of manipulation and positive equilibrium imbalance. Using the analytical results, Table 3 provides the data for comparing the effects on beam solution errors of rounding, rather than truncating, in arithmetic. This table cites the relative error in satisfying equilibrium 8, for the cases of interest. Truncation errors are repeated for convenience. As would be expected, these data support the conclusion that maximum errors will be less with rounding than with truncation. Because the equilibrium errors often vary significantly with the arithmetic mode used, the arithmetic mode should not be overlooked. This emphasizes again that casual selection of computer tests can add to the confusion rather than increase understanding. TABLE3. Truncation and rounding equilibrium errorst Divisor
Truncation B error
1,2,4,8,16 3,6,12 5,lO 7,14 9 11 13 15 17 t Relative errors with a 27 bit mantissa.
Rounding jI error
0
0
2.23 -8 -3.71-s 0 -1.67-s -6.14-s 4.84 -8 0 -1.58-s
2.23 -8 3.71-s -2.61-s -1.67-s 1.96-s -2.42-s 0 6.18-11
1214
ROBERTJ. MELQSH
Eigenvalue analysis
The implications of inherited errors may be better understood by considering their effect on the eigenvalues of the stiffness matrix. Since a is a matrix scalar multiplier, it scales all eigenvalues without changing eigenmodes. Independent errors to coefficients of the ma’trix change both eigenvalues and modes. Assuming truncation errors, small errors to coefficients have little relative effect on zero terms. On a p bit computer with inherited error in the pth bit the elastic modes are effected in the pth bit. Similarly, only the pth bit of each elastic mode eigenvalue is reduced. The zero eigenvalues associated with rigid body modes become non-zero. They may be either small positive or negative value depending on the magnitudes of relative error to each coefficient. Thus, instead of being positive semi-definite, the element stiffness matrix can become indefinite. It is noted that if the forces and moments of each column of the stiffness matrix satisfy macroscopic equilibrium requirements for the element, this is sufficient to insure zero frequency rigid modes. This is not sufficient to avoid errors in the elastic eigenvalues and modes. The energy error due to the inherited error in the element matrix is proportional to the eigenvalues. Thus, expressing the ratio of translational and rotational eigenvalues to that of the lowest elastic mode provides an error indicator. Suppose that each of the coefficients may have a different inherited error with a relative value of eU where the subscripts relate to the row and column numbers of the matrix of equation (1). Then the indicators are:
ell--6e12+6e13-
6el,+4e22-6e2,+4e,,+3e,J
- 6e,, + 4e,,jqL qd
where IT is a relative measure of rigid translation error, ZR is a relative measure of rigid rotation error, qT and q, are rigid mode partrcipation factors, and qxr is the participation factor for the more flexible elastic mode. Note that Zr=ZR =0 if matrix coefficients satisfy equilibrium. These formulas confirm that when the rigid motion of the element is large compared with its elastic deformation, small errors in equilibrium can imply significant energy in ‘rigid’ modes. The indicator for translation projects large errors when very short elements are involved. This situation can be illustrated by the cantilevered beam solution of the test problem with 500 finite elements. For the scaling error rigid translation is exactly represented. Rigid rotation errors, however, would be expected to be particularly large for the tip section. An energy analysis shows the rigid mode energy in the tip element is 0.62 x 10-Z per cent of the total external work even though the relative error to coefficients, eij, has a maximum value of 1.2 x lo- 6per cent. The eigenvalue analysis illustrates that inherited errors change the lowest eigenmodes from energy-free to energy absorbing (or supplying) states. This effect is eliminated if macroscopic equilibrium requirements are met. Even with small inherited errors, this misrepresentation can be expected to result in significant error when elements would have !arge rigid body motion compared with their elastic deformation.
Inherited Error in Finite Element Analysesof Structures 4. MINIMIZING
INHERITED
1215
ERROR EFFECTS
The inherited error includes errors due to the intrinsic inaccuracy of input and errors induced by generation calculations. The latter can be reduced by logic which exercises care in forming coefficients. It has been shown that the effects of the remaining generation errors can be largely eliminated by eliminating equilibrium errors in element stiffness matrices. This section describes the basis for computer logic to accomplish this. Reducing generation errors
The inherited errors due to generation calculations can sometimes be prevented by rejecting self-contradictory input. A common example is the specification of direction cosines to define the orientation of a beam web. Since one parameter is sufficient to define the orientation, use of three direction cosines permits the user to overspecify the problem. Using these data can lead to a nonunitary coordinate transformation matrix. If the user specifies direction cosines to only three digit accuracy the associated inherited error may intrude into the third digit of the matrix coefficients. Another example involves use of logic for describing an axially symmetric structure. Errors implied in the geometric location of a joint on the closure line, because it is located for the first and last elements, may imply errors in equilibrium which are intolerable. A little thought is sufficient to define logic which will preclude errors of this type. If input is acceptable, the strategy for controlling errors should depend on the analysis objective. If maximum accuracy is desired, one strategy should be used; if maximum cost effectiveness, another. Maximum accuracy is usually required in error studies, or when many equations are to be handled compared with the computer precision (e.g. more than 2000 typical loaddeflection equations in single-precision 1108 calculations [9]. Then, assuming exact input, double-precision should be used in all generation calculations before truncating to singleprecision coefficients. This means reading input and converting into double-precision form, performing all generation calculations in double-precision, and storing intermediate results in double precision form. This will require a modest increase in storage requirements during generation and add up to about 15 per cent to calculation time compared with singleprecision operation. Based on worst case errors, [I] this strategy is guaranteed to reduce manipulation errors of generation to a maximum of one part in the last binary place for computers with 12 or more binary places in the mantissa. Only single precision arithmetic can be justified from a cost effectiveness point of view. The analyst rarely knows the structural material elastic constants with less than 2-8 per cent error. Even for simple trusses, manufacturing inaccuracy will induce f24 per cent error in cross sectional areas. Joint positions cannot be expected to be defined with less than H in. error. Double precision is not justified to improve the last 12 bits of the mantissa when the analyst introduces these data which are not accurate beyond the seventh bit. Eliminating equilibrium errors
Minimizing most generation error effects is achieved by enforcing compliance with equilibrium requirements. Even if the double precision is used, it will be advantageous to force this compliance. One way of forcing compliance is to initially represent the stiffness matrix in the computer in the form of its elastic kernel. Then the remaining coefficients of the complete element stiffness matrix can be constructed under the constraint that equilibrium requirements must be satisfied.
1216
ROBERTJ. Mmo.w
The idea and procedure for developing the elastic kernel of the stiffness matrix is not new. (It was used, for example, by Turner et al. [lo] in deriving the stiffness matrix for the triangular membrane element). Not only does it offer a way to enforce equilibrium satisfaction, but it requires fewer calculations to develop the stiffness matrix than other methods. The elastic kernel stiffness matrix (sometimes called the natural stiffness matrix) is expanded to a complete element stiffness matrix using the equilibrium requirements directly. Since as many coefficients are added to each column of the stiffness matrix as there are equilibrium conditions, it is always possible to change the last bit of the new coefficients to enforce equilibrium satisfaction to the precision level used for the elastic part. Note that the kernel matrix should be augmented from the form it assumes just before being combined with other element matrices. Thus all transformations and scaling will have been completed. (Errors in the representation of the elasticity need not be corrected since they will be smaller than the relative errors of input.) Because the augmentation is made to the final form of the stiffness, transformations are imposed on smaller matrices than usual. Thus this approach to generating element stiffnesses avoids equilibrium errors and reduces calculations. Thus gross errors in equilibrium can be eliminated by rejecting contradictory input. By developing the element stiffness matrix by augmenting the elastic kernel, the coefficients can be selected to insure satisfaction of equilibrium. 5. CONCLUSIONS This study of inherited error effects furnishes the following conclusions : 1. Generation errors can induce much large errors in deflection predictions (and hence in stress estimates) than manipulation errors incurred in solving the system loaddeflection equations. (This conclusion implies that if conclusions about manipulation errors are to be drawn from experiments, care must be exercised to avoid or account for inherited errors.) 2. Examination of the difference equations for prismatic beams shows that the largest solution errors due to scaling correlate with the imbalance in moment equilibrium. Depending on the sign of the imbalance, errors may be bounded or unbounded as the number of equations (or joints or elements) increases for a given relative error. When the magnitude of the equilibrium error times the number of equations is small compared with the computer number representation, relative solution errors are directly proportional to equilibrium errors and vary as the square of the number of equations. 3. An eigenvalue/vector analysis of the element stiffness matrix shows that relatively small inherited errors have little effect on changing the eigenvalues or eigenvectors of elastic modes. They significantly change the eigenvalues associated with the ‘rigid’ modes. As a consequence of these changes the bigger the rigid body motion compared with the elastic, and the smaller the element, the greater the solution error induced by the generation error. 4. Most of the solution errors induced by generation errors can easily be eliminated by eliminating the inconsistency implied by the inherited errors. This can be achieved by constructing the equilibrium satisfying element stiffness matrix from the elastic kernel.
Inherited Error in Finite Element Analyses of Structures
1217
Thus the examination shows that relatively large errors are induced by inherited errors that imply an inconsistency in the formulation (violation of equilibrium). Computer logic which eliminates the inconsistency can reduce solution errors to those commensurate with the intrinsic accuracy of input data. REFERENCES [l] R. J. Merow and E. L. PALACOL, Manipulation errors in computer solution of structural equations. NASA Contractors Report, November (1969). [2] R. A. ROSANOFF, J. F. GLOUDEMAN and S. LEVY, Numerical conditioning of stiffness matrix formulations for frame structures. Proc. 2nd Co& on Matrix Metho& in Structural Mechanics, Air Force Flight Dynamic’s Laboratory, AFFDL-TR-68-150, pp. 1029-1060, February (1970). [3] R. J. MELOSH,Characteristics of manipulation errors in solving load-deflection equations. On General Pvpose Finite Element Computer Programs. ASME, New York, pp. 123-142 (1970). [4] A. M. TURING, Rounding-off errors in matrix processes. Q. J. Mech. Phys. 1, 287-308, September (1948). [5] R. J. MELCI~H and E. L. PALACOL, Manipulation errors in 6nite element analysis of structures. National Aeronautics and Space Administration, NASA CR-1385, Washington, D.C., August (1969). [6] R. H. MACNEAL,77reNASTRAN Theoretical Manual. National Aeronautics and Space Administration, NASA SP-221, Washington, D.C., October (1969). [7] W. E. CARROLL and R. M.BARKER,A theorem for optimum tinite element idealizations. ht. J. Solids Struct 9. 883-895. [8] G. M. MCNEICEand P. V. MARCAL,Optimization of tinite element grids based on minimum potential
energy. Brown University Report 7, Div. of Engineering, June (1971). [9] R. J. Me~osrr, Manipulation errors in finite element analysis. In Recent Advances in Matrix Analysis of Structures, pp. 857-877 University of Alabama Press, Alabama, (1971). [lo] M. J. TURNER,R. W. CLOUGH,H. C. MARTINand L. J. TOPP, Stiffness and deflection analysis of complex structures. J. Aero. Sci. 23 (9), 805423 September (1956). (Received 24 February 1972)
1973. Printed in Great Britain
INHERITED ERROR IN FINITE ELEMENT ANALYSES OF STRUCTURES-f ROBERT J. MELOSH 3939 Fabian Way, Palo Alto, California, 94303, U.S.A. A~~~-Mathematicians are quick to point out that round-off and truncation errors induced by the digital computer are only half of the manipulation errors in numerical analyses. The other half am the errors in quantizing the mathematical problem for computer solution: errors inherited for the equation solving process. This paper examines the relevance of inherited errors in structural analyses using the finite element concept and the digital computer. It illustrates error magnitudes using numerical experiments of simple structures. It constructs a theory explaining the errors for these systems. Having identified the most significant controllable computer error source, it describes a process for minimizing its contribution to the inherited error. The paper concludes that orders of magnitude between errors reported by various investigators can be explained by differences in inherited error. The most significant effect of these errors can be identified with inconsistency in problem formuiation. These inconsistencies can be eliminated by exploiting the existence of rigid body states in the finite eIement models. Thereby, solution errors introduced by inherited errors can be reduced to intrinsic errors in parameters defining the geometry, material characteristics, and boundary conditions.
1. INTRODUCTION How CAN a digital computer develop inaccurate data when relatively few calculations are involved? This can occur in the well known GIG0 (garbage in, garbage out) and infamous PICA (precision-in, critical arithmetic) modes. In structural analysis, it is convenient to divide the analysis at the point where the coefficients of the equations of motion are available. The input errors in the GIG0 mode are ‘inherited’ errors to the equation solving calculations. PICA errors, though they may arise in coefficient generation, more often occur in the equation solving calculations. In critical arithmetic, large relative errors are induced in few calculations because the information of importance is embedded in the lower bits of the floating point mantissa. In structural analyses, this can occur in solving the equations of motion when the equation sequence is bad, [1] or, in displacement method analysis, when one or more structural elements is nearly rigid. Though these conditions can induce large errors, these sources of error are now known. Modern structural analysis codes provide safeguards to sense these conditions and take appropriate action. The importance of inherited errors has long been recognized. Rosanoff et af. [2] in 1968 posed a simple problem, involving but two simultaneous equations, for which a 27 bit mantissa is insufficient to obtain two digit solution accuracy. t Presented at the National Symposium on Computerized Structural Analysis and Design at the School of Engineering and Applied Science, George Washington University, Washington, D.C., 27-29 March 1972. 1205
ROBERTJ. MEJBSH
1206
Figure 1 shows the magnitude of inherited errors in a set of structural problems. It illustrates inherited and manipulation errors in computer analysis of a cantilevered beam. Inherited error estimates are predicated on change in the last binary place of the mantissa [3]. These curves show that even though intrusion occurs only into the last binary place, solution error is comparable to the error induced by all the subsequent arithmetic truncation errors in the solving calculations.
I
0001
0
100
Numbers
FIG. 1. Inherited
I
1
200
300
400
of equations,
and manipulation
500
600
N=2J
error growth for beams.
This paper uses the case method to study inherited error rather than a bounding or probabilistic approach. The bounding approach leads to formulas defining the absolute worst case error. These formulas are often costly to evaluate and very conservative compared with expected errors [4]. The probabilistic approach, on the other hand, is not satisfactory because it only indirectly reflects how bad errors can actually become. The case method uses a worst case, but realistic problem, as the basis for evaluating error characteristics and magnitudes. The problem found worst for PICA mode errors is the cantilever beam [5, 61. Because this structure permits systematic accumulation of error in successive steps of the solution process, it will be assumed to be a worst case system for this study. Choice of the cantilevered beam problem as the case also has secondary advantages. It facilitates an exact analysis of error. It permits comparison of the effects of inherited errors with other effects. The problem of predicting response by the displacement method under a static loading will be considered. In this case, equation solving involves developing the solution of a set of linear, simultaneous equations. Since this is the central mathematical problem in structural analysis; be it linear, nonlinear, time dependent or not; the conclusions should have broad application.
Inherited Error in Finite Element Analyses of Structures
1207
For structural analysis, this paper reviews GIG0 mode errors and their effect. Of primary concern are errors that can be automatically reduced or compensated. The next section describes the relative importance of various sources of inherited error. The third section develops simple error estimation formulas. The fourth discusses ~ni~zing the effects of one source of inherited error. The last section summarizes report conclusions. 2. INHERITED
ERRORS
Of the sources of equation solving inherited error, discretization has justifiably received the most attention. Errors in generating equation coefficients, however, also have a significant effect on the accuracy of response predictions, and should be easier to minimize. This section leads to these conclusions by examining the sources and relative ma~itudes of inherited errors. Error sources
Idealization, numerical modeling, and coefficient generation are the principle sources of errors. Idealization errors are incurred in selecting the mathematical model for the system. Numerical modeling errors are implied in particularizing parameters of the mathematical model so it represents the system of interest. Generation errors arise in calculating the c~~cients in the load deflection equations from numerical input data. In finite element analysis of structures, idealization errors are evoked by selecting the distribution of ‘joints’ and ‘cuts’, electing finite element models to be used, choosing material stress-strain relations and failure criteria and selecting joint force and displacement boundary condition equations. Since implementation of the analysis is almost always directed by a computer code, many of these decisions are implicit in the selection of the computer code itself. Idealization errors are attributable to insufficient detail in the equations used to represent the structure. Lack of detail may be due to insu&ient discretization of the structure or omission of terms in the equations. Insufficient discretization occurs, for example, when the uniform load on a non-prismatic beam represented by a cubic displacement function is portrayed by energy equivalent joint forces and moments. An illustration of errordue to omission of terms is the error due to neglecting shear deformations in analyzing a beam. In a broad sense, idealization errors are errors involved in defining the abstract mathematical equations into which the physical system will be mapped. In this sense then, numerical modelIing errors are the errors in defining those parametersin the equations which are to be the knowns. Thus, in finite element analysis these errors are associated with the definition of geometric parameters (joint locations, element thicknesses, joint eccentricities), material parameters (elastic constants, coefficients of thermal expansion, density, strength limits) and boundary condition parameters (magnitudes of prescribed generalized displacements, load distribution coefficients, joint partial restraint values). Numerical modeling errors are introduced by humans. They may involve errors in judgment on the part of the analyst (assignment errors) or errors in communicating the numbers to the computer (tmnscription or keypunching errors). Generation errors may be caused by manipulation or programming errors. Manipulation errors include round-off and truncation errors evoked by limited precision arithmetic. Programming errors include errors in program logic and formula evaluation. In computer impIemented analysis of structures, it is advantageous to define generation errors as all errors induced from the operation of reading input data until the coefficients inherited
RCXERT J.
I 208
MELOSH
of the load-deflection equation are expressed in a floating point representation. Thus in the displacement method these errors include those induced in card reading, input conversion, scaling, nondimensionalizing, coordinate transformations, and evaluation of formulas for stiffness, loading, and stress coefficients. Of the errors discussed, only idealization errors of the discretization type and generation errors of the manipulative type are candidates for completely automatic error reduction or compensation. Though some automatic correction is possible for the other types of errors, those types will not be considered further here. Error magnitudes
The relative importance of discretization and generation manipulation errors is illustrated by the curves of Fig. 2. This figure shows discretization is the more important source of error for coarse meshes. Inherited error is potentially more important with finer meshes. Whereas discretization errors are bounded and will decrease with increasing mesh refinement, generation manipulation errors can induce unbounded errors in the solution and can monotonically increase with mesh refinement.
;, L z 0) z z G n
005 -
/
(diagonals
00055
600
200
Number
FIG. 2.
only)
of equations,
800
2J
Inherited and discretization error growth for beams.
The discretization error leading to the deflection errors of Fig. 2 are due to approximating the loading with a lumped parameter rather than an energy equivalent loading. The inherited error includes only the truncation errors induced in summing element stiffness matrices. A sound basis for automatic control of discretization error has been laid down by Carrel and Barker [7]. Their important work provides the keystone for increasing the understanding of how to distribute mesh points and what refinement to use. Carroll, Barker. McNeice and Marcal [8] have already illustrated some applications in structural analyses. Thus this paper can concentrate on the generation errors and their effects.
Inherited Error in Finite Element Analyses of Structures
1209
The solution errors of the lower curve of Fig. 2 are by no means the maximum generation induced errors that can occur in cantilevered beam analysis. They reflect only the effects of errors on the diagonal of the stiffness matrix. Figure 3 illustrates the error in deflection for a more representative generation error. This figure shows the error in deflection predictions for the beam as a function of the scale factor used for the analysis for a beam represented by 500 finite elements. The element stiffness matrices were divided by the abscissa values to simulate differences induced by using different dimensional units in defining the problem. The smallest errors provide a measure of the manipulation errors in the solving calculations. The largest errors, which are more than 103 x the smallest, reflect the effect of this inherited error on solution error. The dashed curve of Fig. 2 shows how the error for the factor of 11 varies with the number of equations representing the beam. The form of the curves for other scale factors is the same. Comparing this curve with the lower curve of Fig. 2 leads to the conclusion that this generation error source involves solution errors that are about 102 greater than those due to summation errors. Comparing the errors of the dashed curve of Fig. 2 with the manipulation errors in Fig. 1 shows that these inherited errors induce solution errors that are three times greater than those of manipulation in solving the equations. These data expose the importance of exercising care in performing numerical experiments. The variance in error may provide a basis for explaining the divergent conclusions published about error magnitudes. They illustrate that inherited errors may induce solution errors that are much greater than those due to manipulation errors introduced in solving the equations. In seeking an explanation for the variation of these solution errors with scale factor for beams, a basis for minimizing this type of error for more general finite elements will become evident.
2J = 1000 p = 27 bits R
0
.o + 0 Q
.N ;j E
b E
Element stiffness matrix divisor FIG. 3. Error variation with scale factor.
1210
ROBERT J. MELOSH
3. ERROR ESTIMATION To explain the results in Fig. 3, this section examines the difference equations for the beam and energy relations for the finite element. With this basis, the variation in solution errors is attributed to the inaccuracy in representing rigid body modes introduced by inherited errors. Dlyerence equation solution
The load-deflection relation for the beam finite element neglecting shear deformation is expressed by, 2aEf 3-a
Pl
3aG -6C 3aG
where
2a2F
_
-3aG a2F
ml
66 -3aG
.
(1)
p2
2a21:
m2
is the scalar multiplier is Young’s Modulus, is the bending moment of inertia, r is the element axial length, a wt,Pr are the displacement and force acting normal to the neutral axis at the joint at end 1 and @2sm2 are the angular rotation and moment associated with bending of the joint at end 2 of the element.
;
and
F=l+e,. G=l+e,
(2)
with e, and eG the relative errors introduced by scaling. There are only two distinct relative errors because factors of two among coefficients introduce no change in relative error using a binary floating point number representation. Forming the load-deflection relations for the complete prismatic and uniformly loaded system results in equilibrium equations for the beam, --6w,_ rG-3@,_ raG+12w,G--6w,,
tG+3&
laG=P
(3)
3~,_~G+30,_raF+4B,aF--3w,+tG+t3,,taF=M where P is the generalized normal force at each joint, P=pa3/aBl, and M is the corresponding generalized moment (M=ma2/aEI). With the equations in this form it is seen that the effect of a is merely to scale the loading coefficients. Hence errors in this solution, due to a, are inversely proposal to CL For the problems of interest, the boundary conditions are: w(r=O)=O,
e(r=o)=o,
(4)
1211
Inherited Error in Finite Element Analysesof Structures 3w,_ 1G-3wNG+8N_ -6w,_
1G+6~NG-38N_
laF+20+F=A4 laG-3BNaG=P.
If F= G= a = l., the difference equations are exactly represented. equations (3 and 4) for tip deflection is given by
w,,,=--&3N4+12N3-1N2-2N)
The sol$ion
of
(5)
when m =p = a = 1.0 and N is the number of finite elements. Suppose the solution of equations (3 and 4) are expanded in a power series in /3 where /l=6(F-G)/F
(6)
Then, subtracting equation (5) from the fist terms of this expansion and dividing by equation (5) yields a measure of the relative error of the solution, i.e. (7) where the convergence of the expansion is guaranteed when N2/?< 1. Equation (7) suggests that if NzB>O, the error is bounded. If NzpcO, the error IS unbounded. (The complete expansion will confirm this since eT1p can then be expressed, for a given N, as a constant and the cos NJ/? for 8~0 and a constant and cash N,/p for B>O.) For /INz% 1, and /? ~0; equation (7) shows that the tip deflection error increases approximately as the square of the number of finite elements. For /?N2< 1, and b
The solution of the difference equations furnishes the basis for explaining the erratic results depicted in Fig. 3. For this series of problems, the changing of scale factors results in significant changes in the relative importance of manipulation and inherited error effects. The cases considered can be divided into two classes: those which involve negligible equilibrium imbalance as measured by /3, and those where this imbalance is significant. Each case can be classified by calculating a-l and /? for the relevant divisor. These calculations yield the data given in Table 1. For this set of divisors, the equilibrium imbalance is zero. These data were developed on the same computer as used to generate the data of Fig. 3. Table 1 also lists the solution error for these cases.
1212
ROE~T 3. MEL~SH TABLE1. Errorst due to misrepresenting elasticity: Divisor l/3 115 l/7 1,2,4,8,16 7,14 15
a Error
Solution error:
0 0 0 0 7.5-7 3.7-7
622 -6 629 -6 6.27 -6 &Sg-6 6.17 -6 6.28 -6
t Relative errors from experiments with a 27 bit mantissa and truncation arithmetic. $ All errors involve numbers which are less than the exact values.
These data show that for these cases, solution errors are the result of only manipulation errors in the solution process. Solution errors are essentially unaffected by changing c1 errors from zero to nonzero. Since the relative error in the solution is three orders of magnitude greater than the CLerror, the error in misrepresenting stiffness for these cases is due to truncation errors in decomposition, forward, and back substitution. Thus the solution errors for all the small error cases, the circled points on Fig. 3, are due to manipulation errors. Data for the cases in which the equilibrium imbalance is not zero are listed in Table 2. Because the CIerrors for these cases make a negli~ble contribution these c1errors have been omitted. For convenience, the divisors incurring negative and positive /I errors have been grouped separately. Figure 4 provides a log-log plot of the solution errors as a function of the number of equations. Curves are included for four of the sets of divisors listed in Table 2. The
Number
of
equations,
N
FIQ. 4. Solution error due to equilibrium imbalance.
Inherited Error in Finite Element Analyses of Structures
1213
divisors with negative errors are shown by continuous curves; those for 3, 6, and 12 by a dashed curve. The relative manipulation error has been subtracted from the solution error.
Thus results in Fig. 4 pertain to only equilibrium errors. The continuous curves demonstrate that the negative #I cases have solution errors that increase as N2. The errors increase monotonically with 8, as prescribed by the analytical results. The positive /I case shows that as N increases, the error increases monotonically. However, the errors of sign opposite to the manipulation do not necessarily compensate for truncation, probably because they arise in a different way. This is illustrated by examining the data in Table 2. TABLE2. Errorsg due to equilibrium imbalance Divisor
/I Error
17 9 5.10 11 3,6.12 13
-158-S -1.67-8 -3.71-8 -6.14-s 2.23 -7 4.84 -8
Solution Error 1.48 -5 152-s 7.21-3 1.19 -2 150-s 6.18-6
5 Relative errors from experiments with a 27 bit mantissa and truncation arithmetic. 11All errors involve numbers which are less than the exact values.
Thus the large errors in Fig. 3 are caused by ‘negative’ equilibrium imbalances. The small ones involve manipulation errors or ‘negative’ interaction of manipulation and positive equilibrium imbalance. Using the analytical results, Table 3 provides the data for comparing the effects on beam solution errors of rounding, rather than truncating, in arithmetic. This table cites the relative error in satisfying equilibrium 8, for the cases of interest. Truncation errors are repeated for convenience. As would be expected, these data support the conclusion that maximum errors will be less with rounding than with truncation. Because the equilibrium errors often vary significantly with the arithmetic mode used, the arithmetic mode should not be overlooked. This emphasizes again that casual selection of computer tests can add to the confusion rather than increase understanding. TABLE3. Truncation and rounding equilibrium errorst Divisor
Truncation B error
1,2,4,8,16 3,6,12 5,lO 7,14 9 11 13 15 17 t Relative errors with a 27 bit mantissa.
Rounding jI error
0
0
2.23 -8 -3.71-s 0 -1.67-s -6.14-s 4.84 -8 0 -1.58-s
2.23 -8 3.71-s -2.61-s -1.67-s 1.96-s -2.42-s 0 6.18-11
1214
ROBERTJ. MELQSH
Eigenvalue analysis
The implications of inherited errors may be better understood by considering their effect on the eigenvalues of the stiffness matrix. Since a is a matrix scalar multiplier, it scales all eigenvalues without changing eigenmodes. Independent errors to coefficients of the ma’trix change both eigenvalues and modes. Assuming truncation errors, small errors to coefficients have little relative effect on zero terms. On a p bit computer with inherited error in the pth bit the elastic modes are effected in the pth bit. Similarly, only the pth bit of each elastic mode eigenvalue is reduced. The zero eigenvalues associated with rigid body modes become non-zero. They may be either small positive or negative value depending on the magnitudes of relative error to each coefficient. Thus, instead of being positive semi-definite, the element stiffness matrix can become indefinite. It is noted that if the forces and moments of each column of the stiffness matrix satisfy macroscopic equilibrium requirements for the element, this is sufficient to insure zero frequency rigid modes. This is not sufficient to avoid errors in the elastic eigenvalues and modes. The energy error due to the inherited error in the element matrix is proportional to the eigenvalues. Thus, expressing the ratio of translational and rotational eigenvalues to that of the lowest elastic mode provides an error indicator. Suppose that each of the coefficients may have a different inherited error with a relative value of eU where the subscripts relate to the row and column numbers of the matrix of equation (1). Then the indicators are:
ell--6e12+6e13-
6el,+4e22-6e2,+4e,,+3e,J
- 6e,, + 4e,,jqL qd
where IT is a relative measure of rigid translation error, ZR is a relative measure of rigid rotation error, qT and q, are rigid mode partrcipation factors, and qxr is the participation factor for the more flexible elastic mode. Note that Zr=ZR =0 if matrix coefficients satisfy equilibrium. These formulas confirm that when the rigid motion of the element is large compared with its elastic deformation, small errors in equilibrium can imply significant energy in ‘rigid’ modes. The indicator for translation projects large errors when very short elements are involved. This situation can be illustrated by the cantilevered beam solution of the test problem with 500 finite elements. For the scaling error rigid translation is exactly represented. Rigid rotation errors, however, would be expected to be particularly large for the tip section. An energy analysis shows the rigid mode energy in the tip element is 0.62 x 10-Z per cent of the total external work even though the relative error to coefficients, eij, has a maximum value of 1.2 x lo- 6per cent. The eigenvalue analysis illustrates that inherited errors change the lowest eigenmodes from energy-free to energy absorbing (or supplying) states. This effect is eliminated if macroscopic equilibrium requirements are met. Even with small inherited errors, this misrepresentation can be expected to result in significant error when elements would have !arge rigid body motion compared with their elastic deformation.
Inherited Error in Finite Element Analysesof Structures 4. MINIMIZING
INHERITED
1215
ERROR EFFECTS
The inherited error includes errors due to the intrinsic inaccuracy of input and errors induced by generation calculations. The latter can be reduced by logic which exercises care in forming coefficients. It has been shown that the effects of the remaining generation errors can be largely eliminated by eliminating equilibrium errors in element stiffness matrices. This section describes the basis for computer logic to accomplish this. Reducing generation errors
The inherited errors due to generation calculations can sometimes be prevented by rejecting self-contradictory input. A common example is the specification of direction cosines to define the orientation of a beam web. Since one parameter is sufficient to define the orientation, use of three direction cosines permits the user to overspecify the problem. Using these data can lead to a nonunitary coordinate transformation matrix. If the user specifies direction cosines to only three digit accuracy the associated inherited error may intrude into the third digit of the matrix coefficients. Another example involves use of logic for describing an axially symmetric structure. Errors implied in the geometric location of a joint on the closure line, because it is located for the first and last elements, may imply errors in equilibrium which are intolerable. A little thought is sufficient to define logic which will preclude errors of this type. If input is acceptable, the strategy for controlling errors should depend on the analysis objective. If maximum accuracy is desired, one strategy should be used; if maximum cost effectiveness, another. Maximum accuracy is usually required in error studies, or when many equations are to be handled compared with the computer precision (e.g. more than 2000 typical loaddeflection equations in single-precision 1108 calculations [9]. Then, assuming exact input, double-precision should be used in all generation calculations before truncating to singleprecision coefficients. This means reading input and converting into double-precision form, performing all generation calculations in double-precision, and storing intermediate results in double precision form. This will require a modest increase in storage requirements during generation and add up to about 15 per cent to calculation time compared with singleprecision operation. Based on worst case errors, [I] this strategy is guaranteed to reduce manipulation errors of generation to a maximum of one part in the last binary place for computers with 12 or more binary places in the mantissa. Only single precision arithmetic can be justified from a cost effectiveness point of view. The analyst rarely knows the structural material elastic constants with less than 2-8 per cent error. Even for simple trusses, manufacturing inaccuracy will induce f24 per cent error in cross sectional areas. Joint positions cannot be expected to be defined with less than H in. error. Double precision is not justified to improve the last 12 bits of the mantissa when the analyst introduces these data which are not accurate beyond the seventh bit. Eliminating equilibrium errors
Minimizing most generation error effects is achieved by enforcing compliance with equilibrium requirements. Even if the double precision is used, it will be advantageous to force this compliance. One way of forcing compliance is to initially represent the stiffness matrix in the computer in the form of its elastic kernel. Then the remaining coefficients of the complete element stiffness matrix can be constructed under the constraint that equilibrium requirements must be satisfied.
1216
ROBERTJ. Mmo.w
The idea and procedure for developing the elastic kernel of the stiffness matrix is not new. (It was used, for example, by Turner et al. [lo] in deriving the stiffness matrix for the triangular membrane element). Not only does it offer a way to enforce equilibrium satisfaction, but it requires fewer calculations to develop the stiffness matrix than other methods. The elastic kernel stiffness matrix (sometimes called the natural stiffness matrix) is expanded to a complete element stiffness matrix using the equilibrium requirements directly. Since as many coefficients are added to each column of the stiffness matrix as there are equilibrium conditions, it is always possible to change the last bit of the new coefficients to enforce equilibrium satisfaction to the precision level used for the elastic part. Note that the kernel matrix should be augmented from the form it assumes just before being combined with other element matrices. Thus all transformations and scaling will have been completed. (Errors in the representation of the elasticity need not be corrected since they will be smaller than the relative errors of input.) Because the augmentation is made to the final form of the stiffness, transformations are imposed on smaller matrices than usual. Thus this approach to generating element stiffnesses avoids equilibrium errors and reduces calculations. Thus gross errors in equilibrium can be eliminated by rejecting contradictory input. By developing the element stiffness matrix by augmenting the elastic kernel, the coefficients can be selected to insure satisfaction of equilibrium. 5. CONCLUSIONS This study of inherited error effects furnishes the following conclusions : 1. Generation errors can induce much large errors in deflection predictions (and hence in stress estimates) than manipulation errors incurred in solving the system loaddeflection equations. (This conclusion implies that if conclusions about manipulation errors are to be drawn from experiments, care must be exercised to avoid or account for inherited errors.) 2. Examination of the difference equations for prismatic beams shows that the largest solution errors due to scaling correlate with the imbalance in moment equilibrium. Depending on the sign of the imbalance, errors may be bounded or unbounded as the number of equations (or joints or elements) increases for a given relative error. When the magnitude of the equilibrium error times the number of equations is small compared with the computer number representation, relative solution errors are directly proportional to equilibrium errors and vary as the square of the number of equations. 3. An eigenvalue/vector analysis of the element stiffness matrix shows that relatively small inherited errors have little effect on changing the eigenvalues or eigenvectors of elastic modes. They significantly change the eigenvalues associated with the ‘rigid’ modes. As a consequence of these changes the bigger the rigid body motion compared with the elastic, and the smaller the element, the greater the solution error induced by the generation error. 4. Most of the solution errors induced by generation errors can easily be eliminated by eliminating the inconsistency implied by the inherited errors. This can be achieved by constructing the equilibrium satisfying element stiffness matrix from the elastic kernel.
Inherited Error in Finite Element Analyses of Structures
1217
Thus the examination shows that relatively large errors are induced by inherited errors that imply an inconsistency in the formulation (violation of equilibrium). Computer logic which eliminates the inconsistency can reduce solution errors to those commensurate with the intrinsic accuracy of input data. REFERENCES [l] R. J. Merow and E. L. PALACOL, Manipulation errors in computer solution of structural equations. NASA Contractors Report, November (1969). [2] R. A. ROSANOFF, J. F. GLOUDEMAN and S. LEVY, Numerical conditioning of stiffness matrix formulations for frame structures. Proc. 2nd Co& on Matrix Metho& in Structural Mechanics, Air Force Flight Dynamic’s Laboratory, AFFDL-TR-68-150, pp. 1029-1060, February (1970). [3] R. J. MELOSH,Characteristics of manipulation errors in solving load-deflection equations. On General Pvpose Finite Element Computer Programs. ASME, New York, pp. 123-142 (1970). [4] A. M. TURING, Rounding-off errors in matrix processes. Q. J. Mech. Phys. 1, 287-308, September (1948). [5] R. J. MELCI~H and E. L. PALACOL, Manipulation errors in 6nite element analysis of structures. National Aeronautics and Space Administration, NASA CR-1385, Washington, D.C., August (1969). [6] R. H. MACNEAL,77reNASTRAN Theoretical Manual. National Aeronautics and Space Administration, NASA SP-221, Washington, D.C., October (1969). [7] W. E. CARROLL and R. M.BARKER,A theorem for optimum tinite element idealizations. ht. J. Solids Struct 9. 883-895. [8] G. M. MCNEICEand P. V. MARCAL,Optimization of tinite element grids based on minimum potential
energy. Brown University Report 7, Div. of Engineering, June (1971). [9] R. J. Me~osrr, Manipulation errors in finite element analysis. In Recent Advances in Matrix Analysis of Structures, pp. 857-877 University of Alabama Press, Alabama, (1971). [lo] M. J. TURNER,R. W. CLOUGH,H. C. MARTINand L. J. TOPP, Stiffness and deflection analysis of complex structures. J. Aero. Sci. 23 (9), 805423 September (1956). (Received 24 February 1972)