ADAPTIVE FINITE ELEMENT ANALYSIS BASED ON P-CONVERGENCE
B. A. Szabo, P. K. Basu,
and
M. P. Rossow
Washington University
SUMMARY
The results of numerical experiments are presented in which a posteriori
estimators of error in strain energy were examined on the basis of a typical
problem in linear elastic fracture mechanics. Two estimators were found to
give close upper and lower bounds for the strain energy error.
The potential
significance of this is that the same estimators may provide a suitable basis
for adaptive redistribution of the degrees of freedom in finite element models.
INTRODUCTION
One of the most important problems remaining in finite element research is
the development of adaptive finite element software systems, i.e., finite
element computer programs which have a local error estimation capability and a
capability to increase the number of degrees of freedom selectively such that
the quality of approximation is nearly uniform over the entire solution domain
and the error does not exceed some pre-specified tolerance.
Research concerned with the development of adaptive finite element soft-
ware systems has been underway at Washington University for several years.
This work has resulted in the development of an approach for improving th.e
quality of approximation without mesh refinement. In this approach the number
of degrees of freedom is increased by increasing the polynomial orders (p) or
introducing non-polynomial basis functions over fixed finite element mesh divi-
sions.
This process of reducing the error of approximation through the addi-
tion of new basis functions is called "p-convergence" to distinguish it from
the conventional approach (called
"h-convergence") in which the size of finite
elements (h) is reduced while the number and type of basis functions for each
element are fixed.
The efficiency of p-convergent adaptive procedures has been established in
a series of numerical experiments,
reported in references 1 and 2 and it has
been shown that the rate of p-convergence cannot be slower than the rate of
h-convergence (I.
Babuska, private communication). In fact, for the vast
majority of practical problems,
p- convergence is substantially faster.
This and other computational advantages suggest that adaptive finite element
software systems should be based on p-convergence.
The key problem is to find suitable estimators of error which would indi-
cate when and where the number of degrees of freedom should be increased over
the solution domain.
Some estimators have been proposed already:
43
I -
I I II 1llllllI III I
-2-
1.
Babuska and Rheinboldt developed a local, asymptotic, a posteriori error
estimator for h-convergent approximations (ref. 3).
The estimator requires ele-
ment level computations only and measures the error in strain energy associated
with an element. An optimal distribution of the degrees of freedom is obtained
when this error measure is the same for all elements.
2.
Melosh and Marcal (ref. 4) proposed to measure the specific energy differ-
ence (SED), defined as the largest difference over the element domain between
the computed strain energy density function and the same function evaluated at
the origin of the elemental coordinate system.
This measures the effect of
modes higher than those associated with the (generalized) constant strain
states on the distribution of strain energy density within finite elements.
The criterion for mesh refinement is that SED be approximately the same for
each finite element. In practical computations SED is approximated by
the
larg-
est of the strain energy density differences evaluated at quadrature points only.
3.
Peano et al.
(ref. 5) proposed a criterion for p-convergent approximations.
This criterion is based on the rate of change of the total potential energy
with respect to higher order displacement modes, evaluated before the higher
displacement modes are actually introduced.
When the rate of change of the
potential energy exceeds a prescribed tolerance, the stiffness terms corre-
sponding to the higher modes are assembled and the new system of equations is
solved.
The hierarchic structure of the elemental stiffness matrices permits
efficient use of block relaxation procedures in obtaining improved solutions.
In this paper we examine the numerical characteristics of a criterion
similar to that proposed by Babuska and Rheinboldt in ref. 3, but modified for
p-convergent approximations.
Our study is preliminary in nature and is re-
stricted to one specific problem.
ERROR MEASURES
We have examined two measures for the error in strain energy on the basis
of a problem in two:dimensional elasticity containing a geometric singularity.
This problem is typical for a large class of problems in linear elastic frac-
ture mechanics (fig. 1).
The error measures were as follows.
First we define
the ith component of the residual vector, which represents the unbalanced body
force,as
r
i
=G;
i,jj
+ (A + G) ii +x
j,ji i
in which
44
c and h are Lame's constants;
U.
1
is the ith component of the displacement vector
computed by the finite element method.
The subscripts range over 1,2.
'i
is the ith component of the body force vector.
-3-
me of the measures,
to be called the "r-estimator," is defined for the kth
finite element as
(r
2
1
+ r22) dA
(2)
in which
a
is a constant,
to be determined by numerical experiments;
l?k
is the polynomial order of the displacement approximation
over the kth, element;
Ak
is the area'of the kth element.
The other measure, to be called the "t-estimator", is defined over interelement
boundaries and external boundaries on which tractions are specified, as the
square of the unbalanced tractions.
Specifically,
at the boundary of two ele-
ments,
the vector of unbalanced tractions is:
t+> = [S..
1J
(a)(s) - Gijcb)(s)] nj
(3)
in which
; (a>
ij
is the finite element approximation to the stress
tensor for the ath element;
n.
J
is the unit normal to the interelement boundary;
S
is the variable along the element boundary.
At external boundaries the unbalanced traction vector is the difference
between the computed traction vector and the applied traction vector.
When
displacement vector components are specified,
the corresponding unbalanced
traction vector component is zero.
The t-estimator is defined as
(t12 + t22) ds
(4)
in which
B
'k
'k
is a constant,
to be determined by numerical experiments;
is the polynomial order of the displacement approximation
at the kth interelement or external boundary segment;
denotes the kth elemental boundary.
45
I -
-4-
Both the r and t measures were found to give close indications of the
total error in strain energy in p-convergent approximations.
The details are
as follows.
The sample problem represented in fig.
1 does not have a known exact solu-
tion. For this reason it was necessary first to estimate the exact value of
the total strain energy U. This was possible by utilizing the asymptotic
relationship given in reference 6:
IJ=up+c
NDF
(5)
in which
U
P
is the computed (total) strain energy, based on pth
order polynomial approximation;
C
is a constant;
NDF is the (net) number of degrees of freedom.
Extrap?l?ting on the basis of eq.
(5), U was found to be 0.7702 02e2/E
(i- 0.0002 CT R /E) in computations involving two different finite element mesh
dFvisions. The computed dimensionless strain energy values, the corresponding
errors and the values of the two error estimators for parameter values a = 3,
8 = 2 are given in Table I.
The percent changes in the estimators as p is in-
creased tend to bound the corresponding percent change in the strain energy
error with increasing precision such that the change in the r-estimator is
smaller and the change in the t-estimator is larger than that of the strain energy
error. Furthermore,
the values of these percentage changes are monotonically
decreasing for p > 2.
This suggests the possibility that two constants cr and
ct could be found such that a relationship,
C
' Rk(3) < ' - 'P ' =
'k -
C Tk(2)
tk
remains valid for all p-values and that the upper and lower bounds become pro-
gressively closer to the strain energy error as p is increased.
This is, of
course, highly speculative at the present because no theoretical justification
exists,
but consistent with the observations of this numerical experiment.
For
example, if we choose c
= 0.7008 and c
= 0.2657, the two estimators will give
the value of the strainrenergy error attp = 7.
The resulting relationship
between the strain energy error and the two estimators is illustrated in
fig. 2.
The question naturally arises whether the same estimators would bound the
energy error at the element level as well. Clearly, this approach will be use-
ful only if the indicators tend to zero with increasing p at about the same
rate as the error in strain energy does,
not only for the entire solution
domain but for individual finite elements as well.
The presently available
information is sufficient only to indicate trends.
46
-5-
The problem chosen for study does not have a known exact solution; thus
the energy error cannot be computed with precision.
For the global solution
the rate of convergence formula (eq.
5) provided a basis for extrapolation to
the limit, value of the computed strain energy values.
Such formulas are not
available for predicting convergence at the element level.
For this reason an
ad hoc procedure for extrapolation had to be devised.
It was assumed that the
strain energy converges at the element level as
$4 =
3
(k) + 1-1
PV
where u and v are constants, u
0-d
up 04
is the strain energy of the kth element, and
is the strain energy of the kth element computed on the basis of uniform
pth order polynomial approximation over all of the finite elements.
Taking two
consequtive values of
pol
7
nomial orders q and p,
u can be eliminated and the
extrapolated value of u(~
is
v (k) _ qvu (k)
$4 = ’ %
q
PV - qv
where
v
and u
(k)
were chosen such that for q = 5,
p = 6 and q = 6, p = 7 the
value of u(k) was constant.
The resulting estimate of the strain energy for
element 2 is u(k) = 0.0844 02E2/E.
When the p-distribution is uniform over
the entire mesh, the estimators vary over the finite elements by several orders
of magnitude.
As could be expected, their value is the greatest for the crack
tip elements (element numbers 1,3,4) and least for the elements remote from the
crack tip (element numbers 5,6,8). In those elements which are not on the
crack tip,
the estimators apparently approach zero faster than the error in
strain energy.
den the distribution of p is altered such that the element to element
variation in the estimators is reduced,
then the estimators tend to become
closer to the energy error.
Specifically, letting p = 7 for the crack tip
elements and p = 3 for the remote elements,
the variation in the estimators is
reduced somewhat but is still 5 orders of magnitude for the r-estimators and 4
orders of magnitude for the t-estimators. In this case we have complete third
order polynomial approximation in the transition elements 2 and 7 with some ad-
ditional shape functions ranging in order from 4 to 7.
The estimators for ele-
ment 2 were computed as crR2(3) = 0.2 x 10V5, ctT2(2) = 0.2 x 10D4, and the
dimensionless strain energy error as
0.4 x 10-4.
Thus the indications are
that the bounding property of the estimators may be preserved at the element
level only if the p-distribution is such that the estimators do not vary by
more than 2 or 3 orders of magnitude.
In all computational experiments conducted to date the r- and t-estimators
were found to be consistent indicators of the source of energy error, in this
case the crack tip singularity.
47