- Email: [email protected]
Numerical stability of complex finite elements
NUMERICAL STABILITY OF COMPLEX FINITE ELEMENTS F. Mmut Department of Civil Engineering,UniversityCollege, Cardii, Wales
and J. MARSHALQ ~pa~rnent of Civil Engineers, Universityof Strathclyde,Glasgow G4 ONG, !kotland (Received 4 June 1975) Abstract-The use of complex finite elements has been shown to give considerableeconomies in computingtime but there is an inherentdangerthat numericalinstabilitycan arise duringthe formulationof the.stillness matrixfor such elements. This problemis studiedwith particularreferenceto rectangularplane stress elementsfor which any num~r of nodes along a side may be prescribed.It is shown that the occurrenceof numerical~s~~~ is directly related to the degree of the polynomialassumed for the displacementfunction and to the numberof sign&ant figuresused in the calculation.The particularmatrixoperationswhich are susceptibleto errorare indicated.
order polynomialscan lead to numerical instabilitiesin some states of the formulationof the stiffness matrix of matrixof coefficientsof poI~omi~s an element. The occurrence of numerical instability is matrixrelatingstrainand ~lynomial coet?knts illustrated here with reference to rectanguh~ elements matrix relating nodal displacements and polynomial for plane stress analysis. Each node on the element is coefficients restricted to two degrees of freedom but the complexity matrixrelatingstress and strain is increased by increasing the number of nodes on the Young’smodulus element. The method used for derivinga stiffness matrix stiffness matrix number of nodes along one side of an element in the is that proposed by Mishuf41.The choice of displacement ~ly~rni~ and subs~~nt calculations for an x-direction number of nodes atong one side of an element in the element are completely systema~ and depend only on y-diiection the number of nodes on the element. Thus the problem displacementfunction of numerical instability can be related directly to the nodal displacements assumed polynomial and to the elastic and geometric strainmatrix properties of the element. polynomialin x and y The computer program for the finite element analyses subscript referring to component arising from normal using the complex elements was run on three different stress and strain computers-XC1 1904, ICL 470 and CDC 7600. The stress matrix subscript referring to component arising from shear implementation of the program on each of these machinesrequired only minor modificationsin the use of stresses and strains Poisson’s ratio backing store facilities and thus all numericalprocesses were identicalregardlessof t&ecomputer.The significant diierence in these computers is in the numerical NOMRNCUTURE
A tt D ; NX NY u d
i
n
4 S Ir
accuracy with which astir
The continuing study of the finite element method has led to the production of increasingly complex elements. The increase in complexity can take the form of an increase in the number of degrees-of-freedom at a
node[11, an increase in the number of nodes on an eIementf21or a combinationof both these methods.The advantage of these procedures is that the displacement function for the element is capable of givinghigher order deformations over the element than for more simple element types thus fewer elements are required for an analysis. It has been shown[3] that such elements can give economies in computingcosts. The increase in complexity can be achieved by using higher order ~lynomials to represent the ~splacement function. In theory there is no knit to the degree of complexity or the order of the polynomial but higher SPostdoctoral ResearchFellow. SSeniorLecturer.
c~c~ations are evalu-
ated. A range of complex elements with reliablestiffness properties is indicated for each of these computers together with an additionalrange of elements if steps are taken to improve the numerical stability of the mathematicaloperations involved. ~~~
~~~~
A full explanation of the method of deriving an element stiffness matrix is given by Mishu[4] but a brief outline is necessary here to show in particular the composition of terms in the stitlness matrix. In plane stress analysis two ~lyno~~s V. and U, are assumed to represent the ~sp~ement functions in the x and y directions respectively, hence
in which fx and f,. are the terms in the respective 259
F. MISHU and 11. MARSHALL
260
polynomials and A, and A, are the corresponding coefficients. The nodal displacements can be expressed as a function of the coefficients of the terms in the polynomials as
where d, and d, are the nodal displacements in the x and y directions respectively. A solution for the coefficients gives [Al = [ ;]
= [‘ii
;_I][
$1=[C-‘lfd].
Zienkiewiczl21 from the expression
WI=
[C-‘I’ 1 j
WTIWIWIIdv,,
VI.
(6)
The terms in the stiffness matrix may be decomposed into two components K, and KS which represent the contributions from normal and shear stresses respectively. Further the components can be partitioned to correspond to x and y direction stiffnesses as
(3)
The strains can be determined in terms of the nodal displacements by differentiation of eqn (1) and substitution of eqn (3) for the constants to give
.
(7) The following two expressions illustrate the typical form of each of the submatrices in eqn (7).
NuMEIucAL EXAMPLX
The general term B, represents the terms resulting from the partial differentiation of the polynomial f, with respect to direction j. The relationship between the stresses fq] and the strains [e] for plane stress is
(5) From the above equations the sti&ess matrix can be evaluated in the standard manner described by
It is custom~y to test the behaviour of particular finite element on a simple problem for which an accurate solution exists. In this case the problem of a cantilever beam, subjected to a uniformly distributed load of unit intensity, with the dimensions shown in Fig. 1 was adop ted. In all analyses Young’s modulus and Poisson’s ratio were taken as 1.0 and 0.0 respectively. The exact solution to this problem is assumed to be that given by simple beam theory. A slight difhculty arises with the representation of the uniformly distributed load in the finite element analysis which allows only two degrees of freedom at a node. The distributed load has to be applied as a system of concentrated loads at the nodes. Three possible systems are shown in Fig. 1. The results for vertical deflection at “a” and bending stress at “b” are given in Table I. This shows that a reasonable accuracy can be obtained with
w = I Per unit length b
25
2.5
2.5
125
(d)
Fig. 1. Load approximations on cantilever beam.
Numericals~b~ity of
compiex
four loads as shown in Fig. l(d). It should be noted that the variation of deflection along the beam for a uniformly distributed load is quark whereas the equivalent system of concentrated loads leads to a series of discontinuous cubic variations.
The simplicity of the cantilever beam permits a solution using only one element. A possible element possessing five nodes in the x-direction and two nodes in the y-direction is shown in Fig. 2. With this choice of element sufficient flexibility is provided to produce a quartic variation in displacement regardless of the dis~bution of the nodal concentrated loads. The element was analysed with the three load cases shown in Fig. 2. These loadings are identical to those shown in Fig. 1. The solutions, using the ICL 1904 computer, for the deflections at “a” and the stresses at “b” for each case are given in Table 1. A comparison of the results from beam theory and the one element solution show a close agreement. Improvements in the deflections occur as the load is distributed over a larger number of nodes. It could be expected that an increase in the number of nodes on the element in the x-direction (NX) would lead to an increase in accuracy due to the improved dis~bution of concentrated loads which would be possible. The tinite element results also indicate that two nodes in the ydirection (NY) are sufficient for this type of problem. Further analyses were carried out in which NY was held constant at two and NX was varied from four to eleven nodes. In each case concentrated loads were ~pprop~~ely dist~buted to all nodes on the upper surface. The variation in deflections and stresses are plotted against the number of nodes (NX) as shown in Fig, 3 for the case of single precision arithmetic. It can be seen that when NX exceeds seven instability takes place in both deflections and stresses. It is
finite elements
300.--
A.
261
j
-:--“-
*
Beam theory --
* 100x 12 16 Numbzr of nod8esinx - direction kVXf
20
o^ ‘0 ;i
16-
Y
Jq,,
p
15:
B r
l4-
P p
13
E 0
I2
.-.
.-
w
II
. 0
4
8
Beam theory
,x I2
I6
1 20
Number of nodes in x direction (NX
Fig. 3. Solutions using ICL PX.KIcomputer. 0, Single precision arithmetic; x, double precision arithmetic.
thus clear that two regions exist in the present system. The first being a stable region for which NX is less than or equal to seven, and the second being the unstable region for which NX exceeds seven. SOUSCES OF ECU ~~~~ The causes of the ins~b~ity observed may be studied in terms of the parameters involved in the formulation of the element stiffness matrix. There are three possible sources of numerical instability and these are examined in detail.
Table t. Effects of concentrated loads representing a uniformly distributed load Number of points loads
x
Beam theory Stress at “b” DeR’nat “a”
One element Stress at “b” Detl’nat “a”
2.00x IlV
300
2.01x IO’
300
I.51x 104
300
154x 10’
I.50 x IO’
300
303 -
-
cdl
(c) Fig. 2. Load approximations on element.
F. MISHU and J.
262
properties The material properties include Poisson’s ratio and Young’s modulus. The effect of these factors may be studied when the element stiffness matrix is decomposed into the two components arising from normal and shear stresses as shown in eqn (8). It can be seen from eqn (8) that the only effect of Poisson’s ratio is to define a constant multiplying factor which is applied to all terms within a component matrix. The ratios of this factor for shearing and normal components can only vary between 2 : 1 and 3 : 1 approximately, depending on the value of Poisson’s ratio. In this way Poisson’s ratio may be viewed as an external effect. Young’s modulus is external not only to the component stiffness matrices but to the total stiffness matrix. Therefore Poisson’s ratio and Young’s modulus cannot account for the numerical instability of the stiffness matrix. Material
Number of nodes on the element The number of nodes on the element in the x-direction (NX) defines the degree of the polynomial in x which forms part of the total displacement function. The terms of the displacement function are present in all the matrices of eqn (8). To study the effect of the number of nodes on the stiffness matrix implies the study of each of the matrices present in the formulation. By considering eqn (8) which represents the typical operations involved in forming the stiffness matrix, the principal operations involved are summarized by: 1. The formulation and inversion of the matrix [Cl. 2. The evaluation of j[BIT - [ItId,,, which results in a matrix to be denoted as [BB]. 3. The evaluation of the product [C-‘lTIBB][C-‘I. An instability in any of the above matrices or in the operations involved will affect the stability of the stiffness matrix. A very significant influence must have taken place in some or all of these factors to cause the change in the variation of deflection and stress with respect ot NX in the unstable range. Formulation and inverse of the matrix [C] The terms in the polynomial for fx or f, with NX = 11 and NY=2are
MARSHALL
evaluate the product [C-l] x [Cl and compare the result with the unit matrix. This has the advantage that the errors in multiplying two matrices will be small in comparison with those of inversion. This was the procedure which was adopted. An illustration of the result of the product [C-‘1 x [Cl for the case of NX = 11, as given by the ICL 1904 computer using single precision arithmetic is shown in Table 2. This table lists the terms in the fifth row of the matrix. It is obvious that some of the offdiagonal terms are subject to large errors and can have magnitudes greater than the leading diagonal. However in all rows the leading diagonal term was very close to unity. It is inevitable that the evaluation of the terms in the polynomial of eqn (9) will lead at some nodes to large numerical differences. During the inversion process some small numbers are multiplied by a large factor and then subtracted from a number of comparable magnitudes. This can result in a term relying on some of the least significant figures in the calculation and hence becomes inaccurate. On the ICL 1904 computer real numbers are held with eleven significant figures but the difference in magnitudes of numbers in the [C] matrix can be of the order of 10” for NX = 11. This can be explained in the engineering sense in that the occurrence of a very large number in a diagonal can be equivalent to the elimination of that equation. The inversion process was improved by reverting to double precision arithmetic giving twenty significant figures. The values of the terms in the fifth row for the product [C-l] x [Cl are again given in Table 2. This shows a significant improvement. The maximum value of any off-diagonal term in the matrix was of the order of 1O-9and the product is for all practical purposes a unit matrix. Evaluation of IBTBd,,, By way of illustration the particular decomposed submatrix
as shown in eqn (8) will be used. The derivative of the polynomial of eqn (9) to yield B,, followed by the product Bz. B, and subsequent integration leads to terms of the following order:
fx = f(1, y, x, xy, . . .)X8,x*y, x9, x9y, X’O,x’Oy). (9) In general, if NY = 2, then the polynomial contains all terms below those in which x is raised to the power NX There is slight evidence from Fig. 3 of numerical instability at NX = 8 but it is very obvious for NX > 8. In the formulation of [C] the polynomial has to be evaluated at the coordinates of each node. The change from NX= 8 to NX= 9 involves only the additional terms x8 and x*y. The raising of an x coordinate by an extra power does not suggest any complication as the evaluation of the matrix [BB] for NX = 7 contains powers higher than x8 and no instability is evident in that solution. It is apparent that any instability must lie in the inverse of the matrix [Cl. One method of checking the stability of an inversion process is to invert the inverse and compare the result with the original. This greatly magnifies any error and makes it difficult to judge the degree of error in the initial inversion. An alternative is to
x”Y, x”Y*, x19y, x19yZ.
These terms are then evaluated over the area of the element. It is obvious that there are no problems as far Table 2. Evaluationof [C-l] . [C] on Column 1-2 3-4 5-6 7-8 9-10 11-12 13-14 15-16 17-18 19-20 21-22
Single precision 2.074E10 2.820E-8 l.OOOEO 1.211EA 1.004E5 7.129E-5 6.790E-4 5.676E3 5.054E2 3.867El 3.5OOE0
O.OOOEO
ICL 1904 computer
Double precision 1.606E17
O.OOOEO3.815E16 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OMlE0 O.OOOE0 O.OOOE0 O.OOOE0
l.OOOEO 5.866&13 3.206E11 2.008E-9 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0
7.497Ek17 2.165-E-15 8.499~14 4.263E-12 2.533E-10
O.OOOE 0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0
Numerical stability of complex finite elements
as the y-limits are concerned, but the x-limits raised to such powers are subject to rounding-off errors. Each term is liable to various degrees of error which are a function of the power they are raised to. A computer printout of the matrix [BB] revealed that symmetry was obtained in all cases despite errors. This is because the operations are carried out in a specific type of operation. Since the terms occur in exactly the same manner above and below the leading diagonal a symmetric matrix must result. The accuracy of the results obtained depends on the number of significant figures to which any number is held and operated upon. A comparison of the results for the matrix [BB] using single and double precision arithmetic showed that only the eleventh significant figure was altered. It therefore appears that this operation is reasonably stable. Product [C-‘I’ . [BB] . [C-‘1 The improvement on the operations of the inverse of [Cl and the formation of [BB] alone did not lead to a stiffness matrix which was symmetrical. Since the matrix [BB] is symmetrical at all times regardless of the accuracy of the terms, then the above product should lead to a symmetrical stiffness matrix. The product must be evaluated in two parts. As both the matrices [C-‘I and [BB] contain terms with large differences in magnitude the accumulation of the products of these terms will lead to rounding-off errors in the first sub-product [C-‘lT - [BB]. The subsequent post-multiplication by [C-t] will then give an increase over the initial errors. Again it was found that the use of double precision arithmetic led to symmetric stiffness matrices over an increased range of NX. INFLUENCE OF
ARRHMETICAL.
263 Beam theory .-
.-.
X
Number
of nodes in x direction
Number Fig. 4.
of nodes
(NX)
in x dIrectIon
Solutions using ICL 470 computer. 0, Single precision arithmetic; X, double precision arithmetic.
PRFAXION
To some extent the numerical accuracy on a computer can be selected but the degree of selection depends on both the computer and the language used. On the ICL 470 computer real constants are normally held within seven significant figures while double precision constants are held within fifteen significant figures. The significant figures for the ICL 1904 computer for single and double precision are eleven and twenty respectively. On the CDC 7600 the choice lies between fifteen and twenty-nine significant figures. The cantilever beam was solved on each of these computers using both single and double precision arithmetic for the operations which have been studied. Comparisons of the solutions for each computer are shown in Figs. 3-5. A summary of the stable regions is given in Table 3. In the case of the solution on the CDC 7600 using double precision arithmetic a maximum value of NX = 20 was tried, as which point there was no evidence of numerical instability. ASPECT RATIO OF AN ELEMENT
So far the stability of the stiffness properties have been considered with respect to the material properties and the number of nodes on the element. The geometric aspect of the element is of interest and a study was conducted on a cantilever beam subjected to a point load at the free end as shown in Fig. 6. Again a one element solution was used but the number of nodes on the element was kept constant at NX = 5 and NY = 3. The length of the beam was varied from 5 to 80 units. The maximum bending stresses at the fixed end. from solutions using the ICL 1904 computer with single pre-
C
I
4
Number
II
8
I
,
12
I6
of nodes in x direction
’ 4
Number
’
I 6
of nodes
’ 12
’
20
(NH
’ 16
’
’ 20
in x direction
Fig. 5. Solutions using CDC 7600 computer. 0, Single Precision arithmetic;X, double precision arithmetic.
cision arithmetic, are compared with simple beam theory in Fig. 6. Only stresses have beem shown in this case as these are sufficient to indicate any numerical instability. Figure 6 shows that the one-element solution for stress deviates slightly from the theory for span to depth ratios of more than 50, but the maximum error is only of the order of 4%. This gradually changing error is not characteristic of numerical instability in the stiffness matrix. The error is most likely caused by an inadequate displacement function in an extreme situation. It is assumed that the aspect ratio of the element is not a significant contributory factor causing numerical instability.
F. MISHUand J.
264
MARSHALL
Table 3. Maximum values of NX for stable stiffness matrices
Computer ICL 1900 ICL 470 600
CDC 7600
1
500
t
z
400
P G
I
t
2 g 300 v) / 200 L#llj
/ /
100
/ 0
0
IO
20
30
40
50
60
70
80
Span to depth ratio
Fig. 6. Stresses from variation in aspect ratio of cantilever beam. -.-+, beam theory; 0, one element solution.
CONCLUSIONS
The problems studied have been restricted to a one element solution. This has the advantage that the structural stiffness matrix and the element stiffness matrix are identical. Therefore, any numerical instability in the solution is directly a function of the element stiffness matrix. This is not to be confused with the problem of numerical instability associated with a large structural stiffness matrix assembled from many element stiffness matrices in which the individual stiffness matrices may be numerically stable. No special techniques were used to combat numerical instability apart from using double precision arithmetic. The method of inversion was by Gaussian elimination and all integration was carried out explicitly. Table 3 shows clearly that the stability of the stiffness matix is directly related to the numerical accuracy of the operations involved in the formulation. In particular there are methods of inversion available which will provide greater numerical accuracy than Gaussian elimination for a given set of significant figures. It is not known to what extent the range of stability could be increased by such methods in comparision with the use of increased arithmetical precision. The explicit integration over the area of the element produced no significant instability. However, if a numerical integration is used, then it is
Significant figures
Maximum NX
II 20 1 I5 15 29
I5 5 I2 II >20
1
possible that instability could arise at this stage of the calculation. The method of deriving stiffness matrices for complex elements in plane stress analyses which has been used is particularly suitable for the study of numerical stability in general because of the systematic nature of choosing an appropriate polynomial for the displacement function. From the examples given, it is clear that the numerical stability is a function of the degree of polynomial assumed. Had a study been made of a random selection of complex elements then the results could be confused by the use of special techniques to eliminate or add terms in a chosen polynomial. Although these techniques could influence the numerical stability, the results given do indicate the relationship between the degree of polynomial and the numerical accuracy which defines a numerically stable region. It seems reasonable to assume that these relationships are approximately true regardless of the purpose of a finite element. The study has been conducted using the C-inverse techniques but further investigations would be required to establish the range of stable elements when the unit shape function technique is used[5]. It appears likely that the range of elements produced using this approach will be larger because no C-inverse is required, hence eliminating a source of error. The advantages here may be outweighed by the increased computing time and storage requirements although a similar disadvantage also applies when resort has to be made to double precision arithmetic. REFERENCES
I. 1.Holland, The Finire Element Method in Plane Stress Analysis. Finite Element Methods in Stress Analysis, Tapir, Trondheim, Norway (1%9). 2. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill. London (19711. 3. J. Marshall and F. Mishu, An investigation into the economics of using complex elements in plane stress analysis. Inr. Conj. on Finite Element Methods in Eng. Sci. The University of New South Wales, Australia (1974). 4. F. Mishu, Aspects of complex rectangular elements in plane stress analysis. Ph.D. Thesis, University of Strathclyde, Glasgow (1974). 5. 0. C. Zienkiewicz, B. M. Irons, J. Ergatoudis, S. Ahmad and F. C. Scott, Iso-Parametric and Associated Element Families for Two- and Three-Dimensional Analysis. Finite Element Methods in Stress Analysis, Tapir, Trondheim, Norway (1%9).
and J. MARSHALQ ~pa~rnent of Civil Engineers, Universityof Strathclyde,Glasgow G4 ONG, !kotland (Received 4 June 1975) Abstract-The use of complex finite elements has been shown to give considerableeconomies in computingtime but there is an inherentdangerthat numericalinstabilitycan arise duringthe formulationof the.stillness matrixfor such elements. This problemis studiedwith particularreferenceto rectangularplane stress elementsfor which any num~r of nodes along a side may be prescribed.It is shown that the occurrenceof numerical~s~~~ is directly related to the degree of the polynomialassumed for the displacementfunction and to the numberof sign&ant figuresused in the calculation.The particularmatrixoperationswhich are susceptibleto errorare indicated.
order polynomialscan lead to numerical instabilitiesin some states of the formulationof the stiffness matrix of matrixof coefficientsof poI~omi~s an element. The occurrence of numerical instability is matrixrelatingstrainand ~lynomial coet?knts illustrated here with reference to rectanguh~ elements matrix relating nodal displacements and polynomial for plane stress analysis. Each node on the element is coefficients restricted to two degrees of freedom but the complexity matrixrelatingstress and strain is increased by increasing the number of nodes on the Young’smodulus element. The method used for derivinga stiffness matrix stiffness matrix number of nodes along one side of an element in the is that proposed by Mishuf41.The choice of displacement ~ly~rni~ and subs~~nt calculations for an x-direction number of nodes atong one side of an element in the element are completely systema~ and depend only on y-diiection the number of nodes on the element. Thus the problem displacementfunction of numerical instability can be related directly to the nodal displacements assumed polynomial and to the elastic and geometric strainmatrix properties of the element. polynomialin x and y The computer program for the finite element analyses subscript referring to component arising from normal using the complex elements was run on three different stress and strain computers-XC1 1904, ICL 470 and CDC 7600. The stress matrix subscript referring to component arising from shear implementation of the program on each of these machinesrequired only minor modificationsin the use of stresses and strains Poisson’s ratio backing store facilities and thus all numericalprocesses were identicalregardlessof t&ecomputer.The significant diierence in these computers is in the numerical NOMRNCUTURE
A tt D ; NX NY u d
i
n
4 S Ir
accuracy with which astir
The continuing study of the finite element method has led to the production of increasingly complex elements. The increase in complexity can take the form of an increase in the number of degrees-of-freedom at a
node[11, an increase in the number of nodes on an eIementf21or a combinationof both these methods.The advantage of these procedures is that the displacement function for the element is capable of givinghigher order deformations over the element than for more simple element types thus fewer elements are required for an analysis. It has been shown[3] that such elements can give economies in computingcosts. The increase in complexity can be achieved by using higher order ~lynomials to represent the ~splacement function. In theory there is no knit to the degree of complexity or the order of the polynomial but higher SPostdoctoral ResearchFellow. SSeniorLecturer.
c~c~ations are evalu-
ated. A range of complex elements with reliablestiffness properties is indicated for each of these computers together with an additionalrange of elements if steps are taken to improve the numerical stability of the mathematicaloperations involved. ~~~
~~~~
A full explanation of the method of deriving an element stiffness matrix is given by Mishu[4] but a brief outline is necessary here to show in particular the composition of terms in the stitlness matrix. In plane stress analysis two ~lyno~~s V. and U, are assumed to represent the ~sp~ement functions in the x and y directions respectively, hence
in which fx and f,. are the terms in the respective 259
F. MISHU and 11. MARSHALL
260
polynomials and A, and A, are the corresponding coefficients. The nodal displacements can be expressed as a function of the coefficients of the terms in the polynomials as
where d, and d, are the nodal displacements in the x and y directions respectively. A solution for the coefficients gives [Al = [ ;]
= [‘ii
;_I][
$1=[C-‘lfd].
Zienkiewiczl21 from the expression
WI=
[C-‘I’ 1 j
WTIWIWIIdv,,
VI.
(6)
The terms in the stiffness matrix may be decomposed into two components K, and KS which represent the contributions from normal and shear stresses respectively. Further the components can be partitioned to correspond to x and y direction stiffnesses as
(3)
The strains can be determined in terms of the nodal displacements by differentiation of eqn (1) and substitution of eqn (3) for the constants to give
.
(7) The following two expressions illustrate the typical form of each of the submatrices in eqn (7).
NuMEIucAL EXAMPLX
The general term B, represents the terms resulting from the partial differentiation of the polynomial f, with respect to direction j. The relationship between the stresses fq] and the strains [e] for plane stress is
(5) From the above equations the sti&ess matrix can be evaluated in the standard manner described by
It is custom~y to test the behaviour of particular finite element on a simple problem for which an accurate solution exists. In this case the problem of a cantilever beam, subjected to a uniformly distributed load of unit intensity, with the dimensions shown in Fig. 1 was adop ted. In all analyses Young’s modulus and Poisson’s ratio were taken as 1.0 and 0.0 respectively. The exact solution to this problem is assumed to be that given by simple beam theory. A slight difhculty arises with the representation of the uniformly distributed load in the finite element analysis which allows only two degrees of freedom at a node. The distributed load has to be applied as a system of concentrated loads at the nodes. Three possible systems are shown in Fig. 1. The results for vertical deflection at “a” and bending stress at “b” are given in Table I. This shows that a reasonable accuracy can be obtained with
w = I Per unit length b
25
2.5
2.5
125
(d)
Fig. 1. Load approximations on cantilever beam.
Numericals~b~ity of
compiex
four loads as shown in Fig. l(d). It should be noted that the variation of deflection along the beam for a uniformly distributed load is quark whereas the equivalent system of concentrated loads leads to a series of discontinuous cubic variations.
The simplicity of the cantilever beam permits a solution using only one element. A possible element possessing five nodes in the x-direction and two nodes in the y-direction is shown in Fig. 2. With this choice of element sufficient flexibility is provided to produce a quartic variation in displacement regardless of the dis~bution of the nodal concentrated loads. The element was analysed with the three load cases shown in Fig. 2. These loadings are identical to those shown in Fig. 1. The solutions, using the ICL 1904 computer, for the deflections at “a” and the stresses at “b” for each case are given in Table 1. A comparison of the results from beam theory and the one element solution show a close agreement. Improvements in the deflections occur as the load is distributed over a larger number of nodes. It could be expected that an increase in the number of nodes on the element in the x-direction (NX) would lead to an increase in accuracy due to the improved dis~bution of concentrated loads which would be possible. The tinite element results also indicate that two nodes in the ydirection (NY) are sufficient for this type of problem. Further analyses were carried out in which NY was held constant at two and NX was varied from four to eleven nodes. In each case concentrated loads were ~pprop~~ely dist~buted to all nodes on the upper surface. The variation in deflections and stresses are plotted against the number of nodes (NX) as shown in Fig, 3 for the case of single precision arithmetic. It can be seen that when NX exceeds seven instability takes place in both deflections and stresses. It is
finite elements
300.--
A.
261
j
-:--“-
*
Beam theory --
* 100x 12 16 Numbzr of nod8esinx - direction kVXf
20
o^ ‘0 ;i
16-
Y
Jq,,
p
15:
B r
l4-
P p
13
E 0
I2
.-.
.-
w
II
. 0
4
8
Beam theory
,x I2
I6
1 20
Number of nodes in x direction (NX
Fig. 3. Solutions using ICL PX.KIcomputer. 0, Single precision arithmetic; x, double precision arithmetic.
thus clear that two regions exist in the present system. The first being a stable region for which NX is less than or equal to seven, and the second being the unstable region for which NX exceeds seven. SOUSCES OF ECU ~~~~ The causes of the ins~b~ity observed may be studied in terms of the parameters involved in the formulation of the element stiffness matrix. There are three possible sources of numerical instability and these are examined in detail.
Table t. Effects of concentrated loads representing a uniformly distributed load Number of points loads
x
Beam theory Stress at “b” DeR’nat “a”
One element Stress at “b” Detl’nat “a”
2.00x IlV
300
2.01x IO’
300
I.51x 104
300
154x 10’
I.50 x IO’
300
303 -
-
cdl
(c) Fig. 2. Load approximations on element.
F. MISHU and J.
262
properties The material properties include Poisson’s ratio and Young’s modulus. The effect of these factors may be studied when the element stiffness matrix is decomposed into the two components arising from normal and shear stresses as shown in eqn (8). It can be seen from eqn (8) that the only effect of Poisson’s ratio is to define a constant multiplying factor which is applied to all terms within a component matrix. The ratios of this factor for shearing and normal components can only vary between 2 : 1 and 3 : 1 approximately, depending on the value of Poisson’s ratio. In this way Poisson’s ratio may be viewed as an external effect. Young’s modulus is external not only to the component stiffness matrices but to the total stiffness matrix. Therefore Poisson’s ratio and Young’s modulus cannot account for the numerical instability of the stiffness matrix. Material
Number of nodes on the element The number of nodes on the element in the x-direction (NX) defines the degree of the polynomial in x which forms part of the total displacement function. The terms of the displacement function are present in all the matrices of eqn (8). To study the effect of the number of nodes on the stiffness matrix implies the study of each of the matrices present in the formulation. By considering eqn (8) which represents the typical operations involved in forming the stiffness matrix, the principal operations involved are summarized by: 1. The formulation and inversion of the matrix [Cl. 2. The evaluation of j[BIT - [ItId,,, which results in a matrix to be denoted as [BB]. 3. The evaluation of the product [C-‘lTIBB][C-‘I. An instability in any of the above matrices or in the operations involved will affect the stability of the stiffness matrix. A very significant influence must have taken place in some or all of these factors to cause the change in the variation of deflection and stress with respect ot NX in the unstable range. Formulation and inverse of the matrix [C] The terms in the polynomial for fx or f, with NX = 11 and NY=2are
MARSHALL
evaluate the product [C-l] x [Cl and compare the result with the unit matrix. This has the advantage that the errors in multiplying two matrices will be small in comparison with those of inversion. This was the procedure which was adopted. An illustration of the result of the product [C-‘1 x [Cl for the case of NX = 11, as given by the ICL 1904 computer using single precision arithmetic is shown in Table 2. This table lists the terms in the fifth row of the matrix. It is obvious that some of the offdiagonal terms are subject to large errors and can have magnitudes greater than the leading diagonal. However in all rows the leading diagonal term was very close to unity. It is inevitable that the evaluation of the terms in the polynomial of eqn (9) will lead at some nodes to large numerical differences. During the inversion process some small numbers are multiplied by a large factor and then subtracted from a number of comparable magnitudes. This can result in a term relying on some of the least significant figures in the calculation and hence becomes inaccurate. On the ICL 1904 computer real numbers are held with eleven significant figures but the difference in magnitudes of numbers in the [C] matrix can be of the order of 10” for NX = 11. This can be explained in the engineering sense in that the occurrence of a very large number in a diagonal can be equivalent to the elimination of that equation. The inversion process was improved by reverting to double precision arithmetic giving twenty significant figures. The values of the terms in the fifth row for the product [C-l] x [Cl are again given in Table 2. This shows a significant improvement. The maximum value of any off-diagonal term in the matrix was of the order of 1O-9and the product is for all practical purposes a unit matrix. Evaluation of IBTBd,,, By way of illustration the particular decomposed submatrix
as shown in eqn (8) will be used. The derivative of the polynomial of eqn (9) to yield B,, followed by the product Bz. B, and subsequent integration leads to terms of the following order:
fx = f(1, y, x, xy, . . .)X8,x*y, x9, x9y, X’O,x’Oy). (9) In general, if NY = 2, then the polynomial contains all terms below those in which x is raised to the power NX There is slight evidence from Fig. 3 of numerical instability at NX = 8 but it is very obvious for NX > 8. In the formulation of [C] the polynomial has to be evaluated at the coordinates of each node. The change from NX= 8 to NX= 9 involves only the additional terms x8 and x*y. The raising of an x coordinate by an extra power does not suggest any complication as the evaluation of the matrix [BB] for NX = 7 contains powers higher than x8 and no instability is evident in that solution. It is apparent that any instability must lie in the inverse of the matrix [Cl. One method of checking the stability of an inversion process is to invert the inverse and compare the result with the original. This greatly magnifies any error and makes it difficult to judge the degree of error in the initial inversion. An alternative is to
x”Y, x”Y*, x19y, x19yZ.
These terms are then evaluated over the area of the element. It is obvious that there are no problems as far Table 2. Evaluationof [C-l] . [C] on Column 1-2 3-4 5-6 7-8 9-10 11-12 13-14 15-16 17-18 19-20 21-22
Single precision 2.074E10 2.820E-8 l.OOOEO 1.211EA 1.004E5 7.129E-5 6.790E-4 5.676E3 5.054E2 3.867El 3.5OOE0
O.OOOEO
ICL 1904 computer
Double precision 1.606E17
O.OOOEO3.815E16 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OMlE0 O.OOOE0 O.OOOE0 O.OOOE0
l.OOOEO 5.866&13 3.206E11 2.008E-9 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0
7.497Ek17 2.165-E-15 8.499~14 4.263E-12 2.533E-10
O.OOOE 0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0 O.OOOE0
Numerical stability of complex finite elements
as the y-limits are concerned, but the x-limits raised to such powers are subject to rounding-off errors. Each term is liable to various degrees of error which are a function of the power they are raised to. A computer printout of the matrix [BB] revealed that symmetry was obtained in all cases despite errors. This is because the operations are carried out in a specific type of operation. Since the terms occur in exactly the same manner above and below the leading diagonal a symmetric matrix must result. The accuracy of the results obtained depends on the number of significant figures to which any number is held and operated upon. A comparison of the results for the matrix [BB] using single and double precision arithmetic showed that only the eleventh significant figure was altered. It therefore appears that this operation is reasonably stable. Product [C-‘I’ . [BB] . [C-‘1 The improvement on the operations of the inverse of [Cl and the formation of [BB] alone did not lead to a stiffness matrix which was symmetrical. Since the matrix [BB] is symmetrical at all times regardless of the accuracy of the terms, then the above product should lead to a symmetrical stiffness matrix. The product must be evaluated in two parts. As both the matrices [C-‘I and [BB] contain terms with large differences in magnitude the accumulation of the products of these terms will lead to rounding-off errors in the first sub-product [C-‘lT - [BB]. The subsequent post-multiplication by [C-t] will then give an increase over the initial errors. Again it was found that the use of double precision arithmetic led to symmetric stiffness matrices over an increased range of NX. INFLUENCE OF
ARRHMETICAL.
263 Beam theory .-
.-.
X
Number
of nodes in x direction
Number Fig. 4.
of nodes
(NX)
in x dIrectIon
Solutions using ICL 470 computer. 0, Single precision arithmetic; X, double precision arithmetic.
PRFAXION
To some extent the numerical accuracy on a computer can be selected but the degree of selection depends on both the computer and the language used. On the ICL 470 computer real constants are normally held within seven significant figures while double precision constants are held within fifteen significant figures. The significant figures for the ICL 1904 computer for single and double precision are eleven and twenty respectively. On the CDC 7600 the choice lies between fifteen and twenty-nine significant figures. The cantilever beam was solved on each of these computers using both single and double precision arithmetic for the operations which have been studied. Comparisons of the solutions for each computer are shown in Figs. 3-5. A summary of the stable regions is given in Table 3. In the case of the solution on the CDC 7600 using double precision arithmetic a maximum value of NX = 20 was tried, as which point there was no evidence of numerical instability. ASPECT RATIO OF AN ELEMENT
So far the stability of the stiffness properties have been considered with respect to the material properties and the number of nodes on the element. The geometric aspect of the element is of interest and a study was conducted on a cantilever beam subjected to a point load at the free end as shown in Fig. 6. Again a one element solution was used but the number of nodes on the element was kept constant at NX = 5 and NY = 3. The length of the beam was varied from 5 to 80 units. The maximum bending stresses at the fixed end. from solutions using the ICL 1904 computer with single pre-
C
I
4
Number
II
8
I
,
12
I6
of nodes in x direction
’ 4
Number
’
I 6
of nodes
’ 12
’
20
(NH
’ 16
’
’ 20
in x direction
Fig. 5. Solutions using CDC 7600 computer. 0, Single Precision arithmetic;X, double precision arithmetic.
cision arithmetic, are compared with simple beam theory in Fig. 6. Only stresses have beem shown in this case as these are sufficient to indicate any numerical instability. Figure 6 shows that the one-element solution for stress deviates slightly from the theory for span to depth ratios of more than 50, but the maximum error is only of the order of 4%. This gradually changing error is not characteristic of numerical instability in the stiffness matrix. The error is most likely caused by an inadequate displacement function in an extreme situation. It is assumed that the aspect ratio of the element is not a significant contributory factor causing numerical instability.
F. MISHUand J.
264
MARSHALL
Table 3. Maximum values of NX for stable stiffness matrices
Computer ICL 1900 ICL 470 600
CDC 7600
1
500
t
z
400
P G
I
t
2 g 300 v) / 200 L#llj
/ /
100
/ 0
0
IO
20
30
40
50
60
70
80
Span to depth ratio
Fig. 6. Stresses from variation in aspect ratio of cantilever beam. -.-+, beam theory; 0, one element solution.
CONCLUSIONS
The problems studied have been restricted to a one element solution. This has the advantage that the structural stiffness matrix and the element stiffness matrix are identical. Therefore, any numerical instability in the solution is directly a function of the element stiffness matrix. This is not to be confused with the problem of numerical instability associated with a large structural stiffness matrix assembled from many element stiffness matrices in which the individual stiffness matrices may be numerically stable. No special techniques were used to combat numerical instability apart from using double precision arithmetic. The method of inversion was by Gaussian elimination and all integration was carried out explicitly. Table 3 shows clearly that the stability of the stiffness matix is directly related to the numerical accuracy of the operations involved in the formulation. In particular there are methods of inversion available which will provide greater numerical accuracy than Gaussian elimination for a given set of significant figures. It is not known to what extent the range of stability could be increased by such methods in comparision with the use of increased arithmetical precision. The explicit integration over the area of the element produced no significant instability. However, if a numerical integration is used, then it is
Significant figures
Maximum NX
II 20 1 I5 15 29
I5 5 I2 II >20
1
possible that instability could arise at this stage of the calculation. The method of deriving stiffness matrices for complex elements in plane stress analyses which has been used is particularly suitable for the study of numerical stability in general because of the systematic nature of choosing an appropriate polynomial for the displacement function. From the examples given, it is clear that the numerical stability is a function of the degree of polynomial assumed. Had a study been made of a random selection of complex elements then the results could be confused by the use of special techniques to eliminate or add terms in a chosen polynomial. Although these techniques could influence the numerical stability, the results given do indicate the relationship between the degree of polynomial and the numerical accuracy which defines a numerically stable region. It seems reasonable to assume that these relationships are approximately true regardless of the purpose of a finite element. The study has been conducted using the C-inverse techniques but further investigations would be required to establish the range of stable elements when the unit shape function technique is used[5]. It appears likely that the range of elements produced using this approach will be larger because no C-inverse is required, hence eliminating a source of error. The advantages here may be outweighed by the increased computing time and storage requirements although a similar disadvantage also applies when resort has to be made to double precision arithmetic. REFERENCES
I. 1.Holland, The Finire Element Method in Plane Stress Analysis. Finite Element Methods in Stress Analysis, Tapir, Trondheim, Norway (1%9). 2. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill. London (19711. 3. J. Marshall and F. Mishu, An investigation into the economics of using complex elements in plane stress analysis. Inr. Conj. on Finite Element Methods in Eng. Sci. The University of New South Wales, Australia (1974). 4. F. Mishu, Aspects of complex rectangular elements in plane stress analysis. Ph.D. Thesis, University of Strathclyde, Glasgow (1974). 5. 0. C. Zienkiewicz, B. M. Irons, J. Ergatoudis, S. Ahmad and F. C. Scott, Iso-Parametric and Associated Element Families for Two- and Three-Dimensional Analysis. Finite Element Methods in Stress Analysis, Tapir, Trondheim, Norway (1%9).