On reliability in finite element computations

Compurers & S~ucrures Vol. 39, No. 6, pp. 729-134, Printed in Great Britain.

ON RELIABILITY

004s7949/91 $3.00 + 0.00 Pergamon Press plc

1991

IN FINITE

ELEMENT

B. A. Center for Computational

COMPUTATIONS

SZAB6

Mechanics, Washington University, St Louis, MO 63130, U.S.A. (Received 21 May 1990)

Abstract-Direct computation of engineering data from finite element solutions with reasonable guarantee of reliability is discussed with reference to practical engineering problem solving. The discussion focuses on a model problem for which experimental and numerical data are available.

INTRODUCTION

In structural and mechanical engineering practice problems solved by the finite element method are typically posed as follows: Given some topological data; material properties; constraints, and loading, compute a set of data, such as deflections; natural frequencies; stress maxima; stress intensity factors, etc. within a specific error tolerance. Such assignments are typically given by a senior engineer, responsible for decisions relating to design or certification, to a junior engineer who is to perform the analysis. The question is then, given the available numerical or experimental methods, which method should the junior engineer choose so that the required data are computed with reasonable certainty that the error tolerance has not been exceeded and at the smallest overall cost. Most problems in structural and mechanical engineering practice are such that the exact solution is analytic? everywhere within the solution domain and on the boundaries of the solution domain, with the exception of a finite number of points. The points where the solution is not analytic are coincident with the vertices of the finite element mesh. Problems which have this property have been classified in [l, 21 as problems belonging to Class B. On the basis of theoretical and practical considerations it is argued herein that the p-version of the finite element method is the most efficient general numerical procedure for solving problems in this class. The advantages of the p-version are particularly obvious when it has to be shown that the computed data are close to their exact values. DESCRIFTION

OF THE MODEL

PROBLEM

There are very few model problems for which experimental data as well as analytical solutions by various methods are available. Among the few is the problem posed by Floyd [3]. This problem has been studied experimentally [3,4] and solved by varit A function is analytic in a point if it can be expanded into a Taylor series about that point.

ous numerical methods, including the boundary integral method [336]; the h-version of the finite element method [3-71, and the p-version of the finite element method [S, 91. This model problem is concerned with finding stresses in a thick-walled axisymmetric pressure vessel caused by an internal pressure of 2.61 psi (18.0 kPa). The dimensions of the pressure vessel are shown in inches (1 inch = 2.540 cm) in Fig. 1. The z-axis is the axis of symmetry. The plane z = 3.15 is a plane of symmetry. For this reason only half of the generating section is shown. The material is Araldite, the modulus of elasticity (E) is 1435 psi (9.894 MPa) and Poisson’s ratio is 0.49. In the case of this model problem it is known a priori that the exact solution is not analytic in the end points of the circular fillet, and in point F. Since in these points the strains are finite, by the classification given in [ 1,2], this problem is weakly in Category B. The original problem statement in [3] requires computation of the radial, circumferential and axial stress components along line AE. In this paper the problem is solved with a different problem statement formulated to represent more nearly the way similar problems are posed in engineering practice. The problem statement to be discussed in the following is: ‘Compute the location, orientation and magnitude of the largest principal stress and give evidence that the computed values are accurate to within five percent relative error.’ FUNDAMENTALS

AND NOTATION

Corresponding to any properly formulated mathematical model is an exact solution, denoted by uEX, which depends on the choice of formulation and the data which characterize the geometry, material properties loading, and constraints, but is independent of the discretization. Corresponding to a particular choice of discretization is a finite element solution, * wFE,and a number N, called the number of degrees of freedom, which is the number of linear algebraic equations one has to solve in order to obtain i&. One is typically interested in certain functionals @Jo, (i = 1,2, . . .). Since uEXis not known, only @&I,) can be computed. If the problem is well

129

730

B. A.

%ABC!l

element solution is that function from 3 which minimizes the potential energy

-7 E

L-6.0”--’

IT(&)

b ;;:

-6 ,: Alr

Fig. I. The generating section and the boundary conditions.

formulated and the discretizations are properly selected then @&I,) + @i(u,) as N + co. The computations are reliable if it is known with a reasonable degree of certainty that

where ti are small numbers, typically in the range 0.01-0.0.5, i.e. l-5%. In engineering computations accurate estimates of 7, are not essential but it is important to know whether 2, are large or small. If only one finite element solution is available then it is impossible to know in general whether ri are large or small. In order to estimate zi a sequence of finite element solutions, corresponding to a sequence of discretizations, must be obtained. In the following discussion the displacement formulation is considered. The choice of a particularly discretization involves the creation of a set of functions S on the solution domain ft. The set S, called the finite element space, is characte~~d by the mesh A, the mapping functions Q, that is the functions which map standard finite elements onto the elements of the mesh, and the function spaces defined on the standard elements, called standard spaces. Most commonly the standard spaces are polynomial spaces characterized by the polynomial degree p. Hence, formally; S = s@, A, Q, p). By definition, ~h;r~Q?$2’~ . . . IQ~(~))~ P’ {PI.&, . . .PICl(A)) represents the transformation x = @)({, r~), y = @)(c, TV)between a standard element (e.g. triangle or quadrilateral) and the kth element of the mesh; M(A) is the number of elements in the mesh. Typically the standard polynomial space characterized by degree p is either the space Y’Por the space YJ’Qwhich are defined in [2] for various standard elements. For example, on the standard quadrilateral element the space 9~’ spans all monomial terms of degree p plus the term &I for p = 1 and the terms {Pq, &’ for p 2 2. In contrast, the space YpJ spans the sets of monomials { 1, r, t*, . . , rp} (4 d, * . , f), and their products. In [8] the space YP was used whereas in [9] the space Yfip was used. The subset of functions in S which satisfies the kinematic boundary conditions is denoted by 3. The number of degrees of freedom N is the number of linearly independent functions in 5. The finite

= nltl fl(ic)

(2)

If a hierarchic sequence of finite element spaces is constructed, that is S, c S2c . . I , then n(u,) converges to U(u,) monotonically. This monotonic convergence can be exploited to estimate n(u,) and, using the estimated value of n(u,), error estimates in energy norm can be computed. By definition, the energy norm of a displacement function u is the square root of the strain energy of II, which is usually denoted by 11 u //E, The relative error in energy norm is closely related to the relative error in the root mean square measure of stresses [2, lo]. Denoting the estimate of n(u,) by Il’*, an estimate for the error in energy norm can be computed from

This estimate is as~ptoticaIIy correct and has been shown to work well in the preasymptotic range [2, 10, I 1J. Sequences of discretization can be constructed in various ways: by mesh refinement; changing the mapping; changing the standard spaces; increasing the degree of the standard polynomial space, or any combination of these. For reasons of implementation most commonly the mesh is refined, or the degree of the standard polynomial space is increased, or mesh refinement is combined with increase of the polynomial degrees. These approaches are respectively call h-extension, p-extension and up-extension. When aspects of implementation rather than reduction of error are discussed then the word version instead of extension is used. SOLUTION OF THE MODEL PROBLEM

The solution presented in the following was obtained with the p-version finite element computer program MSC/PROBE. This program uses the standard space Yp and p may range from 1 to 8. Mapping is by the blending function method ]2]. The finite etement mesh, consisting of 17 elements, is shown in Fig. 2. In designing this mesh certain CIpriori knowIedge about uEX was taken into consideration: It is

