1、 Chapter 7 Error Analysis There are various sources of errors in finite element analysis. Error analysis has become indispensable in todays engineering application of finite element software. In this chapter, we will discuss the mature basis of error analysis for linear static analysis using finite
2、elements. 1 Sources of error The error of finite element results may arise at any procedures of finite element analysis. A good and mature user of finite element method should be aware of the all possible error sources in the finite element analysis. Generally, possible major sources of error can be
3、 put into the following categories 1) Modeling error It refers to the difference between a physical system and its mathematical model. It should be kept at an acceptable degree. Otherwise, the mathematical model need be revised. Usually, simplifications are introduced into the mathematical model. In
4、 the model, for instance, probably small holes and other geometric irregularities in the structure are neglected; loads are simplified; boundary conditions are idealized as by treating a component with relatively large stiffness as rigid. Quite often a two-dimensional problem is studied, ignoring it
5、s three-dimensional characteristic; a static analysis is performed, ignoring the dynamic character. 2) Discretization error The mathematical model has infinite number of dof, but finite dof are used in finite element analysis. The finite element solution is affected by the number of elements in the
6、model, the dof per node, the numerical integration rule, and etc. Actually, this kind of error has been discussed at the first section of last chapter. 3) Truncation error or rounding error It refers to loss of information due to truncation or rounding of numbers to fit a finite computer word length
7、. 4) Accumulation error This kind of error arises when global equations are solved repeatedly as in nonlinear or dynamic problems. The term numerical error is the combined result of truncation (or rounding) error and accumulation error. In this chapter, emphasis is mainly placed on discretization er
8、ror and numerical error. The above classification excludes possible user error and bugs. User error refers to mistakes made users after an effective model has been established. Bugs are programming errors embedded in finite element analysis software. Some of them are debugged during the preprocessin
9、g of finite - - 223element analysis. Dangerous bugs are those that affect the accuracy of finite element solution in an unobvious or unnoticeable manner. Although error is inevitable in finite element analysis, it cannot be regarded as excuse for gratuitous (unnecessary) errors. 2 Ill-conditioning T
10、he term “ill-conditioned” is often used in discussions of numerical error. As applied to a matrix, the term means that small numerical changes in matrix coefficients produce large changes in computed results. A coefficient matrix that is ill-conditioned for equation-solving is almost singular. This
11、can be illustrated in the following simple example. P E 1 A 1 , l 1 E 2 A 2 , l 2 x, u 1 2 Fig. 7.1 Non-prismatic bar under axial tension Denote 11 1 1 EA k l = , 22 2 2 EA k l = (7.1) The two segments are represented by two-node bar elements. The global equilibrium equation is given as 122 1 222 0
12、kk k u kku P + = or 12 2 2 1 22 2 1 () kk k u u kk uP u 0 + = +=(7.2a)(7.2b) The solutions are shown in Fig. 7.2. In the figure, the solid lines represent the exact description of the elements and the shaded bands indicate the numerical representation of the elements taking into account a finite num
13、ber of digits in computer memory. The solution is obtained at the intersection of the two lines. The intersection region is very small for well-conditioned situation if while the region is much larger in the ill-conditioned situation if . In the ill-conditioned situation, small change of can lead to
14、 significant changes in the solution of the system. 1 kk 2 1 k 2 k 1 kConsider an elimination solution of the system in (7.2). The displacement at node 1 is determined by 1 12 () P u kkk = 2 + (7.3) In the case of , the denominator approaches zero and significant error is produced if a 2 k 1 k - - 2
15、24small change of occurs. Physically, this implies that left segment has virtually no constraints on the right segment. An alert should be sounded in the solution process, indicating that the coefficient matrix is almost singular. 1 kIntersection region (6.2b) (6.2a) u 1 u 2(a) Well-conditioned ( )
16、12 kk Intersection region (6.2b) (6.2a) u 1 u 2(b) Ill-conditioned ( ) 21 kk Fig. 7.2 Example for matrix condition According to the classification in the last section, this is ascribed to truncation error. Mathematically, it has to do with the condition number of the coefficient matrix. The ill-cond
17、itioned situation in the above example can be summarized as the case where a stiff member is connected to a much more flexible member. A further plane frame example is given in Fig. 7.3. The axial stiffness of the horizontal member is much larger than the sideway stiffness of the two columns. As a r
18、esult, the computed values of the horizontal displacements at the two rigid joints may lack sufficient accuracy since they both are small quantities. It would be even worse if we use the difference of the two small quantities to compute the axial stress in the horizontal member, i.e., - - 225BC C B
19、BC () / Eu u L = . Actually, this stress is much smaller than the stresses due to bending. It may simply be taken as zero so that the ill-conditioned source is eliminated. P A B C DFig. 7.3 A possible ill-conditioned plane frame A way to circumvent ill-conditioned coefficient matrix is to choose dif
20、ferent dof. For instance, we may choose the relative dof 21 r uuu = (7.4) in the case shown in Fig. 7.1. Thus, eqns. (7.2) are changed into 1 1 2 0 0 r k P u ku P = (7.5) which is well-conditioned for any values of and . Other examples have already been mentioned in the construction of hierarchical
21、elements, since the added hierarchical dof achieve the improvement of solution accuracy without large increase of condition number of coefficient matrix. 1 k 2 k 3 Computation of condition number The condition number of a matrix is defined by max min = (7.6) where max and min are the largest and sma
22、llest eigenvalues of the matrix. It can be used as an index for prediction of possible loss of accuracy during solution process. The actual loss of accuracy during solution process is not as great as the estimate based on the condition number. Therefore, exact computation of condition number of a co
23、efficient matrix is not necessary. In addition, the result in the above expression sometimes may be misleading. For instance, the eigenvalues of the coefficient matrix of (7.2) can be found as - - 22622 121 max 24 2 kkkk + = 2 , 22 121 min 24 2 kkkk + = 2 . (7.7) According to the definition, the con
24、dition number is 22 22 12121212 22 12 1212 12 1 2 21 1 2 24( 24 4 24 ,;4, . kkkkkkkk kk kkkk kk k k kk k k + + = + 2 ) 2 1 k(7.8) It is large for both the case and the case , yet only the latter leads to significant error during equation-solving process. To identify the truly ill-conditioned situati
25、on, a scaling technique can be used. With a scaling matrix which is a diagonal matrix given as 1 kk 2 k S 11 22 1 1 1 nn K K K = S (7.9) Then, the condition number is calculated using the scaled coefficient matrix S = KS K S (7.10) with unity diagonal entries rather than the original coefficient mat
26、rix . The scaled matrix for the one in (7.2) is K 12 12 122 22 22 212 212 1010 01 01 11 S kk kk kk k kk kk kkk kkk + + = + = + K(7.11) The eigenvalues of the scaled matrix are found as 2 max 12 1 S k kk =+ + , 2 min 12 1 S k kk = + . (7.12) The condition number of the scaled matrix is obtained as ()
27、 2 22 12 2 12 12 11 S kk kk kk kk kk =+ = + + 1 ,. (7.13) which is large when the matrix in (7.2) is truly ill-conditioned. Suppose that digits are used to represent each number. The truncation error in the solution process of equations can only maintain the accuracy of coefficients to digits. With
28、the condition number obtained from the scaled matrix, a bound of the accuracy of coefficients matrix can be given as d acc d acc loss ddd = , with 10 log loss S d (7.14) For instance, if 14 d = and , then . A reason that may 8 10 S = 6 acc d 10 log S - - 227overestimate the accuracy loss is the effe
29、ct of force vector in the solution process of . The loss of accuracy of order is only materialized in special form of . f Kd = f 10 log S fThe mesh parameters can also influence the condition number significantly. A formula has been proposed, 21 2/ max min m mn e h bn h(7.15) = where is a positive n
30、umber that is proportional to the expression of Poissons ratio b 1( 1 2) , are maximum and minimum node spacings in the mesh, is the number of elements in the finite element model, 2 is the differential equation order and n is the dimensionality of the finite element mesh. max , h min h e n mEqn. (7
31、.15) is useful in predicting the effect of mesh parameters on the condition number. For example, if the ratio is changed from 1/1 to 1/10 while other parameters are not changed, increases by a factor of 1000 in a beam problem ( max / h min h 2 m 4 = ). If only the number of elements is doubled, incr
32、eases by a factor of 17. 4 Equation residuals The solution of global equations Kd = f (7.16) may not satisfy the equations exactly due to numerical error. The vector of residuals are obtained as R=f-K d (7.17) which should be a zero vector if no numerical error is present. Obviously, the vector of r
33、esiduals is a measure of error. A scalar form of the error measure is defined as = T T dR df(7.18) e Unfortunately, this error measure may not serve the purpose well since residuals tend to be small whether solution vector is accurate or not when direct solution method is used. Worse of all, it is p
34、ossible for the less accurate of two solution vectors to have the smaller residuals. Therefore, a small is not a reliable indicator of good accuracy of solution. However, a large is a reliable indicator that something goes wrong in the solution process. d R RTo minimize numerical error in the soluti
35、on of equations, one may adopt iterative solution technique. If a set of equations is seriously ill-conditioned, it is better to build a new finite element model than to continue improving a poor solution. - - 228 5 Mesh refinement and extrapolation Assume that is a quantity obtained from a finite e
36、lement solution. It can be a displacement component or a stress component. Let 1 ( m Oh) +be the order of error of , where is a characteristic dimension of element and may be regarded as the order of the element. A plot of h m versus 1 m h +is a straight line as shown in Fig. 7.4. Let 1 and 2 be val
37、ues of the quantity for two typical meshes and , respectively. A linear extrapolation gives 1 h 2 h 11 12 21 11 21 mm mm hh hh + + + = (7.19) which corresponds to the limit case as . The above equation is known as Richardsons extrapolation formula. 0 h 1 m h + 2 1 1 1 m h + 1 2 m h + Fig. 7.4 Richar
38、dsons extrapolation The application of Richardsons extrapolation formula is subject to several restrictions. In each refinement, nodes and interelement boundaries of the coarser mesh are preserved while new nodes, elements, and interelement boundaries are added. Corner nodes remain corner nodes and
39、side nodes remain side nodes as shown in Fig. 7.5. Fig. 7.5 Regular mesh refinement - - 229 When incompatible elements are used, the convergence is not monotonic. Besides, the value of is affected by distortion of mesh. Therefore, it is usually not appropriate to take as the exact value of . But it
40、can be used to get an estimate of the percentage error of a finer mesh, i.e., 2 100% e = (7.20)The refinement of mesh usually is put into the following three categories: 1) h refinement An h refinement is the increase of finite element mesh density by adding the same type elements to the mesh. 2) p
41、refinement A p refinement is referred to the increase the polynomial order in one element of the mesh without changing the mesh density. 3) r refinement An r refinement is referred to the relocation of some of the nodes in the finite element mesh without changing the number of nodes and element type
42、. A pictorial illustration of the three basic refinement types is given in Fig. 7.7. As dof are increased, p refinement has a high convergence rate than h refinement. Also, in hierarchical form a p refinement produces a global stiffness matrix of lower condition number than an h refinement having th
43、e same number of dof. The r refinement allows only limited improvement because the number of dof is unchanged. Although, optimization of node location is possible, the effort is expensive and rarely worthwhile. Of fair interest is the combination of the above three basic types of refinement in pract
44、ice. It is referred to as hr refinement when h refinement is accompanied with relocation of nodes; hp refinement is the case when the increase of mesh density with the same type of element and the addition of higher polynomial terms are concurrent. The hr and hp type refinements are relatively effec
45、tive in improving convergence as dof are added to a coarse mesh. - - 230y,v x,u (a) (b) (c) (d) P Fig. 7.6 Three basic refinement types: (a) initial mesh; (b) A possible h refinement; (c) A possible p refinement; (d) A possible r refinement. 6 Patch recovery technique It has discussed in last chapte
46、r that stresses computed at Gauss points are more accurate than those computed directly from nodal dof () e =EBd (7.21) The extrapolation of stresses from Gauss points is the improvement of stress accuracy at an element level. At the boundary nodes of the element, nodal averaging technique is used t
47、o smoothen the discontinuous stresses between elements. In this section, a technique is introduced to improve stress computation. First, a patch, which consists of a number of elements, is chosen. The patch may be one of the following three types (see Fig. 7.7): (1) An adjacent group of elements sur
48、rounding a particular node; (2) All elements adjacent to a particular element; (3) The adjacent elements that form a surface or a side - - 231 (a) A node surrounded by elements (b) An element surrounded by elements Boundary of domain(c) Elements along a portion of boundary Fig. 7.7 Three types of patches used in patch recovery of stresses A stress component over the patch is expressed in terms of some generalized parameters, i.e., * = P (7.22) where * is any stress component, contains terms of a polynomial, and P - - 232contains generalized parameters to
