- Email: [email protected]
Theoretical and numerical convergence rates for some conforming finite elements
Computers& Sfructures, Vol. 4, pp. I 185-1206. Pergamon Prw IY74. Printed in Great Britain
THEORETICAL AND NUMERICAL CONVERGENCE RATES FOR SOME CONFORMING FINITE ELEMENTS GEORGEE. RAMEY~and NATARAJAN KRISHNAMURTHY: Auburn University, Auburn, Alabama 36830, U.S.A. Abstract-One of the major r~uimments of an approximate solution method is to rapidly converge to the exact solution. This paper is concerned with the rapidity of soiution convergence for many conforming finite element models. Convergence rate are determined by theoretical analyses and checked by numerical investigation. Solution error plots are made for each element and in most cases these plots indicate a numeric& convergencerate whichis in agreementwith that predicted by theory. The investigation~&ders only discretizationerrors; and the solution quantities referred to are for the systeme.nergyquantities which correspond to the eigentiues for eigenvalue problems. The convergence rates determined am applicable to e&value and static problems devoid of stress singularities. NOMENCLATURE
4 B torsion problem bar dimensions flexural rigidity of plate, Eh3/ 12(1- u2) modulus of elasticity shearing rigidity strain energy normality or subsidiary condition in isopcrimetric variational formulation 4.i integer subscripts K constant k buckling coefficient L length of plate n integer number denoting mesh subdivision; derivative order directional cosines mesh pattern designations P,Z P number which designates the rate of convergence Pi surface stresses uniform load R; remainder in finite Taylor series expansion S surface boundary of an elastic body S part of S where surface stresses are prescribed S” part of S where displacements are prescribed S.S. simple-supported boundary u, Y, )l’ exact displacement functions 61*, (b* finite element approximate displacement functions D
E G I J
ff*,
v*,
w*,
t Assistant Professor of Civil Engineering. : Professor of Civil Engineering.
GEORGEE. RAMEY and NATARAJANKRISHNAMURTHY
designation for Taylor series expansion of the quantity two-dimensional coordinates body forces finite element solution quantity designation undetermined parameters material constant,
E (1 +u)(l-20)
; lrequency parameter
rth eigenvalue rth eigenvalue using finite elements Poisson’s ratio comparison functions potential energy mass density per unit plate area real domain; area of domain for two-dimensional first variation second variation change in potential energy.
problems
1. INTRODUCTION IF ONE reviews the major requirements of an approximate to be as follows [l]:
solution method he finds them
(I) Provide convergence to the exact solution (preferably at a rapid rate); (2) Provide a means of improving the approximate solution; (3) Provide an error estimate of the approximate solution. Recent literature [l-7] has given mathematical convergence proofs for conforming finite elements for static and eigenvalue problems and also requirements for monotonic convergence. In these proofs solution convergence is for the idealized problem, i.e. for the mathematical model of the problem. If this model does not accurately represent the real problem, then obviously the finite element approximate solution will not converge to the real solution. Use of conforming elements will accomplish requirements (1 and 2) above. McLay [3] was the first to present a formal analysis of the rate of solution convergence. His work was for static cases and the solutions mentioned were the problem energy quantities (potential and strain energies). Cowper et al. [8], using the method proposed by McLay, conducted analyses for rates of convergence of strain energy quantities for static cases employing two different conforming plate bending elements. They concluded that the convergence rate depends on the particular element being utilized and the stress singularities existing in the exact solution. Additionally, they performed extensive numerical evaluations of the eigenvalue convergence rates for these same two elements. Reference [9] previously reported the results of a theoretical analysis of the convergence rates for eigenvalue solutions employing four different conforming plate bending elements. For this class of problem the eigenvalues correspond to the system strain energy quantities and hence the convergence rates quoted are those for the strain energies and eigenvalues.
Theoretical and Numerical
Convergence.
1187
Rates
This paper is concerned with the rapidity of convergence mentioned in requirement (1). More specifically, It is concerned with theoretically and numerically investigating the solution convergence rate for various conforming elements These rates are of value in assessing the accuracy of solutions and in comparing the merits of different finite element models. Idealization errors are not considered in this paper. In the theoretical treatment of rates of convergence, manipulation errors are nonexistent and hence only discretization errors are considered. When numerical results are later utilized, manipulation errors also come into play: however, no attempt has been made herein to quantify this error source.
2. BRIEF DISCUSSION
OF CONVERGENCE
PROOFS
Since convergence proofs for static and eigenvalue problems are readily available in the literature, detailed proofs will not be given here. However, a brief discussion of these proofs will be given to set the stage for theoretical determinations of convergence rates. Convergence proofs of finite element solutions have been centered around static cases where mathematical proofs have been constructed primarily through the use of variational principles, i.e. the minimum potential energy principle (I-4, 61. In the variational approach all trial functions must satisfy certain smoothness or continuity requirements along with the idealized problem kinematic boundary conditions. ‘These requirements are imposed by the necessity of maintaining proper integrals in the energy functional. Elements satisfying these requirements are identified as conforming elements and for these elements the first variation of the energy functional goes to zero. The continuity conditions normally required are: (1) The displacements and all their derivatives up to order (n- 1) (n is the order of the highest derivative entering the energy functional) must be continuous and finite in the domain. (2) The nth derivative must remain finite in the domain such that the energy functional does not contain improper integrals [3]. When using conforming elements, the question of convergence moves to examining the second variation. Recently, Lynn and Dhillon [5] extended the theoretical convergence proof to plate eigenvalue problems. They accomplished this by transforming plate eigenvalueeigenfunction problems into corresponding isoperimetric variational problems. From this transformation, the rth eigenvalue A,*calculated from a finite element analysis is shown to be the minimum of the plate strain energy with respect to the approximate rth set of eigenfunctions under the normalization and orthogonality conditions. Static problem
Consider the case of an isotropic elastic body of region R with surface boundary S(S=S’+S”), in equilibrium under body forces (Xi, A’,, A’,) and prescribed external surface stresses (p,, p2, pJ on S’ (see Fig. 1). On the surface S”, the displacement boundary conditions are prescribed. Denoting the exact solution as u,, the potential energy of the system is &)=strain
energy+ potential energy of body forces
+ potential energy of surface stresses or, in index notation:
1188
GEORGIZ E. RMEY and NATARAJANKRISHNAMURTHY
FIG. 1. An isotropic elastic body.
where E and u denote moduhrs of elasticity and Poisson’s ratio respectively, and E
G= &, 2(1+ 0)
1=(1+“)(1-2”)’
Let at be the conforming finite element displacement function. where N =number
Recall that:
of elements.
(2)
Denote the difference between u,* and u, by ~.=l&uu,
or
Uof=U,+5a.
(3)
The change in the potential energy of the system due to the increment of r, is: Alc=n(u:)-A@,),
(4)
or after integration by parts and use of Green’s formula, AF==SlC++6% where the first variation, 6n, is
(5)
Theoretical and Numerical Convergence Rates
1189
and the second variation, (5’71,is
with n, denoting the directional cosines of the outer-normal to the boundary surface S. For conforming elements, it can be shown [3-51 that the first variation, &r, is equal to zero. Consequently, we have: Ax=+S’~({.)r0,
(6)
since d2n is positive definite. Equation (6) indicates tbat the approximate potential energy is an upper bound to the exact value. The question of convergence begins with the inequality of equation (6). Convergence is proved if one can show that the second variation, S2x({,), can be made arbitrarily small by reducing the size of the finite elements. Based on equations (5b and 6) we have: (7) By use of the triangular and Schwarz inequalities of norms of functions, we find :
Note, from the above we can see that it is sufficient that the functions r, (i.e. IL.*)possess piecewise continuous first order partial derivatives which converge to zero in the L2 sense, Showing that the upper bound of d2x (the right hand side of the inequality of equation 8) and hence An can be made arbitrarily small by reducing the element size is closely related to the element being utilized. The convergence proof for various elements is shown in the next section along with evaluations of the rate of convergence. Eigenvalue problem Consider the case of the free vibration of a thin elastic plate. This problem may be formulated as an isoperimetric variational problem. Thus, instead of solving the differential equation with appropriate boundary conditions, we may seek the rth eigenfunction, w;, such that it minimizes the plate strain energy I=
[(~,x+w,~~)~ - X1 - u)(w,xxw,,.,,- w,$)] dx dy,
(9)
under the subsidiary condition J= and the (r- 1) orthogonality
conditions
n
w2dxdy=1,
WV
1190
GEORGE E. RAMEY
W,Wj
and
dX dy =O
NATARAJAN KRISHNAMURTHY
forJ=l,
2,. . . ) (r-l),
where w is the lateral deflection of the middle-plane of the plate, ~~~~ represents previous eigenfunctions, Q is the plate area, and D is the plate flexural rigidity. It can be shown that the true eigenvalue is given by the minimum value of the strain energy under the stated restrictions [4. For the finite element approximating system, which replaces the real plate, similar conclusions may be reached, i.e. A: = I(w:).
(12)
Subject to the conditions
J(wF)= 1,
(13)
and w:w; dx dy=O
i=l,
2,. . . ) (I’-_),
n
(14)
where w: is the rth minimizing finite element approximate eigenfunction defining the rth approximate eigenvalue A,*. it can also be shown that $ provides an upper bound to 1, [S], i.e. (15) or A,*- 1, = f(wF)- I( w,) 2 0.
(16)
For convergence of J: to & it is only necessary to estimate the upper bound of the left hand side of the inequality of equation (16) in terms of element size or number of elements and then show this bound can be made as small as desirable. This procedure is demonstrated in the next section where convergence is shown and convergence rate determined for various elements.
3. THEORETICAL
CONVERGENCE
RATES
In this section we will theoretically investigate the convergence (second variation going to zero) and convergence rate of various conforming elements. The approach taken here parallels that given in [9] for convergence and convergence rate of eigenvalues; however, the application has been expanded to include more elements. Application to static type problems (for problem energy quantities) has also been made and shows results identical to those for eigenvalue problems regarding element convergence and convergence rate. Eigenvalue problem
The solution convergence (and rate) for this class of problem will be demonstrated here by means of a conforming 16degrees-of-freedom rectangular thin plate bending element [IO, I I, 191. For this particular element four nodal displacements (w; H’,.~:w.,; II’,,).)are
Theoretical
and Numerical Convergcncc Rates
assigned at each corner node, and the ith element trial displacement taken as O~‘*)i=a1 +a2x+a3y+a,s2+asxy+a,y2
1191
function (IV*), is
+a,x3+a,x2y+agy’+a,,y3+a,,~3y
+~12.vy3+a,3xzy2+a14x3y2+aI~x2y3+a16x3y3.
(17)
It should be noted that (\v*)~ contains a complete cubic polynomial in x and y. Since (,~,*)i is defined only in the ith element and zero elsewhere, the finite element approximating function for the entire plate domain is defined as
where N is the total number of elements. We now introduce a comparison function in the form of a third degree (since the trial function is a complete cubic) finite Taylor series expansion (in x and y) of the exact solution IV,about a point 0 within the ith element, i.e. ~w~)i~w~l~+(x~xxO)W~,~~~+(Y~Y~)W~.~(~+
- . . +
:,(Y-_Y,) 3W,,wy(o. .
. (19)
Let (R); be the remainder term of the expansion, i.e.
(RI,=
ai[(x-x,)‘W~.~~~=I~+ ***+(Y-Yo)*W, yyyylb’J*
(20)
Note, the remainder is of fourth order in x and y, and we are assuming in equations (I9 and 20) that up through the fourth order partial derivatives of the exact solution w, exist and are bounded (no stress singularities). The notations lo and 16 designate that the quantities to the left of the signs are evaluated at points 0 and 0 respectively. The point 0 is a point on the line joining points (xc, ye) and (x, y). Accordingly, the third degree Taylor series comparison function for the entire plate domain is
(21) where (Q,)i is defined in the ith element and is zero elsewhere. Next we note that among all the third degree polynomials in x and y, e.g. w: and Gr, w: minimizes the integral I under the subsidiary conditions. Hence in view of equation ( 16) we have
~t%kf(w:)i2f(w,) Or
(22)
GEORGE E.
1192
RAMEY
and NATARAJAN KRISHNAMURTHY
When 9, in equation (23) is replaced by 8,= w,-
(24)
the left hand side of the equation becomes [4, 51
(25)
where for simplicity the subscript i, which denotes the quantities belonging to the ith element, on the remainder R has been taken out of the parentheses, and Ri denotes the area of the ith element. Equation (25) can be expanded and put in the form
I(%)-Z(w,) =
fl{;j-jn@ - Wfx, +R&J+ @Lx + &J2 + 2(1- u)R,$] dx dy
I
20.
Consider now the term u(R,, + R,J’, under the integral sign in equation (26). triangular and Schwarz inequalities or norms of functions, we find
OS
ss R,
(26)
By the
@,,+K,,J2dx dy~C(IR,,l(n,+IIR,~vlln,l2=IIR,:xlI~,+2lIR,xxlln,(IR~vrlIn~ + 1IR,~#L.
(27)
Upon substituting equation (27) into (26) and simplifying, we find
Recall that the remainder (R)i is of fourth order in x and y, i.e. Ri=f(E4, . . . , J4),
(29)
The maximum value of Z and 7 is of the order L/n (length where 3=x-x, and J=y-y,. of domain side divided by number of elements along the side). The derivatives of Ri are R i,,=m2,
* * - 3 J2)
RSyy =h(9,
. . . , Y2)
Ri,xy =_&(I’, . . . , 7’).
(30)
Theoretical and Numerical ConvergenceRates
1193
Hence,
2=j (Rx,)’ dQ=&d~,. . . , J4)dQ, IIn IIR‘XY i-l s
(31)
where the subscript i on R and R has been dropped for convenience. We can see that the above intergrands approach zero as n+co at a rate of (l;/n)4. Hence from equation (28) we have (32) and on redefining the K, we have
01(A~),-j2,1zw4,
(33)
where n is a measure of the number of elements along one side (or the mesh size). Accordingly, the convergence rate of the eigenvalue obtained from this element is nm4. This is in agreement with the value determined numerically by Lindberg and Olson [l I]. Static problem For completeness an analysis similar to that above will now be performed for the potential energy convergence rate in a static type problem. For convenience, we choose one of the simplest elements, a 6-degrees-of-freedom triangular in-plane element with trial displacement functions (U)i’c(r + a2x + a3y (34) To show that the second variation goes to zero and hence convergence with reduction of element size we proceed as before. That is, we first introduce comparison function in the form of linear finite Taylor series expansions (in s and y) of the exact solutions u and v about a point 0 within the ith element, i.e. (“>i=u10
+
(x-~xO)u~xJO +
(Y-YOb4ylO
(v)I=v(o+(x--so)v~,~o+(Y-Yo)v~~~o.
(35)
Let (R,)i and (R,), be the second order remainder terms of the expansion, i.e. (R,)i=3[(-u-XO)2U,,,IB+
* * * +
(Y -YO)2u9yylbl
(R”)i=~[(x-x,)21’,,,l + . * - +(Y-Yo)~v,~~ITI*
(36)
1194
GEORGE E. RAMEY and NATARAJAP* KRISHNAMURTHY
As before, in equations (35 and 36) we assume that up through the second order partial derivatives of the exact solution exist and are bounded (hence no stress singularities). The linear Taylor series comparison functions for the entire domain are
Among all the complete linear functions in .Y and y, U* and Y* minimize the potential energy of the problem. Hence we conclude that
or, in other words,
where
and
We now use the upper bound to S*rr({,, 5,) to show the convergence and convergence rate for the element. Based on equation (8) in Section 2, we have
From the definitions of u’ and < we see that
(41)
Hence, equation (40) may be rewritten as
OS i l-1
4CllR..,II’+2ll%II Ik.Il+
{ ~[ll&,,J)*+llRo,,ll’J+
IR...i12-$ (42)
Recall that the remainders (R,,), and (R,), are of the second order; hence, the first derivatives of these expressions are linear in x and y. Therefore,
Theoretical and Numerical Convergence Rates
1195
(43)
Hence, from equation (40) we have Osn(ti,
i;)-I[(u,
v)=~~J&, &
KII-~
w
or Oln(u*,
v*)-n(u,
v)~Kn-~.
(45)
Accordingly, the asymptotic convergence rate of the potential energy obtained from this element is Knv2. If one is confident of convergence when employing a particular element and is interested in the rate of mean convergence, the procedure for obtaining this rate may be greatly simplified. For example, in the above case, (g-u) and (g-v) are of order (n-‘); and the first derivatives of these quantities are of order (n-r). Since the second variation, and hence Aii, involves the square of these first derivatives, then this quantity is also of order (nW2). The area of the domain remains fixed as the mesh is refined; and ,therefore, the summation over all the elements leads to the conclusion that Aii (which bounds Arc+) is of order (PI-‘). Similarly for the thin plate rectangular bending element discussed previously, (r? - w) was of order (ns4); and, therefore, the second derivatives of displacements were of order (nW2). Since the strain energy expression involves the square of the second derivatives, then this quantity and hence the approximate eigenvalue is of order @I-~). In general, the asymptotic convergence rates of the eigenvalues and energy quantities may be regarded as Kemp. It can be noted that the partial derivatives of the remainder terms, (Z& are at least of order (n-l) for convergence and that the square of the norms of these derivatives will be at least of order (nm2). Hence, in the general formulation of the convergence rate, i.e. Knmp,p will be greater than or equal to two. Exceptions to this arise when there are unbounded higher order derivatives of the exact solutions which are exhibited by stress singularities. From the theoretical convergence studies, we see that the rate of solution convergence depends on the particular finite element used, the character of any singularities possessed by the exact solution; and it is independent of the mesh pattern This is assuming the problem’s essential boundary conditions are satisfied, and as previously, mentioned, neglects idealization errors. Theoretical convergence rate analyses similar to the two demonstrated here were performed for various other conforming elements. The results of these analyses are summarized in Table 1. The trial displacement functions for the elements considered are given in the Appendix.
GEORGE E. RAMEY
1196
and NATARAJAN KRJSHNAMURTHY
TABLE1. Conforming displacement finite element energy and eigenvalue theoretical convergence rates
Symbol
Type element
Element shape
T6 R8 RP16 CT20 R16 Tl8 W? LDT TT3 lT6
In-plane In-plane In-plane In-plane Thin plate Thin plate Thin plate Thick plate Torsion Torsion
Triangle Rectangle Rocttigle Triangle RectangIe Triangle Quadrilateral Triangle Triangle Triangle
Element degrees-offreedom 6 8 Y 16 :i 9 63
Theoretical convergence rate Kn-2 G-2 G-4 Kfl-6
iyn-4 Kn-6 Kn-2 Kn-2 Kn-2 Kn-4
References D4, I61 WI f:;
IS1
[IO, 11, 191 (81 WI [S, 131 W, 181 U6, 181
4. NUMERICAL INWSTIGATION Numerical results for the elements listed in Table 1 were studied with the emphasis being on the rate of solution convergence. For the most part, published numerical data were utilized and log-log plots of solution error or relative error versus mesh size were made for the various elements. For those cases where exact solutions were available the error terms shown plotted are absolute errors. Relative errors are plotted for the remainder of the cases by estimating (from results available) exact solutions and using these estimates to deermine the relative error. The format of these plots is identical to that presented by Cowper et al. [8]. It provides an excellent graphical means of checking tbe numerical verification of the theoretical convergence rate. In general, the error plots tend to be asymptotic as n increases. The siopc of the asymptotic line is the convergence rate which is easily compared to the slope of a straight line representing the theoreti~lly determined rate. For convenience of discussion the elements have been categorized as: (1) in-plane elements, (2) plate bending elements, (3) torsion elements. A separate plot for each element is presented in order that an unencumbered display of the numerical convergence rate for each element can be easily compared with the predicted theoretical rate. In all figures presented herein, the theoretical convergence rates are indicated by the dashed lines. Lighter dash-dot lines are provided in some instances for convenience in assessing the numerical solution convergence rate. In-plane elements
Error or relative error plots for in-plane elements are shown in Figs. 2-5. The elements considered are: 26, R8, RPl6 and 0’20. Several of these elements are applied to a variety of problems. For the cases considered, Figs. 2 and 3 offer good numerical evidence of verifying the theoretical convergence rates of fi’” and Kn-4 for the T6 and R8 elements respectively. Figure 4 indicates, for the small n values utilized, a numerical convergence rate for element RP16 of approximately 3.0 rather than the theoretically predicted rate of 4.0. Similarly, Fig. 5 indicates a numerical rate for element CT20 of 5.0 as compared to 6.0 predicted by theory. For higher n values, elements RR16 and CT20 may show results which are in closer agreement with theory.
Theoretical and Numerical Convergence Rates
r(Blaua
U!DJ+S
IJ!
anl+ojau
JOJJS
R
Tx
c
\
‘,
;;
tt6JaUe
U!OJ+S~
U!
JOJJ8
eA!+Dlau
c
1198
GEORGEE. RAMBYand NATARAJANKRISHNAMURTHY
j “:
0
,
,
n
0
”
0
'0 -
b d6J8US
U!DJ+S
U!
1:
,
JOJJS
SA!+DlS#
a
'0 -
1199
Theoretical and Numerical Convergence Rates
Comparing Figs. 2-5, one can observe the superior convergence rate of the higherorder elements. However, in the case of the RP16 element, it is doubtful that its slightly improved convergence rate relative to the RS element will justify its additional degrees-offreedom. Hence, for a given n size the R8 element will probably be more efficient than the higher-order RP 16 element. Plate bending elements
Error or relative error plots for plate bending elements are shown in Figs. 6-9. The elements considered are: R16, T18, CFQ and UK Both static and eigenvalue problems are considered for a variety of load and boundary conditions. For small n values, the numerical curves vary depending on the particular problem. However, most of the error plots indicate a numerical convergence rate that is in agreement with the theoretical rate for that element.
IO
I
k
10.
L_
1
I
I d
7
fvleosure
of number
, c
skm
of elements
FIG. 6. Relative error of finite element solutions-thin
per
side,
n
plate bending W6 elements.
1200
GEORGE E. RAMEY and NATARAJANKRISHNAMLJRTHY
The buckling coefficient for the square plate with a free edge shown in Fig. 6 exhibits numerous fluctuations such that a trend in the numerical results cannot be ascertained. The fluctuation of the buckling coefficient for the clamped square plate shown on that same figure is so great that non-monotonic convergence is indicated. This is understandable since it occurs at an n value of three which is not in the n=2, 4, 8, 16 etc. sequence. For these two buckling coefficient curves, it is believed that with increases in n, the numerical rates will be in agreement with theory as are the other curves on Fig. 6.
SYm I-Clomped
square
(n =4
shown)
\ \\
loadlpattern
Measure
Fro.
of
number
of
elements
7. Relative error of iinite element solutions-thin
per
side,
n
plate bending TlS elements.
Numerical results taken from [8] and shown plotted in Fig. 7 indicate the excellent convergence characteristics of the T18 element and offer good numerical confirmation of the theoretically predicted convergence rate of fi -6. Cowper et al. [8] concluded from their numerical study that solution accuracy depends in a large measure on the element arrangement (i.e. mesh pattern). Their numerical data indeed indicate this; however, if one makes an error plot utilizing their data, it can be seen that the rate of solution convergence for the problems considered does not depend on the mesh pattern. This can be seen in Fig. 7 where results for the two different mesh patterns P and Q are shown.
Theoretical and Numerical Convergence Rates
1201
. l
A,,-S. S. \ k-S.S.
S
\
plate
.
rectangular plate under pure shear
. \
\
\ \ k-Clamped square plate under biaxial compressIan \ \ \
.
.
.
\
.
\ l
.
\ k-S.S. square plate under A unlarial compression
\
0
\ \
.
\ \
\
’ \ “-2\ \ \
A- Frequency
\ \ ’
parameter coefficient
Buckling
‘\.
. \
k-
square
\
“\. \ \ \
1
I 2 Measure
of
4 number
of
I
i/I
6
8
elements
FIG. 8. Relative error of finite element solutions-thin
II IO 12 per
16 side,
20 n
plate bending CP’Q elements.
Figure 8 shows some eigenvalue errors resulting from use of CFQ elements. For the cases shown it appears that the order of convergence for the curves is Kn” which is in agreement with the theoretical rate The curve for 3,, 1 is converging at a faster rate; however, it shows fluctuations and the values of n for the curve are small. At higher n values, this curve will also probably converge at a rate of Kn- 2. Two of the buckling coefficient curves in Fig. 8 are for problems shown earlier in Fig. 6. Comparing these corresponding curves from Figs. 6 and 8, we can observe the superior convergence rate of the R16 element relative to that of the CFQ elements. The conforming thick plate element LDT exhibits a numerical rate in Fig. 9 which is m agreement with the Knq2 predicted by theory.
Torsion elements Torsional rigidity error plots for torsion elements TT3 and TT6 are shown in Figs. 10 and 11 respectively. In each case, the numerical convergence curve approaches the theoretical rate quite closely. It should be noted that it would be highly desirable to have
1202
GEORCXE.RAMEY and NATARAJANKRISHNAMURTH~
----
3
2
Measure
of
number
s s
square
Clamped plate
d of
5 elements
FIG. 9. Relative error of finite element solutions-thick
Plate
s~uore
6
7 per
8 side,
n
plate bending L DT elements.
at least one more n value data point for the 77’6 element to verify that element’s numerical convergence rate. The superior convergence rate of the higher-order TT6 element is evident from the figures; however, the additional degrees-of-freedom resulting from its use may result in its being less efficient than the lower-order TT3 element.
5. CONCLUSIONS The use of finite Taylor series expansions along with compatible form comparison functions is an effective means of determining conforming finite element theoretical convergence rates. These rates depend only on the particular finite element and are the same for both static and eigenvalue problems and are independent of the mesh pattern. This is assuming the idealized problem’s essential boundary conditions are satisfied. The theoretical convergence rate determination is seen to be invalid in the presence of stress singularities.
Theoretical and Numerical Convergence Rates
n
N
b
b Je+euoJod
r(+!p!b!J
IDUO!SJO+
U!
JOJJ3
1203
1204
Gsoxo~ E. F~AMEYand NATARAJANKRISHNAMUR~
In general, the asymptotic convergence rates may be regarded as Knsp, where p will be greater than or equal to two. From Table 1, one may conclude that in general as the degrees-of-freedom increase for a particular type (i.e. bending, torsion, in-plane) element the convergence rate also increases. The element convergence rates are of value in accessing solution accuracy and in comparing the relative merits of various finite element models. TABLE2. Summary table of element convergence rates Element
Theoretically predicted rate
T6
Kn -2 Kn-2
RP16 CT20 R16
Kn-4
R8
T18
CFQ
LDT lT3 TT6
Kn-6 Kn-4
Kn-6
Kn-2
Kn-2 Kn-2 Kn-4
Numerically determined rate Kn-2
Kn-2 Kn-3 K?l-5 Kn-4 Kn-6 Kn-2 Kn-2 Kn-2 Kn-4
Most of the error plots indicate a numerical convergence rate that is in agreement with the theoretical rate for that element. Relatively larger error fluctuation was exhibited for small values of n when comparing different elements, or when comparing the results for different problems while utilizing the same element. However, as n increased, the fluctuations decreased and the numerical curves appear to become asymptotic to the theoretical rate. A few curves exhibited so much fluctuation that a trend in numerical results cannot be ascertained. It is felt that for larger values of n, these curves will also become asymptotic at a rate in agreement with theory. A summary of the element convergence rates, theoretically predicted and numerically determined, is shown in Table 2.
REFERENCES [l] S. W. REY, A convergence investigationof the direct stiffness method. Ph.D. Thesis. University of Washington, Seattle, Washington (1966). [2] M. W. JOHNSON,JR. and R. W. MCLAY,Convergence of the finite element method in the theory of elasticity. J. Appl. Mech. ASME, 274-278, June (1968). [3] R. W. MCLAY.Completeness and convergence properties of finite element displacement functionsa general treatment. AIAA Paper No. 67-143, AIAA 5th Aerospace Sciences Meeting, New York, [4] ?P?z% A. P. BORESI,Kinds of convergence and improved convergence of conforming finite element solutions in plate bending. Nucl. Em Design 11,159-176 (1970). [5] P. P. LYNNand B. S. DHILLON,Convergence of eigenvalue solutions in conforming plate bending finite elements. Int. J. Num. Meth. Engng 4 (2). 217-234 (1972). [6] P. TONOand T. H. H. PUN, The convergence of finite element method in solving linear elastic problems. Znt. J. SW& Struct. 3 (5), 865-879 (1967). [7j J. E. WALZ,R. E. FULTON and N. J. Cvuus, Accuracy and convergence of finite element approximations. Proc. 2nd Con/. on Matrix Methoak in Structural Mechanics. Wright-Patterson Air Force Base, Ohio, pp. 995-1025, October (1968). [8] G. R. Cowpe~, E. Kosuo, G. M. LINDBERGand M. D. OLSON,A high-precision triangular platebending element. National Research Council of Canada Aeronautical Report LR-514, Ottawa, December (1968).
Theoretical and Numerical Convergence Rates
1205
(91 P. P. LYNNand G. E. RAMEY,A method to generate both upper and lower bounds to plate eigenvalues by conforming displacement finite elements. Proc. 1st Intertwtionai Conference on Structural Mechanics in Reactor Technology (1971).
[lo] F. K. BOONER,R. L. Fox and L. A. SCHMIT,JR., The generation of the inter-element compatible stiffness and mass matrices by the use of interpolation formulas. Conf on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 39743, October (1965). [ll] G. M. LINDBERG,M. D. OLSONand H. A. TULL~CK,Closed form, finite element solutions for plate vibrations. National Research Council of Canada Aeronautical Report LR-518, Ottawa, February (1969). [12] R. W. CL~UGHand C. A. FELIPPA,A refined quadrilateral element for analysis of plate bending. Proc. 2nd Conf. on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 399-440, October (1968). [13] B. S. DHILLON,Triangular finite elements for the bending analysis of thick elastic plates. Ph.D. Thesis, Dept. of Civil Engineering, University of Colorado, Boulder, Colorado (1970). [14] R. W. CLOUGH,The finite element method in structural mechanics. In Stress AnaIysis (Edited by 0. C. ZIENKIEWICZ and G. S. HOLISTER), Chapter 7. Wiley, New York (1965). [15] R. W. CL~UGHand J. L. TOCHER,Finite element stiffness matrices for analysis of plate bending. Conf. on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 515-545, October (1965). [16] 0. C. ZreMuEwrcz, The Finite Element Method in Structral and Continuum Mechanics. McGraw-Hill, New York (1967). [17] G. R. COWPER,G. M. LINDB~RG and M. D. OWN, A shallow shell finite element of triangular shape. ht. J. Solids Struct. 6, 1133-l 156 (1970). [18] G. E. RAMEY, A study of the convergence characteristics and accuracy of conforming finite element solutions. Ph.D. Thesis, Dept. of Civil and Environmental Engineering, University of Colorado, Boulder, Colorado (1972). [19] W. G. CARSONand R. E. NEWTON,Plate buckling analysis using a fully compatible finite element. AZAA J. 7 (3), 527-529 (1969). Received December 1973
APPENDIX Following is a listing of the code names and trial displacement functions for the conforming finite element models considered in this article. R16: 16-degrees-of-freedom
rectangular thin plate element
A plate bending element reported in [8, 11, 191, but generally credited to [lo]. (w*)i=cCl +a~x+a~~+a,.~2+a,x~+a,y2+a,x3+a,x2y+agxy2+aloy3+~11x3~ +a~2xy3+a13x2y2+a,,x3y2+a15x2y3+aIgx3y3. CFQ:
Clough and Felippa’s thin plate quadrilateral element
This is a plate bending element presented in [12]. Clough and Felippa first proposed the LCCT-12 triangular element with 12 degrees-of-freedom, and eaploying a compatible cubic displacement field. They further improved it to give the quadrilateral element designated Q-19 (referred to as the CFQ element in this paper), retaining the 12-degreesof-freedom for the element, the explicit displacement functions for which are very complicated and are not given here. T18 : Cowper et al. thin plate triangular element
A plate bending element reported in [S]. (~~*)i=a,+~2x+a3~+a,.~2+x5x~‘+a,~’+~,x3+a,x2y+agx~2+alOY3+a11x4 +qIlx 3~~+~sc,,x-2~2+a,4.~~3+cc,,~~4+a,,ssfa,,.u3y2+~ 1s.~2y3+y:19x~4+a20~s.
1206
GF_ORGE E. RAMEYand NATARAJAN KRKSHNAMLIRTHY
LDT:
Lynn and Dhillon’s thick plate triangular element
A thick plate bending element presented in [5 and 131.
T6 : 6-degrees-of-freedom
(W*)i=c(l
~azX+a,y+Cr,X2+~5~y+S.CLsS..’
(8*)i=c(7
+ ct8x +
agy
triangular in-plane element
The classical two-dimensional constant strain element given in [ 161 (U*)i=cIt f (v*)i’c14.
RS : g-degrees~f-freedom
CQX+ a,y
f c(gX + a6ym
rectangular in-plane element
The classical two-dimensional
plane element given in 1161.
(U*)i=tLi +a2x+a3y+a4xy (y*)i=a,+aex+a7y+awxy. RR16 : M-degrees-of-freedom
rectangular in-plane element
A higher-order two-dimensional plane element given in [16]. (u*) i=a~+a2x+a3y+a4x2+asxy+afjy2fa7x2ytagxy2 (v*)~=ag+aloxfa,,y+at2x2+a,3xy~a,,y2+a1S.~2y+a16_~y2. TT3 : 3-degrees-of-freedom
triangular torsion element
A first order (for warping function w) torsion element reported on in [16]. (w*>i=al +a,x+a,y. TT6: 6-degrees-of-freedom
triangular torsion element
A higher-order (for warping function w) torsion element given in [16]. (w*)j’ai CT20:
fazx+xSy+&$.~2
+a5xy+fC(6y2.
modified Cowper et aI. shallow shell element of triangular shape
In its ~rn~d~ed form, the element is a curved triangular element designed for shelI analyses and is reported [I?]. CT20 is the limiting case of a fiat plate, or the plane stress element contained within the total element. (u*),=a,+azx+a,y+a,x2+a,.uyfa,y2+a,x3t-a,x2y+agxy2+a,,y3 (v*)j=a,,+al2x+a,,yfa,,x2+a,,xy+a,~y2+a,,x3+ff,,x2yfa,~xy2+a2,y3.
THEORETICAL AND NUMERICAL CONVERGENCE RATES FOR SOME CONFORMING FINITE ELEMENTS GEORGEE. RAMEY~and NATARAJAN KRISHNAMURTHY: Auburn University, Auburn, Alabama 36830, U.S.A. Abstract-One of the major r~uimments of an approximate solution method is to rapidly converge to the exact solution. This paper is concerned with the rapidity of soiution convergence for many conforming finite element models. Convergence rate are determined by theoretical analyses and checked by numerical investigation. Solution error plots are made for each element and in most cases these plots indicate a numeric& convergencerate whichis in agreementwith that predicted by theory. The investigation~&ders only discretizationerrors; and the solution quantities referred to are for the systeme.nergyquantities which correspond to the eigentiues for eigenvalue problems. The convergence rates determined am applicable to e&value and static problems devoid of stress singularities. NOMENCLATURE
4 B torsion problem bar dimensions flexural rigidity of plate, Eh3/ 12(1- u2) modulus of elasticity shearing rigidity strain energy normality or subsidiary condition in isopcrimetric variational formulation 4.i integer subscripts K constant k buckling coefficient L length of plate n integer number denoting mesh subdivision; derivative order directional cosines mesh pattern designations P,Z P number which designates the rate of convergence Pi surface stresses uniform load R; remainder in finite Taylor series expansion S surface boundary of an elastic body S part of S where surface stresses are prescribed S” part of S where displacements are prescribed S.S. simple-supported boundary u, Y, )l’ exact displacement functions 61*, (b* finite element approximate displacement functions D
E G I J
ff*,
v*,
w*,
t Assistant Professor of Civil Engineering. : Professor of Civil Engineering.
GEORGEE. RAMEY and NATARAJANKRISHNAMURTHY
designation for Taylor series expansion of the quantity two-dimensional coordinates body forces finite element solution quantity designation undetermined parameters material constant,
E (1 +u)(l-20)
; lrequency parameter
rth eigenvalue rth eigenvalue using finite elements Poisson’s ratio comparison functions potential energy mass density per unit plate area real domain; area of domain for two-dimensional first variation second variation change in potential energy.
problems
1. INTRODUCTION IF ONE reviews the major requirements of an approximate to be as follows [l]:
solution method he finds them
(I) Provide convergence to the exact solution (preferably at a rapid rate); (2) Provide a means of improving the approximate solution; (3) Provide an error estimate of the approximate solution. Recent literature [l-7] has given mathematical convergence proofs for conforming finite elements for static and eigenvalue problems and also requirements for monotonic convergence. In these proofs solution convergence is for the idealized problem, i.e. for the mathematical model of the problem. If this model does not accurately represent the real problem, then obviously the finite element approximate solution will not converge to the real solution. Use of conforming elements will accomplish requirements (1 and 2) above. McLay [3] was the first to present a formal analysis of the rate of solution convergence. His work was for static cases and the solutions mentioned were the problem energy quantities (potential and strain energies). Cowper et al. [8], using the method proposed by McLay, conducted analyses for rates of convergence of strain energy quantities for static cases employing two different conforming plate bending elements. They concluded that the convergence rate depends on the particular element being utilized and the stress singularities existing in the exact solution. Additionally, they performed extensive numerical evaluations of the eigenvalue convergence rates for these same two elements. Reference [9] previously reported the results of a theoretical analysis of the convergence rates for eigenvalue solutions employing four different conforming plate bending elements. For this class of problem the eigenvalues correspond to the system strain energy quantities and hence the convergence rates quoted are those for the strain energies and eigenvalues.
Theoretical and Numerical
Convergence.
1187
Rates
This paper is concerned with the rapidity of convergence mentioned in requirement (1). More specifically, It is concerned with theoretically and numerically investigating the solution convergence rate for various conforming elements These rates are of value in assessing the accuracy of solutions and in comparing the merits of different finite element models. Idealization errors are not considered in this paper. In the theoretical treatment of rates of convergence, manipulation errors are nonexistent and hence only discretization errors are considered. When numerical results are later utilized, manipulation errors also come into play: however, no attempt has been made herein to quantify this error source.
2. BRIEF DISCUSSION
OF CONVERGENCE
PROOFS
Since convergence proofs for static and eigenvalue problems are readily available in the literature, detailed proofs will not be given here. However, a brief discussion of these proofs will be given to set the stage for theoretical determinations of convergence rates. Convergence proofs of finite element solutions have been centered around static cases where mathematical proofs have been constructed primarily through the use of variational principles, i.e. the minimum potential energy principle (I-4, 61. In the variational approach all trial functions must satisfy certain smoothness or continuity requirements along with the idealized problem kinematic boundary conditions. ‘These requirements are imposed by the necessity of maintaining proper integrals in the energy functional. Elements satisfying these requirements are identified as conforming elements and for these elements the first variation of the energy functional goes to zero. The continuity conditions normally required are: (1) The displacements and all their derivatives up to order (n- 1) (n is the order of the highest derivative entering the energy functional) must be continuous and finite in the domain. (2) The nth derivative must remain finite in the domain such that the energy functional does not contain improper integrals [3]. When using conforming elements, the question of convergence moves to examining the second variation. Recently, Lynn and Dhillon [5] extended the theoretical convergence proof to plate eigenvalue problems. They accomplished this by transforming plate eigenvalueeigenfunction problems into corresponding isoperimetric variational problems. From this transformation, the rth eigenvalue A,*calculated from a finite element analysis is shown to be the minimum of the plate strain energy with respect to the approximate rth set of eigenfunctions under the normalization and orthogonality conditions. Static problem
Consider the case of an isotropic elastic body of region R with surface boundary S(S=S’+S”), in equilibrium under body forces (Xi, A’,, A’,) and prescribed external surface stresses (p,, p2, pJ on S’ (see Fig. 1). On the surface S”, the displacement boundary conditions are prescribed. Denoting the exact solution as u,, the potential energy of the system is &)=strain
energy+ potential energy of body forces
+ potential energy of surface stresses or, in index notation:
1188
GEORGIZ E. RMEY and NATARAJANKRISHNAMURTHY
FIG. 1. An isotropic elastic body.
where E and u denote moduhrs of elasticity and Poisson’s ratio respectively, and E
G= &, 2(1+ 0)
1=(1+“)(1-2”)’
Let at be the conforming finite element displacement function. where N =number
Recall that:
of elements.
(2)
Denote the difference between u,* and u, by ~.=l&uu,
or
Uof=U,+5a.
(3)
The change in the potential energy of the system due to the increment of r, is: Alc=n(u:)-A@,),
(4)
or after integration by parts and use of Green’s formula, AF==SlC++6% where the first variation, 6n, is
(5)
Theoretical and Numerical Convergence Rates
1189
and the second variation, (5’71,is
with n, denoting the directional cosines of the outer-normal to the boundary surface S. For conforming elements, it can be shown [3-51 that the first variation, &r, is equal to zero. Consequently, we have: Ax=+S’~({.)r0,
(6)
since d2n is positive definite. Equation (6) indicates tbat the approximate potential energy is an upper bound to the exact value. The question of convergence begins with the inequality of equation (6). Convergence is proved if one can show that the second variation, S2x({,), can be made arbitrarily small by reducing the size of the finite elements. Based on equations (5b and 6) we have: (7) By use of the triangular and Schwarz inequalities of norms of functions, we find :
Note, from the above we can see that it is sufficient that the functions r, (i.e. IL.*)possess piecewise continuous first order partial derivatives which converge to zero in the L2 sense, Showing that the upper bound of d2x (the right hand side of the inequality of equation 8) and hence An can be made arbitrarily small by reducing the element size is closely related to the element being utilized. The convergence proof for various elements is shown in the next section along with evaluations of the rate of convergence. Eigenvalue problem Consider the case of the free vibration of a thin elastic plate. This problem may be formulated as an isoperimetric variational problem. Thus, instead of solving the differential equation with appropriate boundary conditions, we may seek the rth eigenfunction, w;, such that it minimizes the plate strain energy I=
[(~,x+w,~~)~ - X1 - u)(w,xxw,,.,,- w,$)] dx dy,
(9)
under the subsidiary condition J= and the (r- 1) orthogonality
conditions
n
w2dxdy=1,
WV
1190
GEORGE E. RAMEY
W,Wj
and
dX dy =O
NATARAJAN KRISHNAMURTHY
forJ=l,
2,. . . ) (r-l),
where w is the lateral deflection of the middle-plane of the plate, ~~~~ represents previous eigenfunctions, Q is the plate area, and D is the plate flexural rigidity. It can be shown that the true eigenvalue is given by the minimum value of the strain energy under the stated restrictions [4. For the finite element approximating system, which replaces the real plate, similar conclusions may be reached, i.e. A: = I(w:).
(12)
Subject to the conditions
J(wF)= 1,
(13)
and w:w; dx dy=O
i=l,
2,. . . ) (I’-_),
n
(14)
where w: is the rth minimizing finite element approximate eigenfunction defining the rth approximate eigenvalue A,*. it can also be shown that $ provides an upper bound to 1, [S], i.e. (15) or A,*- 1, = f(wF)- I( w,) 2 0.
(16)
For convergence of J: to & it is only necessary to estimate the upper bound of the left hand side of the inequality of equation (16) in terms of element size or number of elements and then show this bound can be made as small as desirable. This procedure is demonstrated in the next section where convergence is shown and convergence rate determined for various elements.
3. THEORETICAL
CONVERGENCE
RATES
In this section we will theoretically investigate the convergence (second variation going to zero) and convergence rate of various conforming elements. The approach taken here parallels that given in [9] for convergence and convergence rate of eigenvalues; however, the application has been expanded to include more elements. Application to static type problems (for problem energy quantities) has also been made and shows results identical to those for eigenvalue problems regarding element convergence and convergence rate. Eigenvalue problem
The solution convergence (and rate) for this class of problem will be demonstrated here by means of a conforming 16degrees-of-freedom rectangular thin plate bending element [IO, I I, 191. For this particular element four nodal displacements (w; H’,.~:w.,; II’,,).)are
Theoretical
and Numerical Convergcncc Rates
assigned at each corner node, and the ith element trial displacement taken as O~‘*)i=a1 +a2x+a3y+a,s2+asxy+a,y2
1191
function (IV*), is
+a,x3+a,x2y+agy’+a,,y3+a,,~3y
+~12.vy3+a,3xzy2+a14x3y2+aI~x2y3+a16x3y3.
(17)
It should be noted that (\v*)~ contains a complete cubic polynomial in x and y. Since (,~,*)i is defined only in the ith element and zero elsewhere, the finite element approximating function for the entire plate domain is defined as
where N is the total number of elements. We now introduce a comparison function in the form of a third degree (since the trial function is a complete cubic) finite Taylor series expansion (in x and y) of the exact solution IV,about a point 0 within the ith element, i.e. ~w~)i~w~l~+(x~xxO)W~,~~~+(Y~Y~)W~.~(~+
- . . +
:,(Y-_Y,) 3W,,wy(o. .
. (19)
Let (R); be the remainder term of the expansion, i.e.
(RI,=
ai[(x-x,)‘W~.~~~=I~+ ***+(Y-Yo)*W, yyyylb’J*
(20)
Note, the remainder is of fourth order in x and y, and we are assuming in equations (I9 and 20) that up through the fourth order partial derivatives of the exact solution w, exist and are bounded (no stress singularities). The notations lo and 16 designate that the quantities to the left of the signs are evaluated at points 0 and 0 respectively. The point 0 is a point on the line joining points (xc, ye) and (x, y). Accordingly, the third degree Taylor series comparison function for the entire plate domain is
(21) where (Q,)i is defined in the ith element and is zero elsewhere. Next we note that among all the third degree polynomials in x and y, e.g. w: and Gr, w: minimizes the integral I under the subsidiary conditions. Hence in view of equation ( 16) we have
~t%kf(w:)i2f(w,) Or
(22)
GEORGE E.
1192
RAMEY
and NATARAJAN KRISHNAMURTHY
When 9, in equation (23) is replaced by 8,= w,-
(24)
the left hand side of the equation becomes [4, 51
(25)
where for simplicity the subscript i, which denotes the quantities belonging to the ith element, on the remainder R has been taken out of the parentheses, and Ri denotes the area of the ith element. Equation (25) can be expanded and put in the form
I(%)-Z(w,) =
fl{;j-jn@ - Wfx, +R&J+ @Lx + &J2 + 2(1- u)R,$] dx dy
I
20.
Consider now the term u(R,, + R,J’, under the integral sign in equation (26). triangular and Schwarz inequalities or norms of functions, we find
OS
ss R,
(26)
By the
@,,+K,,J2dx dy~C(IR,,l(n,+IIR,~vlln,l2=IIR,:xlI~,+2lIR,xxlln,(IR~vrlIn~ + 1IR,~#L.
(27)
Upon substituting equation (27) into (26) and simplifying, we find
Recall that the remainder (R)i is of fourth order in x and y, i.e. Ri=f(E4, . . . , J4),
(29)
The maximum value of Z and 7 is of the order L/n (length where 3=x-x, and J=y-y,. of domain side divided by number of elements along the side). The derivatives of Ri are R i,,=m2,
* * - 3 J2)
RSyy =h(9,
. . . , Y2)
Ri,xy =_&(I’, . . . , 7’).
(30)
Theoretical and Numerical ConvergenceRates
1193
Hence,
2=j (Rx,)’ dQ=&d~,. . . , J4)dQ, IIn IIR‘XY i-l s
(31)
where the subscript i on R and R has been dropped for convenience. We can see that the above intergrands approach zero as n+co at a rate of (l;/n)4. Hence from equation (28) we have (32) and on redefining the K, we have
01(A~),-j2,1zw4,
(33)
where n is a measure of the number of elements along one side (or the mesh size). Accordingly, the convergence rate of the eigenvalue obtained from this element is nm4. This is in agreement with the value determined numerically by Lindberg and Olson [l I]. Static problem For completeness an analysis similar to that above will now be performed for the potential energy convergence rate in a static type problem. For convenience, we choose one of the simplest elements, a 6-degrees-of-freedom triangular in-plane element with trial displacement functions (U)i’c(r + a2x + a3y (34) To show that the second variation goes to zero and hence convergence with reduction of element size we proceed as before. That is, we first introduce comparison function in the form of linear finite Taylor series expansions (in s and y) of the exact solutions u and v about a point 0 within the ith element, i.e. (“>i=u10
+
(x-~xO)u~xJO +
(Y-YOb4ylO
(v)I=v(o+(x--so)v~,~o+(Y-Yo)v~~~o.
(35)
Let (R,)i and (R,), be the second order remainder terms of the expansion, i.e. (R,)i=3[(-u-XO)2U,,,IB+
* * * +
(Y -YO)2u9yylbl
(R”)i=~[(x-x,)21’,,,l + . * - +(Y-Yo)~v,~~ITI*
(36)
1194
GEORGE E. RAMEY and NATARAJAP* KRISHNAMURTHY
As before, in equations (35 and 36) we assume that up through the second order partial derivatives of the exact solution exist and are bounded (hence no stress singularities). The linear Taylor series comparison functions for the entire domain are
Among all the complete linear functions in .Y and y, U* and Y* minimize the potential energy of the problem. Hence we conclude that
or, in other words,
where
and
We now use the upper bound to S*rr({,, 5,) to show the convergence and convergence rate for the element. Based on equation (8) in Section 2, we have
From the definitions of u’ and < we see that
(41)
Hence, equation (40) may be rewritten as
OS i l-1
4CllR..,II’+2ll%II Ik.Il+
{ ~[ll&,,J)*+llRo,,ll’J+
IR...i12-$ (42)
Recall that the remainders (R,,), and (R,), are of the second order; hence, the first derivatives of these expressions are linear in x and y. Therefore,
Theoretical and Numerical Convergence Rates
1195
(43)
Hence, from equation (40) we have Osn(ti,
i;)-I[(u,
v)=~~J&, &
KII-~
w
or Oln(u*,
v*)-n(u,
v)~Kn-~.
(45)
Accordingly, the asymptotic convergence rate of the potential energy obtained from this element is Knv2. If one is confident of convergence when employing a particular element and is interested in the rate of mean convergence, the procedure for obtaining this rate may be greatly simplified. For example, in the above case, (g-u) and (g-v) are of order (n-‘); and the first derivatives of these quantities are of order (n-r). Since the second variation, and hence Aii, involves the square of these first derivatives, then this quantity is also of order (nW2). The area of the domain remains fixed as the mesh is refined; and ,therefore, the summation over all the elements leads to the conclusion that Aii (which bounds Arc+) is of order (PI-‘). Similarly for the thin plate rectangular bending element discussed previously, (r? - w) was of order (ns4); and, therefore, the second derivatives of displacements were of order (nW2). Since the strain energy expression involves the square of the second derivatives, then this quantity and hence the approximate eigenvalue is of order @I-~). In general, the asymptotic convergence rates of the eigenvalues and energy quantities may be regarded as Kemp. It can be noted that the partial derivatives of the remainder terms, (Z& are at least of order (n-l) for convergence and that the square of the norms of these derivatives will be at least of order (nm2). Hence, in the general formulation of the convergence rate, i.e. Knmp,p will be greater than or equal to two. Exceptions to this arise when there are unbounded higher order derivatives of the exact solutions which are exhibited by stress singularities. From the theoretical convergence studies, we see that the rate of solution convergence depends on the particular finite element used, the character of any singularities possessed by the exact solution; and it is independent of the mesh pattern This is assuming the problem’s essential boundary conditions are satisfied, and as previously, mentioned, neglects idealization errors. Theoretical convergence rate analyses similar to the two demonstrated here were performed for various other conforming elements. The results of these analyses are summarized in Table 1. The trial displacement functions for the elements considered are given in the Appendix.
GEORGE E. RAMEY
1196
and NATARAJAN KRJSHNAMURTHY
TABLE1. Conforming displacement finite element energy and eigenvalue theoretical convergence rates
Symbol
Type element
Element shape
T6 R8 RP16 CT20 R16 Tl8 W? LDT TT3 lT6
In-plane In-plane In-plane In-plane Thin plate Thin plate Thin plate Thick plate Torsion Torsion
Triangle Rectangle Rocttigle Triangle RectangIe Triangle Quadrilateral Triangle Triangle Triangle
Element degrees-offreedom 6 8 Y 16 :i 9 63
Theoretical convergence rate Kn-2 G-2 G-4 Kfl-6
iyn-4 Kn-6 Kn-2 Kn-2 Kn-2 Kn-4
References D4, I61 WI f:;
IS1
[IO, 11, 191 (81 WI [S, 131 W, 181 U6, 181
4. NUMERICAL INWSTIGATION Numerical results for the elements listed in Table 1 were studied with the emphasis being on the rate of solution convergence. For the most part, published numerical data were utilized and log-log plots of solution error or relative error versus mesh size were made for the various elements. For those cases where exact solutions were available the error terms shown plotted are absolute errors. Relative errors are plotted for the remainder of the cases by estimating (from results available) exact solutions and using these estimates to deermine the relative error. The format of these plots is identical to that presented by Cowper et al. [8]. It provides an excellent graphical means of checking tbe numerical verification of the theoretical convergence rate. In general, the error plots tend to be asymptotic as n increases. The siopc of the asymptotic line is the convergence rate which is easily compared to the slope of a straight line representing the theoreti~lly determined rate. For convenience of discussion the elements have been categorized as: (1) in-plane elements, (2) plate bending elements, (3) torsion elements. A separate plot for each element is presented in order that an unencumbered display of the numerical convergence rate for each element can be easily compared with the predicted theoretical rate. In all figures presented herein, the theoretical convergence rates are indicated by the dashed lines. Lighter dash-dot lines are provided in some instances for convenience in assessing the numerical solution convergence rate. In-plane elements
Error or relative error plots for in-plane elements are shown in Figs. 2-5. The elements considered are: 26, R8, RPl6 and 0’20. Several of these elements are applied to a variety of problems. For the cases considered, Figs. 2 and 3 offer good numerical evidence of verifying the theoretical convergence rates of fi’” and Kn-4 for the T6 and R8 elements respectively. Figure 4 indicates, for the small n values utilized, a numerical convergence rate for element RP16 of approximately 3.0 rather than the theoretically predicted rate of 4.0. Similarly, Fig. 5 indicates a numerical rate for element CT20 of 5.0 as compared to 6.0 predicted by theory. For higher n values, elements RR16 and CT20 may show results which are in closer agreement with theory.
Theoretical and Numerical Convergence Rates
r(Blaua
U!DJ+S
IJ!
anl+ojau
JOJJS
R
Tx
c
\
‘,
;;
tt6JaUe
U!OJ+S~
U!
JOJJ8
eA!+Dlau
c
1198
GEORGEE. RAMBYand NATARAJANKRISHNAMURTHY
j “:
0
,
,
n
0
”
0
'0 -
b d6J8US
U!DJ+S
U!
1:
,
JOJJS
SA!+DlS#
a
'0 -
1199
Theoretical and Numerical Convergence Rates
Comparing Figs. 2-5, one can observe the superior convergence rate of the higherorder elements. However, in the case of the RP16 element, it is doubtful that its slightly improved convergence rate relative to the RS element will justify its additional degrees-offreedom. Hence, for a given n size the R8 element will probably be more efficient than the higher-order RP 16 element. Plate bending elements
Error or relative error plots for plate bending elements are shown in Figs. 6-9. The elements considered are: R16, T18, CFQ and UK Both static and eigenvalue problems are considered for a variety of load and boundary conditions. For small n values, the numerical curves vary depending on the particular problem. However, most of the error plots indicate a numerical convergence rate that is in agreement with the theoretical rate for that element.
IO
I
k
10.
L_
1
I
I d
7
fvleosure
of number
, c
skm
of elements
FIG. 6. Relative error of finite element solutions-thin
per
side,
n
plate bending W6 elements.
1200
GEORGE E. RAMEY and NATARAJANKRISHNAMLJRTHY
The buckling coefficient for the square plate with a free edge shown in Fig. 6 exhibits numerous fluctuations such that a trend in the numerical results cannot be ascertained. The fluctuation of the buckling coefficient for the clamped square plate shown on that same figure is so great that non-monotonic convergence is indicated. This is understandable since it occurs at an n value of three which is not in the n=2, 4, 8, 16 etc. sequence. For these two buckling coefficient curves, it is believed that with increases in n, the numerical rates will be in agreement with theory as are the other curves on Fig. 6.
SYm I-Clomped
square
(n =4
shown)
\ \\
loadlpattern
Measure
Fro.
of
number
of
elements
7. Relative error of iinite element solutions-thin
per
side,
n
plate bending TlS elements.
Numerical results taken from [8] and shown plotted in Fig. 7 indicate the excellent convergence characteristics of the T18 element and offer good numerical confirmation of the theoretically predicted convergence rate of fi -6. Cowper et al. [8] concluded from their numerical study that solution accuracy depends in a large measure on the element arrangement (i.e. mesh pattern). Their numerical data indeed indicate this; however, if one makes an error plot utilizing their data, it can be seen that the rate of solution convergence for the problems considered does not depend on the mesh pattern. This can be seen in Fig. 7 where results for the two different mesh patterns P and Q are shown.
Theoretical and Numerical Convergence Rates
1201
. l
A,,-S. S. \ k-S.S.
S
\
plate
.
rectangular plate under pure shear
. \
\
\ \ k-Clamped square plate under biaxial compressIan \ \ \
.
.
.
\
.
\ l
.
\ k-S.S. square plate under A unlarial compression
\
0
\ \
.
\ \
\
’ \ “-2\ \ \
A- Frequency
\ \ ’
parameter coefficient
Buckling
‘\.
. \
k-
square
\
“\. \ \ \
1
I 2 Measure
of
4 number
of
I
i/I
6
8
elements
FIG. 8. Relative error of finite element solutions-thin
II IO 12 per
16 side,
20 n
plate bending CP’Q elements.
Figure 8 shows some eigenvalue errors resulting from use of CFQ elements. For the cases shown it appears that the order of convergence for the curves is Kn” which is in agreement with the theoretical rate The curve for 3,, 1 is converging at a faster rate; however, it shows fluctuations and the values of n for the curve are small. At higher n values, this curve will also probably converge at a rate of Kn- 2. Two of the buckling coefficient curves in Fig. 8 are for problems shown earlier in Fig. 6. Comparing these corresponding curves from Figs. 6 and 8, we can observe the superior convergence rate of the R16 element relative to that of the CFQ elements. The conforming thick plate element LDT exhibits a numerical rate in Fig. 9 which is m agreement with the Knq2 predicted by theory.
Torsion elements Torsional rigidity error plots for torsion elements TT3 and TT6 are shown in Figs. 10 and 11 respectively. In each case, the numerical convergence curve approaches the theoretical rate quite closely. It should be noted that it would be highly desirable to have
1202
GEORCXE.RAMEY and NATARAJANKRISHNAMURTH~
----
3
2
Measure
of
number
s s
square
Clamped plate
d of
5 elements
FIG. 9. Relative error of finite element solutions-thick
Plate
s~uore
6
7 per
8 side,
n
plate bending L DT elements.
at least one more n value data point for the 77’6 element to verify that element’s numerical convergence rate. The superior convergence rate of the higher-order TT6 element is evident from the figures; however, the additional degrees-of-freedom resulting from its use may result in its being less efficient than the lower-order TT3 element.
5. CONCLUSIONS The use of finite Taylor series expansions along with compatible form comparison functions is an effective means of determining conforming finite element theoretical convergence rates. These rates depend only on the particular finite element and are the same for both static and eigenvalue problems and are independent of the mesh pattern. This is assuming the idealized problem’s essential boundary conditions are satisfied. The theoretical convergence rate determination is seen to be invalid in the presence of stress singularities.
Theoretical and Numerical Convergence Rates
n
N
b
b Je+euoJod
r(+!p!b!J
IDUO!SJO+
U!
JOJJ3
1203
1204
Gsoxo~ E. F~AMEYand NATARAJANKRISHNAMUR~
In general, the asymptotic convergence rates may be regarded as Knsp, where p will be greater than or equal to two. From Table 1, one may conclude that in general as the degrees-of-freedom increase for a particular type (i.e. bending, torsion, in-plane) element the convergence rate also increases. The element convergence rates are of value in accessing solution accuracy and in comparing the relative merits of various finite element models. TABLE2. Summary table of element convergence rates Element
Theoretically predicted rate
T6
Kn -2 Kn-2
RP16 CT20 R16
Kn-4
R8
T18
CFQ
LDT lT3 TT6
Kn-6 Kn-4
Kn-6
Kn-2
Kn-2 Kn-2 Kn-4
Numerically determined rate Kn-2
Kn-2 Kn-3 K?l-5 Kn-4 Kn-6 Kn-2 Kn-2 Kn-2 Kn-4
Most of the error plots indicate a numerical convergence rate that is in agreement with the theoretical rate for that element. Relatively larger error fluctuation was exhibited for small values of n when comparing different elements, or when comparing the results for different problems while utilizing the same element. However, as n increased, the fluctuations decreased and the numerical curves appear to become asymptotic to the theoretical rate. A few curves exhibited so much fluctuation that a trend in numerical results cannot be ascertained. It is felt that for larger values of n, these curves will also become asymptotic at a rate in agreement with theory. A summary of the element convergence rates, theoretically predicted and numerically determined, is shown in Table 2.
REFERENCES [l] S. W. REY, A convergence investigationof the direct stiffness method. Ph.D. Thesis. University of Washington, Seattle, Washington (1966). [2] M. W. JOHNSON,JR. and R. W. MCLAY,Convergence of the finite element method in the theory of elasticity. J. Appl. Mech. ASME, 274-278, June (1968). [3] R. W. MCLAY.Completeness and convergence properties of finite element displacement functionsa general treatment. AIAA Paper No. 67-143, AIAA 5th Aerospace Sciences Meeting, New York, [4] ?P?z% A. P. BORESI,Kinds of convergence and improved convergence of conforming finite element solutions in plate bending. Nucl. Em Design 11,159-176 (1970). [5] P. P. LYNNand B. S. DHILLON,Convergence of eigenvalue solutions in conforming plate bending finite elements. Int. J. Num. Meth. Engng 4 (2). 217-234 (1972). [6] P. TONOand T. H. H. PUN, The convergence of finite element method in solving linear elastic problems. Znt. J. SW& Struct. 3 (5), 865-879 (1967). [7j J. E. WALZ,R. E. FULTON and N. J. Cvuus, Accuracy and convergence of finite element approximations. Proc. 2nd Con/. on Matrix Methoak in Structural Mechanics. Wright-Patterson Air Force Base, Ohio, pp. 995-1025, October (1968). [8] G. R. Cowpe~, E. Kosuo, G. M. LINDBERGand M. D. OLSON,A high-precision triangular platebending element. National Research Council of Canada Aeronautical Report LR-514, Ottawa, December (1968).
Theoretical and Numerical Convergence Rates
1205
(91 P. P. LYNNand G. E. RAMEY,A method to generate both upper and lower bounds to plate eigenvalues by conforming displacement finite elements. Proc. 1st Intertwtionai Conference on Structural Mechanics in Reactor Technology (1971).
[lo] F. K. BOONER,R. L. Fox and L. A. SCHMIT,JR., The generation of the inter-element compatible stiffness and mass matrices by the use of interpolation formulas. Conf on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 39743, October (1965). [ll] G. M. LINDBERG,M. D. OLSONand H. A. TULL~CK,Closed form, finite element solutions for plate vibrations. National Research Council of Canada Aeronautical Report LR-518, Ottawa, February (1969). [12] R. W. CL~UGHand C. A. FELIPPA,A refined quadrilateral element for analysis of plate bending. Proc. 2nd Conf. on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 399-440, October (1968). [13] B. S. DHILLON,Triangular finite elements for the bending analysis of thick elastic plates. Ph.D. Thesis, Dept. of Civil Engineering, University of Colorado, Boulder, Colorado (1970). [14] R. W. CLOUGH,The finite element method in structural mechanics. In Stress AnaIysis (Edited by 0. C. ZIENKIEWICZ and G. S. HOLISTER), Chapter 7. Wiley, New York (1965). [15] R. W. CL~UGHand J. L. TOCHER,Finite element stiffness matrices for analysis of plate bending. Conf. on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 515-545, October (1965). [16] 0. C. ZreMuEwrcz, The Finite Element Method in Structral and Continuum Mechanics. McGraw-Hill, New York (1967). [17] G. R. COWPER,G. M. LINDB~RG and M. D. OWN, A shallow shell finite element of triangular shape. ht. J. Solids Struct. 6, 1133-l 156 (1970). [18] G. E. RAMEY, A study of the convergence characteristics and accuracy of conforming finite element solutions. Ph.D. Thesis, Dept. of Civil and Environmental Engineering, University of Colorado, Boulder, Colorado (1972). [19] W. G. CARSONand R. E. NEWTON,Plate buckling analysis using a fully compatible finite element. AZAA J. 7 (3), 527-529 (1969). Received December 1973
APPENDIX Following is a listing of the code names and trial displacement functions for the conforming finite element models considered in this article. R16: 16-degrees-of-freedom
rectangular thin plate element
A plate bending element reported in [8, 11, 191, but generally credited to [lo]. (w*)i=cCl +a~x+a~~+a,.~2+a,x~+a,y2+a,x3+a,x2y+agxy2+aloy3+~11x3~ +a~2xy3+a13x2y2+a,,x3y2+a15x2y3+aIgx3y3. CFQ:
Clough and Felippa’s thin plate quadrilateral element
This is a plate bending element presented in [12]. Clough and Felippa first proposed the LCCT-12 triangular element with 12 degrees-of-freedom, and eaploying a compatible cubic displacement field. They further improved it to give the quadrilateral element designated Q-19 (referred to as the CFQ element in this paper), retaining the 12-degreesof-freedom for the element, the explicit displacement functions for which are very complicated and are not given here. T18 : Cowper et al. thin plate triangular element
A plate bending element reported in [S]. (~~*)i=a,+~2x+a3~+a,.~2+x5x~‘+a,~’+~,x3+a,x2y+agx~2+alOY3+a11x4 +qIlx 3~~+~sc,,x-2~2+a,4.~~3+cc,,~~4+a,,ssfa,,.u3y2+~ 1s.~2y3+y:19x~4+a20~s.
1206
GF_ORGE E. RAMEYand NATARAJAN KRKSHNAMLIRTHY
LDT:
Lynn and Dhillon’s thick plate triangular element
A thick plate bending element presented in [5 and 131.
T6 : 6-degrees-of-freedom
(W*)i=c(l
~azX+a,y+Cr,X2+~5~y+S.CLsS..’
(8*)i=c(7
+ ct8x +
agy
triangular in-plane element
The classical two-dimensional constant strain element given in [ 161 (U*)i=cIt f (v*)i’c14.
RS : g-degrees~f-freedom
CQX+ a,y
f c(gX + a6ym
rectangular in-plane element
The classical two-dimensional
plane element given in 1161.
(U*)i=tLi +a2x+a3y+a4xy (y*)i=a,+aex+a7y+awxy. RR16 : M-degrees-of-freedom
rectangular in-plane element
A higher-order two-dimensional plane element given in [16]. (u*) i=a~+a2x+a3y+a4x2+asxy+afjy2fa7x2ytagxy2 (v*)~=ag+aloxfa,,y+at2x2+a,3xy~a,,y2+a1S.~2y+a16_~y2. TT3 : 3-degrees-of-freedom
triangular torsion element
A first order (for warping function w) torsion element reported on in [16]. (w*>i=al +a,x+a,y. TT6: 6-degrees-of-freedom
triangular torsion element
A higher-order (for warping function w) torsion element given in [16]. (w*)j’ai CT20:
fazx+xSy+&$.~2
+a5xy+fC(6y2.
modified Cowper et aI. shallow shell element of triangular shape
In its ~rn~d~ed form, the element is a curved triangular element designed for shelI analyses and is reported [I?]. CT20 is the limiting case of a fiat plate, or the plane stress element contained within the total element. (u*),=a,+azx+a,y+a,x2+a,.uyfa,y2+a,x3t-a,x2y+agxy2+a,,y3 (v*)j=a,,+al2x+a,,yfa,,x2+a,,xy+a,~y2+a,,x3+ff,,x2yfa,~xy2+a2,y3.