I

Y

Fig. 2. Finite element mesh consisting of 17 elements.

731

Finite element computations known that if the fillet radius were zero then there would be a strong corner singularily in the point r = 5.4, z = 1.5 and a geometric mesh, that is a mesh graded in geometric progression with a common factor of 0.15, would be used in conjunction with p-extension[2, 10, 111. If, on the other hand, the solution were very smooth (for example the fillet radius would be large) then a nearly uniform mesh would be best. The mesh design shown in Fig. 2 is based on the expectation that, because the fillet radius is small, the solution is likely to be closer to the singular case than to the smooth case. The convergence of the strain energy U(ic,) and the estimated relative error in energy norm is shown in Table 1. It is seen that as p is increased the functional U(u,,) converges to a limit value which is U(u,,). In addition to the estimated relative error in energy norm the estimated rate of convergence, which is by definition the slope of the log U(u, - uFE) versus log N curve, is shown in the column 2fi. The fact that the value of 2/l is decreasing with increasing p indicates that the mesh is underrefined for the higher p values. In hp-extensions the mesh is nearly optimal for each p-level and the rate of convergence is exponential. In that case 28 increases with p. In the original problem statement [3,4], the goal is to compute the stress components along line AE. Since the solution corresponding to this problem statement was already presented in [8], only a few remarks concerning aspects of releability are made in the following: In point A the exact value of the minor principal stress is known to be -2.61 psi and the orientation of the major principal stress 8, can be computed from the topological data to be 56.9 counterclockwise from the r-axis. The computed values of r~,, a>, a3 and 0, are shown in Table 2. The relative error in oZ is also shown. The convergence of the major principal stress is evidently stronger, hence it is reasonably certain that the value computed for p = 8, (a, z 92.9 psi) is within 5% (and probably closer to 1%) of its exact value. This demonstrates a simple test of consistency of the computed stress components. Such tests are more stringent than the estimated relative error in energy norm, which is a global measure of error. This is particularly true when Poisson’s ratio is close to 0.5, as in this case. The large consistency error at p = 4 is Table 1. p-Convergence of the strain energy and estimated relative error in energy norm (percent) P

N

1 48 2 130 3 216 4 336 5 490 6 678 7 900 8 1156

(lbf

Est. 28

Est. rel. error

0.60621 1.75370 1.84831 1.88290 1.88808 1.88927 1.88959 1.88979

2.25 2.33 3.95 3.26 2.48 1.48 1.48

82.42 26.87 14.89 6.22 3.36 2.25 1.82 1.51

Wu,) in)

Table 2. Convernence of the stress comnonents in noint A P

N

~1

4

336

100.2

asso

*2

5.31

4

(303.4%j

46.1

57.0

(0.2%)

5 490 95.0

-4.42

‘(69.3%j

40.8

57.0 (0.2%j

6 678 7 900 8 1156

-2.41 -2.78 -2.72

(7.7%) (6.5%) (4.2%)

39.0 38.6 38.7

57.0 (0.2%) 56.9 (0.0%) 56.9 (0.0%)

cc

93.2 92.9 92.9

00

-2.61

56.9

(0.0%)

(0.0%)

a manifestation of locking which occurs at low p-levels when the material is nearly incompressible [12]. The computed values of the location and magnitude of the maximum principal stress are shown in Table 3. Computation of stress maxima is performed by subdividing the standard element into an arbitrary (user-specified) uniform grid, known as the ‘data mesh’. In this case a 12 x 12 data mesh was used. The stress values were computed in each grid point and the maximum value was selected from the resulting sets computed for each element. The stress contours, shown in the neighborhood of the fillet in Fig. 3, were computed similarly. On interelement boundaries the stresses were computed twice; one for each element, hence discontinuities would be visible. No significant stress discontinuity can be observed in Fig. 3. As seen in Table 3, the location, magnitude and orientation of the first principal stress exhibit strong convergence with respect to increasing p. It can be easily checked that the location of the maximum stress (r = 5.387, z = 1.588) is on the boundary of the circular fillet. This implies that the exact value of a2 is -2.61. Knowing that the orientation has to be tangential to the circular boundary, the data can be checked for consistency. In particular, the orientation of tangent to the fillet at r = 5.387, z = 1.588 is 65.7” measured from the r-axis. This is close to the computed orientation of 65.6”. The results of consistency tests, that is tests based on the assumption that r = 5.387, z = 1.588 are accurately known, are shown in the last row of Table 3. The problem statement for this particular model problem requires ‘evidence that the computed values are accurate to within 5% relative error’. The following evidence has been developed: (1) The estimated Table 3. Convergence of the location, orientation, and magnitude of the maximum principal stress (12 x I2 data mesh) P

N

3 216 336 4 5 490 6 678 7 900 8 1156 co

co

r

Z

64

(in)

5.356 5.362 5.392 5.387 5.387 5.387

1.544 1.550 1.603 1.588 1.588 1.588

(5.387)

(1.588)

01 (Psi) 129.4 104.6 97.0 94.1 93.9 93.7

c2

@4 51.69 14.48 2.86 -2.55 -2.40 -2.57 (-2.61)

8,

(deg.) 47.0 49.4 71.8 65.6 65.6 65.6 (65.7)

B. A. &Ad

732

RBcO-

2.0202*00 0,1202+02 0.222002 O.B20&02

EFtnI -

O.woE~2 0.52oE+Q2 0.2OoE*o2 O.'IoOE+QI 0.2202*02

J-

0.000L+O2

I(=

O.lOOE*OS

___---

Fig. 3. Mesh detail and contours of the first principal stress (psi) in the vicinity of the fillet (p = 8).

relative error in energy norm was shown to be less than the specified error tolerance for p = 5,6,7,8. (2) All required data, that is the location, magnitude and orientation of the maximum principal stress, were found to converge strongly to their limit values with respect to increasing p. (3) The orientation of the major principal stress passed a consistency test. (4) There are no significant discontinuities in the first principal stress in interelement boundaries. THE COST OF RELIABILITY

Obviously, obtaining a converging sequence of solutions is more expensive than obtaining one solution only. Nevertheless, it is shown in the following that in the case of p-extensions the additional cost is small.

quently, the operations for which the times are shown in brackets were not actually performed. It is seen that stiffness generation and solution for p=l,2,.., 8 was 203 CPU sec. The overhead, that is all other operations, required 42 CPU sec. Extraction, the computation of the error estimates and the data shown in Table 3, required 85 CPU sec. Of the total time of 330 CPU set only 43 CPU set were required to obtain the sequence of solutions from 1 to 7. This is about 15% of the time required to obtain the solution for p = 8 alone. Thus, the difference in the costs of obtaining a sequence of solutions, and obtaining a single solution for the highest p only, which is the cost of reliability, is not large. By comparison, in adaptive h-extensions, with the goal Table 4. Computer time for various operations (CPU set)

Computer resource requirements

P

A VAX station 3100, model 40, VMS-operating system and a VT 100 terminal were used to solve this problem. The CPU times needed for the various operations are shown in Table 4. The stiffness generation times for p ranging from 1 to 7 are shown in brackets. This is because hierarchic basis functions were used, hence the stiffness matrices for p = 1 to p = 7 were embedded in the stiffness matrices corresponding to p = 8. The element stiffness matrices were computed for p = 8 only, the other stiffness matrices were obtained simply by deleting rows and columns from the stiffness matrix for p = 8. Conse-

1 48 2 130 3 216 4 336 5 490 6 678 7 900 8 1156

N

Stiffness

1 I:,’ I;; (15) ‘,:;; 127

Stiffness + solution: Overhead: Extraction: Total:

Solution

Sum

I

1 2 3 6 11 19 33

: 3 6 11 19 160 203 42 85 330

133

Finite element computations

to reduce the error in energy norm to within a tlxed tolerance, the total cost of computation is typically three times the cost of computation on the final mesh P31. Memory requirements in terms of disk space were as follows: Batch input file size: 11 blocks; database: 590 blocks; temporary space, needed only for the duration of the run; 4281 blocks. Cost estimate

In estimating the cost of solving this model problem, one must consider engineering time; computer time, and memory requirements. The cost estimate presented here is based on the following assumptions: Cost of engineering time: $40.00 per hour. Cost of computer time: $0.005 per CPU sec. This estimate is based on 10% average utilization. Cost of disk space: $0.007 per block per year. This problem was solved by an experienced user of MSC/PROBE with average typing skills. The data entry using a VT 100 terminal, required 35 min. Upon completion of the data entry a check run was performed at p = 1. This required 1 min. Waiting time for the solution to be completed for p ranging from 1 to 8 was 7min, and computation of engineering data from the finite element solution, that is computation of the data in Tables 1 and 3, required 4 min. Thus, the total engineering time was 47 min. In estimating the cost of disk utilization it was assumed that the batch input file remains on the disk for 6 months; the database for 1 month and the temporary space is used for the duration of the run only. The cost estimates are summarized in Table 5. Of course, it could have happened that the desired accuracy could not be demonstrated for the initial mesh within the available range of p-extension (p = 1,2,. . . ,8). In that case the mesh would have had to be refined further and the sequence of solutions computed again, nearly doubling the cost. Most of the additional expense would have been related to the analyst spending time to make the necessary changes in the topological section of the database. (Revisions of the mesh are guided by local indicators of error based on tests of equilibrium and continuity 121.This process can be automated, however automated mesh refinement is not available at present.) This underlines the importance of proper inital mesh design based on a priori information about uEX[ 11.

Total

41 min 330 CPU set 11 blocks 590 blocks 428 1 blocks

REMARKS

In this relatively simple model problem approximately 95% of the cost of analysis was associated with engineering time. This high percentage of cost related to engineering time is typical for finite element analyses in general. In engineering practice the time available for performing finite element analyses is generally limited not only by the cost of analysis but also by scheduling considerations. Often analysts have time to produce input data for one finite element mesh only. Although an increasing number of conventional finite element programs offer some estimate of the error in energy norm, the reliability of these estimates has not been fully established; the indicated error may not be closely related to errors in the engineering data of interest, and the time and cost associated with repeating analysis with revised meshes tend to be high. In p-extensions the finite element mesh is usually much simpler than in h-extensions and the time required for data preparation is substantially less. Therefore, the overall cost of analysis is substantially less. Importantly, p-extension produces a strongly converging sequence of solutions at a small marginal cost and in a small amount of additional time. This makes reliable estimation of the relative error in energy norm possible. The accuracy and reliability of the engineering data of interest as ascertained by observation of the convergence of the quantities of interest to their limit values and by performing simple consistency tests. The importance of proper initial mesh design, guided by a priori information about the exact solution of the problem, has been noted. It is not essential that the initial mesh be close to optimal, however. This is because discretization errors are not very sensitive to mesh design when p-extensions are used. Thus, p-extensions compensate for shortcomings in initial mesh design to a certain extent. In this paper reliability was considered under the assumption that the problem was properly formulated and no errors occurred in the preparation of the input data. It often happens, however, that the problem is improperly formulated or some errors are present in the input data. The presence of errors caused by improper problem formulation and the presence of many kinds of input error can be detected by p-extensions failing to converge. This is an important additional consideration favoring the choice of p-extensions. Acknowledgement-This work was supported in part by the National Science Foundation under grant No. DMC-

Table 5. Estimated cost of analysis Engineering time Computer time Disk space Batch input file Database Temporary

CLOSING

S31.33 Sl.65 $0.04 $0.08 $0.00 $33.10

8606533.

REFERENCES

1. B. A. Szabo, The use of a priori estimates in engineering computations, Comput. Meth. Appl. Mech. Engng, in press (1990).

B. A.

734

2. B. A. Szabo and I. BabuSka, Finite Element Analysis. John Wiley, New York (to be published, 1991). 3. C. G. Floyd, The determination of stresses using a combined theoretical and experimental analysis approach. Proc 2nd Int. Conf on Comp. Meth. and Exp. Meth.

(Edited by C. A. Brebbia). Springer, Berlin (1984). 4. P. W. Dickenson and C. G. Floyd, A comparative study of different techniques in the stress analysis of a nuclear component. Proc. 8th Int. Conf. Structural Mechanics and Reactor Technology, Brussels, paper B8/6 (1985). 5. P. W. Dickenson, Results from a pump bowl analysis. A nossible finite element test model. NAFEMS memorandum (1985). 6. C. A. Brebbia and J. Trevelyan, On the accuracy and convergence of boundary element results for the Floyd pressure vessel problem. Comput. Struct. 24, 513-516 (1986).

7. T. Sussman and K.-J. Bathe, Studies of finite element procedures-on mesh selection. Comput. Struct. 21, 257-264 (1985).

SZAB~

8. J. Schiermeier and B. Taylor, An alternate algorithm for FEA. Comout. Mech. Engng 5. 45-51 (1987).

9. K. Kato and K.-J. Bathe, &‘the use-of higher-order finite elements in structural mechanics. Finite Element Research Group. report 89-4, MIT (1989). 10. B. A. Szabo, Mesh for thep-version of the finite element method. Comput. Meth. Appl. Mech. Engng55, 181-197 (1986).

11. B. Szabo, Estimation and control of error based on P-convergence. In Accuracy Estimates and Adaptive Refinements in Finite Element Computations (Edited by I. BabuSka, J. Gago, E. R. de A. Oliveira and 0. C. Zienkiewicz), pp. 61-78. John Wiley (1986). 12. B. A. Szabo. I. BabuSka and B. Chayapathy, Stress computations for nearly incompressible materials by the p-version of the finite element method. Int. J. Numer. Methods Engng 28, 2175-2190 (1989). 13. I. BabuSka and D. Yu, Asymptotically exact a postriori error estimator for biquadratic elements. Finite Elements Anal. Des. 3, 341-354 (1987).