Error estimation and adaptivity in elastodynamics

Computer Methods in Applied Mechanics and Engineering North-~~~~~~~

1.01 (IW2) 369-395

CMA 285

The paper discusses adaptive procedures for etastodynamic problems using semi-discrete finite elements. The adaptivity in space may be done by mesh refinement fh-versinn), by increase of order of the a~~~ox~rnat~~n po~y~~mia~s f p-versionf or a combi~~tjon of these ~~~~-~ers~o~~*In time, the integration is either made by mode superposition of characteristic functions or by direct integration. The adaptation in time may either be made by change of the giobai time-step (h-version) or increase of the urder of ap~rox~m~~~on ~I~nornia~~ (p-version). Crwial parts of the analysis are the error estimatw and the error indicator, which are based either on interpolation theory or on tots! energy obtaincaS from ~stprocessed stresses. For the mesh generation and regener~~inn~ an isoline technique is utilized by which the mesh can be smoothed by inter~oi~~~n theory. The isolines are cre&ed for some characteristic scsk function. For time discretization, a new simple a posteriori local error estimator for the Newmark scheme is described. It is derived by a post-processing technique without soh&g new equations so the addit~o~~~ computational costissmall and the imp~em~~t~tion is convenient. The paper also describes a sequence of nested time-integration methods for dynamic probtcms formulated as two first order equations in time. These nested ~nte~ation methods in time based WI a hier~~h~ca~ ~~rmuiation in time are A-stable schemes of order 21%and L-stabk schemes of order 2% - 1. A considerable ~rnpro~em~~t in terms of accuracy as we11 as effectiveness is obtained, compared to currently available methods.

The ultimate goal in adaptive eiement prccedures is [I] to effkiendy and e~~~~rn~~a~~yfind a d~s~r~t~~~~~~~~ so that the dis~ret~~~tj~~ errur is g~~~~~lywithin a given toferance, and focaffy is uniformly distributed, Le. to find a so-called optimal discretization, A reliable estimation of the dis~r~t~~~ti~n error and an efficient ~~~tstr~~ti~n of the dis~r~ti~~ti~n are cruciaf issues for the dev~~~prn~nt of the ~d~p~~~~ finite clement pr~~~d~res. Despite that, ~u~sid~rab~~ effort has been devoted to these issues and, indeed, signikmt advances have been achieved for some pr~~~~rns~ these issues are still ~~~~~le~ging= During the early days of finite etements, the errors were estimated from the jumps of the stresses between the efements. Where the jumps were tao high the mesh was refined. Possibly also, the stress Geld was smoothed by averaging. The mathern~t~~a~ analysis fez effiptic probtems was studied with the assumption that the mesh was refined in a regular fashion. Asymptotic error estimates were then obtained which gave i~f~~m~tion about the rate of

370

N.-E. Wiberg et al., Error estimation and adaptivity

convergence and the dependence upon the regularity of the true solution, but it could not be used directly for a practical error estimation. In the early 198Os, two new concepts were developed which were both of significant importance for the subsequent development of adaptive methods. One was the p-method with finite element functions of hierarchical type allowing an easy improvement of the solution which also gave an a posteriori estimation of the error (see [2-41). The other idea was introduced by BabuSka [5]. He calculated the residual error consisting of the error in the equilibrium at the interior of the elements and the interelement errors and used this for an adaptive h-method refinement to obtain a small and even error over the region. It was later found that a proper combination of h- and p-refinement gave the best convergence [S-8]. A very important step forward was taken in 1987 by Zienkiewicz and Zhu [9] when they found a very simple way to obtain both local and global error estimates for linear elements. For problems in elasticity they found that the exact stress fields could be approximated with surprising accuracy by smoothing the discontinuous stress fields obtained from the finite element solutions. With the aid of an efficient generator for unstructured meshes, they succeeded in developing an economic adaptive procedure of high accuracy. They later extended this concept to higher order elements [lo]. Another important development was started in 1988 by Eriksson and Johnsson [ll], who transformed the a priori asymptotic estimates to a posteriori form. However, the estimates include a constant that is problem dependent. Its determination requires some effort. This estimate can with advantage be used for error smoothing to construct optimal constructed meshes [12]. Extensive studies have been made for the h-, p- and hp-versions for static problems. A summary of different procedures is found in [.5]. The h-version redefines the mesh either by a remeshing (a new unstructured mesh) or by a refinement [12]. In the p-version, the increase of polynomial order is made by use of hierarchical functions to be applied in space [3] as well as in time [4, 131. As fine meshes for the h-version are numerically unstable, the p-version has no dramatic increase of condition number which makes iterative solution procedures for the system equation possible. For the dynamic problems, few studies have been made for h, p and hp in space and h and p in time. In [14], Friberg et al. propose an error estimate and an adaptive finite element procedure of the p-version for eigenpair computations. In [15], JOO and Wilson study structural dynamic problems by a mode superposition of load-dependent Ritz vectors. They use the error estimate of residual type directly as an error measurement for each Ritz mode needed. Refinements are performed until the errors calculated for all the necessary Ritz modes are small. Then they perform the mode superposition based on the final mesh to compute transient responses. The accuracy control of time integrations in dynamic analysis has been studied for some years. Park and Underwood [16] proposed methods for the time-step size control in the central difference schemes, where the step size is predicted in terms of an ‘apparent highest frequency’. Following their work, Bergan and Mollestad [17] proposed an automatic timestepping algorithm in the implicit integration schemes. To estimate the local error when a single step method is used, a traditional method is to compare the results when two different step sizes arc used, or to compare the results given by two integration methods of different order [18]. Recently, Zienkiewicz and Xie [19] presented a simple local error estimate and an adaptive time-stepping procedure for the integration schemes of Newmark type, in which the

N.-E.

371

Wiberg et al., Error estimation and adaptivity

local error estimate is derived using a Taylor series. The adaptivity in the time domain can also be performed by a hierarchical p-version in time [4, 181. A comprehensive study of error estimation and adaptivity in FE analysis of static and dynamic problems has been performed by Zeng [20], Wiberg and Zeng [21], and Zeng et al. [22, 231. Recently, adaptive procedures using space-time finite elements have been proposed for some parabolic and first order hyperbolic problems and very interesting and promising results have been obtained [24]. For elastodynamic problems, Bajer [25] studied space-time finite elements with moving spatial nodes in two dimensions, in which he used the error estimate of interpolation type as an error indicator to change the locations of spatial nodes from time to time to simulate the motion of steep stresses.

2. Basic equations and FE formulation for dynamic problems Elastodynamical

problems

are governed by the differential

p$+cpKdh=J(t),

u = u(x, t)

)

equation

x = x, y,

2

3

(1)

where u is the displacement, f the strain operator, Kd the elasticity matrix, c the damping, p the mass density, f the time-dependent body force and t the time. By the approximation

2(X, t) = MU

= u(x, t)

(2)

)

with basis functions 4 in space and nodal parameters U and applying the Galerkin procedure, the weak form of (1) leads to a linear system of ordinary differential equations of semidiscrete form Mij+Cti+KU=F,

0(O) = tie )

U(0) = U, >

(3)

where M is the mass matrix, C the damping matrix, K the stiffness matrix and F the forcing vector. The integration in time is made either by mode superposition or by direct integration. Mode superposition is based on the solution of an undamped system C = 0 with harmonic behavior U(t) = eiw’fi ,

giving the eigenproblem (K-hM)&=O,

A=6J2.

For the direct integration, we use the Newmark equivalent set of first order equations. The Newmark scheme is defined by

(4)

scheme or higher order schemes

for an

372

N.-E.

Wiberg et al., Error estimation and adaptivity

MO,+,.,, + ~i(+~, + Ku,+~, = FI+Al , tit+,, = 0, + [(I -

ut+A,,=

U, + ir, At + ii; At2/2,

iif” = [(l - 2p)ii,

Y

)ii, + Y~~+~~I At ,

+ 2p0c+A,] ,

@a, b) (5~3 4

where At is the time-step, p and y are two parameters, U,+Ar, ot+Ar and o,+A.r are the approximations of U, ii and 0 at time station t + At, respectively, and 0: is the approximation of the acceleration. Time integration can also be made on an equivalent set of first order equations. We thus change the second order equation into two first order equations in time,

or j,+Ay=f,

W)

Y=[y.

If we make an approximation r(t) = 4(t) y

in the time domain (7)

7

and use a weighted residual approach with weighting functions &, we obtain a time integration method of the type TY-j,

Tii =

+ A4)j dt I,:$ii.c(i,

3

f”= j$;f.dt.

(8)

Different adaptive procedures for reducing the spatial discret~zation error can be obtained by changing the mesh (h-version), changing the degree of polynomial (p-version) by hierarchical functions or by the combined @-version. Static problems are obtained when M = 0 and C = 0. Similar procedures of h-, p- or @-version in time are also possible for reducing the time discretization error.

3. Spatial discretization 3. I. h-version

The h-version means a rede~nition of the mesh using the same element types. The procedure is straightforward and frequently used for static problems. For dynamic problems, the interpolation of function values (and its derivatives) between successive meshes is very important [26] and must be taken into consideration. The solution procedure is very much dependent on an appropriate error estimation and an effective mesh (re)generator [12]. 3.2. p-version For the p-version,

the basis functions are selected in a hierarchical

way so that

N.-E.

Wiberg et al., Error estimation and adaptivity

where & denotes basic basis functions, &, the hierarchical basis functions, and U, and U, the corresponding nodal variables. For hierarchical formulation in space, (3) obtains the partitioned form if C = 0,

The system of equations

corresponding

to the p-version

for static problems,

U, = K;;F,

(direct) ,

K,,U,

WG>

= F,, - K,,U,

WV

>

was due to the positive-definiteness and the limited condition number [4] solved by a combination of a direct method for the first equation and a preconditioned conjugate gradient (PCG) method for the second equation. For the p-version of the eigenvalue problem, (4) has the form [14]

The system matrix is not always positive-definite a direct method.

or well conditioned

so (10) must be solved by

4. Time discretization 4.1. h-version The h-version means that the global time-step At is changed. The integration Newmark scheme (5) giving

m+a, =C,A, 7

l? = K + M//3

At'

+ yC/p

At

,

F=

F(P,

is made by the

y, U,,

ir,,ii,,At). (11)

4.2. p-version For the first order system (6), we now increase the order of the approximation in time by adopting polynomials of higher order in order to reduce the local truncation error. The p-version of hierarchical form (9) can also be applied in the time domain [4, 181, in a hierarchical way so that the basis function &(ti+,) = 1 and 4(t,+,) = 0 for i > 1, where 4j is a polynomial of order i. For convenience, the following presentation is restricted to discretizations such that i < 4. A possible set of base functions is given by [18]

374

N.-E.

$1 =

0 3

sb, = im(f

& =

Wiberg et al., Error estimation

e(1 - e> )

- e)( $ - e)( j - e)

and adaptivity

6 = eu - w:

- e) , (12)

,

where 8 = r/At and r is a local time variable. By use of the basis functions, we can expand y as Y = Mt)Y=

?Cbh!Yb + &Y, = &Y, + &Yz + #%Y, + 64Y, f

Yb = y, f

03)

where Yb contains the unknown displacement Y1, connected to the basic linear approximation, at the end of the time interval, and Y,, contains the hierarchical variables. We see that Y, represents the increment of the solution over the time interval. Subsequently, we select four weight functions +i so that four equation systems are obtained from (8). These systems may be assembled to a single system: TY=f,

T, =

+:(I& -+-A+,) de 3 i,j=l,2

,...,

pn,

(14)

We observe that each node variable in space is associated with m-variables in time in each time-step. So in the general case, T may be partitioned into 16 N X N matrices, where N is the number of spatial variables. 4.2.1. A-stable schemes We now select weight functions $i. We note that the highest order of accuracy for m = 1 is 2 which is featured by the Crank-Nicholson method, which may be derived with #i = I. In an attempt to recover the properties of this method, we therefore make the same selection here; the remaining functions are selected in a Galerkin-like fashion. Thus, we choose (l/i =

f3-l ,

(15)

i=1,...,4.

This set of weight functions gives the coefficient matrices 0

T=

-&+;A

-&+$A

1 40h

Tm.m9m=3 (16)

shoufd be emphasized that the system (14) describes a ‘hierarchicaf’ apprux~mat~on in the sense that the vector Y, (cf. (13)) alone gives the irxrement of the sofution between times b, and ti+l. The remaining inefficient vectors Y,, j > 1 may be includes if one wants to increase the poly~~~~a~ order of the time discretizatron (13) of the unknown vector y, The outlined weighted residual ap~~~~b is thus anafaguus to so-called hiera~~~j~~~ ~-re~n~m~~t methods that have been derived for space d~s~~~ti~at~uns of ~~~~~~~pti~a~ problems [4]. ff the concept is utilized in conjunction with a simuftaneous space-time discretization approach, it should thus be possible to obtain a hierarchical finite element for~~~~lation for r~~~~rnents in both time and space domains. By studying, over a time interval, the ~rnp~i~~atio~ factor a(z) corresponding ,4, we find that a(z) is the Pade appro~~~~t~~n P,.fn for P- ‘, and the t* LJ8 form=1,... method ob~~~ed is A-stable and of order 25~ (see [I$]).

It

We derive t-stable rn~~~~ds in an analogous rna~~~~. We note that the t-stable backwarddifference method is obtained for a linear dis~retj~~tion of y in conjunction with a Dirac d~str~but~~~~&(I - 8) rentred at B = 1 as a weight fu~~ti~~. Hence, it might be feasible to apply the weight functions

t/jr= S(1 - e> 1

iL;“P$

i=2,3,4.

(17)

We now can establish and solve (14) fur different vafues of m = 1 3 . . + 9 4, As for the A-stable methods, the respective ~rnp~i~~ati~~ functions a(z) can be determined_ The ~~s~retizatio~ implied by (12) in conjunction with the set of weight functions (X7), result in a sequence of nested integration rneth~~s that are L-stable and of order 2m - 1.

The error can either be fcrr the fuu~~i~~ itself or its derivatives e=

ak-

M ~

~~~~~~~~~~~-~

L

w

The error comes both ,from the spatial discretization and the time ~is~retization~ This means that for an ~ptirna~ solutions both the spatial mesh and the time-step ought to be refined in an adaptive way, based on global and local error control,

The spatiat discretizatian error is either estimated by use of i~t~~~~~at~o~ theory or by energy considerations. In interpolation theory [lltt an a priori error estimate in the maximum norm of the spatial discretization with mesh size h using linear trianguiar elemems is

In dynamics it might be expropriate

to Control the analysis by studying the error in energy.

376

N.-E. Wiberg et al., Error estimation and adaptivity

For these problems, the total energy E, is written as the sum of strain energy E, and kinematic energy E,, (20a), and the error of the total energy AE, is denoted as a sum of the error in strain energy AE, and kinematic energy AE,, (20b) , E, = E, + E, = &[a@, u) + (ci, pci)] ,

AE, = AE, + AE, .

(2Oa, b)

By introducing first order approximations of uh it can be easily shown for lD-problems that AE, S C,h2 and AE, d C,h4, where Cj are constants. We find that AE, is two orders smaller than AE, so AE, -AE,

,

(21)

which means that the same error estimation can be used both for static and dynamic problems

WI-

The interpolation theory can be used to smooth the error. We can, by adjusting the element size, make the right-hand side of (19a) almost equal for all elements. If isolines are created for, [VU] (or ]a]), the distance between the isolines are inversely proportional to the gradient of ]Vu]. By grading elements with sizes equal to these distances, an equally distributed error (19b) is obtained [12]. The energy formulation of the error can be used both for error smoothing and stop criteria. 5. I. 1. Error indicator, h-version In the h-version, an error estimate of postprocessing the energy norm for linear triangles,

ME = (1”e’?Kd?e

type is applied by studying the error in

dV )“’ = (2 AE,)1’2 = [ 1” (ah - a)Kd-‘(ah

- a) dV]“’

, (22)

where u is the exact solution. Obviously, if a stress or strain solution obtained by postprocessing is of one order higher accuracy than the one obtained from the finite element analysis, then the error can simply be estimated by substituting this post-processed solution into (22). Let (T” be a stress field given by a post-processing on the same mesh as used for the finite element analysis, we then have an a posteriori error estimate

llellE(R)=

- a*)do (ah- a*)‘D-‘(CT’

112

.

(23)

This error estimate was first proposed by Zienkiewicz and Zhu [9], where they used a ‘projection’ technique to recover a continuous stress field from the discontinuous stress field given by the finite element analysis. This error estimate is asymptotically convergent and computationally practical and will not overestimate the error. For this reason, they suggest using some empirical correction factors to improve the effectivity of the error estimate. For linear triangular elements, they proposed the factor 1.3 and for quadratic triangular elements

317

N.-E. Wiberg et al., Error estimation and adaptivity

the factor 1.4. Our experience is that this is not needed if the iterations in the recovery procedure are driven far enough. An error estimate of this type provides not only an estimation of the global energy norm error but also a good evaluation of the local errors, which is of importance for efficient mesh refinements. Note that the global error (23f can be obtained by summing up element contributions, M

’ Ilell2,,,= c (Il4lcck,>*

(24)

k=l

where II4 E(k)represents

the local error of element k and M is the number of elements. can define a relative error by

We

(2% The post-processing of finite element solutions has been studied for some years and much work has been done with the aim of achieving a superconvergent approximation [27,28]. Among the proposed techniques, some simple methods that are widely used in practice are: evaluating the stresses at numerical quadrature points and determining their values at the displacement nodes by extrapolation; averaging and smoothing the stresses based on some projection technique; using iterative techniques to construct continuous stresses and displacements. The simplest one may be the so-called nodal averaging technique, which calculates the stresses at a node by averaging the stresses of those elements connected to this node and then interpolating the stresses based on the shape functions used for the displacements. 5.12. Error indicator, p-version The error indicator often used in the p-version of static problems (26a) is based on the change of potential energy when a new extra hierarchical basis function j is added (26b), rt =

= 2(AEj)* llel~~~~l~ll~ Y 7: = (fi - Ki,~~~)‘/~;,,

For the eigenvalue

problem,

GM b)

an error indicator similar to (26b) is used [14],

vi = ([z$, - A”Mi,,]~j”‘)2/(kj”)[Kj,j

- h,‘“‘Mj,j]) )

(27)

which approximates the relative change of an eigenvalue i when one new hierarchical function i is added. In (27), ki(“)is the modal stiffness and (Al”‘, c$~“‘)is the eigenpair i with y1degrees of freedom. 5.2. Time discretization error, h-method 5.2.1.

Error indicator, ~-version

Different

methods

for automatic

time-step

control have been described

in the literature.

378

N.-E. Wiberg et al., Error estimation and adaptivity

Generally, the global error is difficult to obtain, and a local error method is often used in order to estimate the time discretization error and to control the step size. Traditionally, the local error estimates were obtained either from the difference of approximate solutions between a full step and two half steps or from the comparison of the results by two different methods. They are often impractical because of the substantial additional computational cost. Considering computational efficiency and convenience of implementation, Zienkiewicz et al. [29] introduced a simple expression to indicate the local error. Unfortunately, it does not always give an accurate approximation of the real local error. Recently, Zienkiewicz and Xie [19] proposed another local error estimator for the Newmark scheme by comparing the numerical solution with the exact solution expanded in Taylor series. Recently, Zeng et al. [22] derived a similar expression in a more simple and intuitive way based on the so-called ‘post-processing’ technique. This error estimator is simple and asymptotically convergent for the original Newmark scheme, however, it is invalid for the integration schemes which assume that the accelerations vary linearly within each time interval. Furthermore we propose a new simple and reliable local error estimator which can overcome this problem [30]. It can be obtained either by the post-processing or the Taylor expansion. An error estimate of postprocessing type [22] can be obtained if we assume that the acceleration 0: is constant within each time interval. A one order better approximation is a linear variation of ii,

(28) A post-processed

Integrating reads

error can be written as

twice and letting T be At gives the local error estimate,

which in a certain norm

This is exactly the same estimate as given in [19], where the derivation was made by Taylor expansion. For the constant-average-acceleration method (/3 = l/4), we obtain the error ]]e”]] = (At*/12)](0,+,, - ii,) I]. This estimate cannot be used for the linear acceleration method (p = l/6) because UN is linear and ii* in (28) is then not a higher order approximation. For many parameters, the Newmark scheme is of second order, so the local error should be of order At’ and not At2 as in (30). In order to obtain a local error of order At3, we go back to the basic assumptions of the Newmark scheme. The Newmark scheme contains an approximation of the derivative of the acceleration which can be seen if we further expand (5c, d),

u

r+br

=

U, + ti, At + ii,

At’f 2 + ti; At3t6,

(31)

N.-E.

Wiberg et al., Error estimation and adaptivity

379

The approximation of 0: depends on an average derivative of the acceleration fif which keeps constant within each time-step. A higher order approximation is then a linear variation of the derivative of the acceleration during the time-step, ...

...

u;=Ut+;(u,,,,-ti,),

(33)

TEO,At.

Figure 1 shows the variation of the derivative of acceleration for the case of the linear acceleration method. Naturally, the derivatives of accelerations in the linear acceleration method ( p = l/6) keeps constant within each time interval, which means that t?f” = up. Now we construct a post-processed error of the derivatives of accelerations [29],

Since the values at time station t are assumed to be exact, integrating (34) three times and letting T be At, we obtain the local error estimate which is given in a certain norm either as (35a) or (35b), ((eNI =

IkNII

& At3(((12P - l)ut

= &

When /3 = l/6,

At3(1(6p

- f)tit+At

(35a)

- c,(] ,

- (;

-

6p) tit11.

(35b)

we obtain the error (36a) for the linear acceleration

IleLII =

method

as

& At311up - tit11.

(36)

...

u

,+&

_-------_-_----.

.I.

.

.

Ur = 6/Y&’q-----_--___//... u,.,

/ ----I

-

-

- 2 4

=

$(i;;+N + ii,)

1 I I

Fig. 1. The derivative of acceleration: linear acceleration post-processed continuous derivative (---).

method (-),

general Newmark scheme (- - -) and the

380

N.-E.

Wiberg et ai., Error estimation and adaptivity

When p = l/4, we obtain the error estimator (37a) for the constant-average-integration method, which for ii, = 0 equals the error estimate Ile”“ll in (37b) and for 6: = o{ equals the error estimate I/e”// in (37c). lIeAl) =

k At”/l20: -

Ile”ll =

A At"(ll?,+,, - ii,11.

tirll ,

IleoOJl= i At211~,+At -

ii,11,

(3% b) (37c)

As shown in [30], the simple local error estimator (35) may also be obtained either by Taylor expansion or by the difference between the solution of (5) and a solution based on the assumption (33). 5.2.2. Error ~nd~~~torp-version We define the error vector e as the difference approximation yh, similar to (18). As a convenient interval [ti, tit,], we use the &-norm

between measure

the exact solution y and the of the magnitude of e in any

(38%b) This error measure gives a ‘global’ estimate of the error within the considered time interval. The total error is obtained from a summation of contributions ei denoting the jth component of the error vector,

llejllL,= (I:” efdt]“‘;

llellLz= (i lle~lli2]“2. j=l

(3% b)

To estimate the error in the approximate solution, we consider two consecutive solutions, one of order k and the other of order k + 2, of the system (14). The difference between the two solutions is denoted i. We may expect it to be a reasonable approximation to the exact error e in the lower order method, as the higher order method embraces two more terms in the power series of the exponential expansion of y. If Sj denotes the jth element of the vector 6, one thus has the ‘local” error estimate

IlijllL, =

{ire’ 2:

dt}"*

= lIejIlL,.

6. Adaptive procedures 6.1. Error bounds

In an adaptive process in space and time, we need the error of some problem dependent parameters (of some important physical quantity) described by some norm //e/l, e.g. the

N.-E.

Wiberg et al., Error estimation and adaptivity

381

maximum norm, &-norm or energy norm. We can prescribe the tolerance of the error itself lje[l (24) or the error in a relative form n (29,

where y1 and y2 are two given parameters, E is an absolute error tolerance and ;i is a relative error tolerance. It is often appropriate to set yZ= 1. If the error bounds are not passed, we need in space to change the mesh size (h-version) or approximation (p-version) and in time to change the time-step length or the approximation (‘-version). To study the effectiveness of the error estimates discussed, an effectivity index is defined as @_

Il4””

0 > 1

Il4””

0 < 1

where the superscript ones.

’

(overestimation of error) , (underestimation of error) ,

(42)

‘ex’ denotes exact error or exact relative errors, and ‘es’ the estimated

6.2. Spatial mesh adaptation Suppose that we have accepted numerical solutions at time t on a spatial mesh. To find numerica solutions at time t + At, we use the same spatial mesh as the first trial spatial mesh. In our refinement strategy, we aim to generate an optimal mesh, on which the condition (41a) or (41b) is satisfied and the error is uniformly distributed on each element. We therefore define a local refinement parameter for element k on the current mesh, 5: = l/el/E(KJ k efn ’

5< 1, 5: =

1 ,

5>1,

enlarge element size, optimal size, reduce element size,

(43)

where Zm is a permissible local error given in an average sense for the current mesh. ff an absolute error toIerance is specified, ~7~is given by (44a) and the local refinement parameter 5, clearly provides information on Iocal accuracy achieved on each element and gives a convenient guidance for remeshing. For an element, if 5 > 1, the element size should be reduced and if 5 < 1, the element size should be enlarged. A mesh for which all elements possess 5 = 1 is an optimal one, in the sense that the local error in strain energy norm is uniformly distributed. If a relative error tolerance is specified, tYmis given by (44b).

2E,(uh, zi”) + ~lell~~*~ 1’2 e;, = +j > , M

(

where M is the number of elements on the current mesh.

Pa M

382

N.-E.

Wiberg et al., Error estimation and adaptivity

We also notice that, if the system possesses no dynamic effects, i.e. Ek = 0, (25) and (44) degenerate into formulae commonly used for the elastostatic problems. When a remeshing is necessary, we use the current mesh as a background mesh and predict new element sizes based on the local refinement parameter tk. Noting that the convergence rate is O(hP) and by a comparison of the error on the old mesh with size h, and the new ‘optimal’ mesh with size hi we obtain the prediction h;

=

t;liPh,

.

(45)

The element size predicted in (45) is given for each element on the current mesh. To obtain nodal values of predicted element sizes, we take the mean value of predicted element sizes around each node. In this study, linear triangular elements (p = 1) are used. When the position of the system is almost in neutral, the stresses are nearly zero all over the spatial domain and, thus, a coarse mesh can be used for the stress energy part of the problem. In such stages, the accuracy of the numerical solution depends mainly on the time integration, rather than the spatial mesh. Therefore, the change of the spatial mesh is not necessary and should be avoided for these stages. In the present procedure, we identify such stages by comparing the current strain energy Et with the maximum strain energy (Et),,, recorded during the computation. When the current strain energy is only a small amount of the maximum strain energy, i.e. (46) where (Y is a small threshold parameter, the spatial mesh is not updated if the accuracy requirement (41) is violated. Numerical experiments indicate that it is adequate to set cv = 0.1. In practice, it is also often necessary to specify both the upper and lower limits for the element size, i.e. to set (47)

The knaxis specified so that the element size will not exceed a generally acceptable size. The hmin is specified in order to keep the stability condition of a time integration scheme, and to avoid numerical difficulties caused by very small elements. Another important issue must also be addressed. As the spatial mesh is updated from time to time, a finite element solution accepted on a mesh for a time station, say t, has to be redefined on an other mesh newly generated for the next time station t + At. This process should be efficient and economic. To speed up this process, it is important to develop a searching algorithm which, given a point, finds in an efficient way an element on the mesh used for time t, where the point is located. An overall search is time consuming and should be avoided. In this paper, we combine a quadtree technique with a searching technique using area coordinates [26]. Remeshing between time-steps means that the displacement field becomes discontinuous, and energy dissipation will take place. This matter requires further study.

N.-E.

6.3. 6.3.1.

Wiberg et al., Error estimation

383

and adaptivity

Time-step adaptation h-version

An adaptive analysis should, ideally, find a discretization economically, such that the local error is uniformly distributed and the global error is within a given tolerance. For time dependent problems, however, the global error estimation is generally difficult to obtain. Adaptive algorithms for dynamic problems are usually designed based on the control of the local error in a time-step. Here, the aim becomes to adjust the time-step size or the p-level in the time-step in an efficient and economic way so that, for each step, the local error is roughly equal to a prescribed error tolerance F or 77as in (41), using the L,-norm. When condition (41) satisfied, the solution is accepted and the time integration proceeds to the next time-step without changing the step size. The right-hand side of (41a), i.e. y2F, is an ‘upper error limit’. If the estimated error is larger than this upper limit, the solution is not accepted and the step size should be reduced. The left-hand side of (41a), i.e. y,E, is a threshold, rather than a lower error limit. When the estimated error is less than this threshold, the solution is accepted, however, the step size should be enlarged before stepping to the next time-step. It is from the efficiency point of view that we do not set a lower error limit, but set this threshold instead. In practice, an implicit algorithm cannot be efficient if the step size varies too frequently, as the effective stiffness matrix has to be factorized when a new step size is used. Therefore, the step size should not be enlarged until ]]e]]< ylE has been registered consecutively for K0 time-steps [19]. A proper choice of this parameter K, can improve the computational efficiency. On the other hand, y1 could be set so K, = 1 and thus unnecessary. Usually, it is inconvenient to specify the absolute error (41a) or the relative error (41b), so we may define a relative error qm. (48)

wherewll>,ax is the maximum value of the corresponding norm of the displacement solution recorded during the computation. When a step size should be updated, the prediction of the new step size has to be made such that the prescribed accuracy can be achieved with the least cost. It is known that, for the Newmark integrations, the rate of convergence of the local error achieves O(At’) as shown above. Suppose that the current step size is At and an estimation of the local truncated error has been obtained. Then, we may write ]]e]] and ((e]]’ with a new step size At’; we aim at obtaining Ilell’ = C,(At’) = i ,

lIeI\= Cl At3,

namely, the local error is equal to the given tolerance. the new step size -

At

or

At’=

0 77 77

(49) If C, = C,, we obtain a prediction

for

113

At,

At’=

fi (

At. -

i”‘”

Pa, b)

To determine the d~strib~tj~~ of the error, we may utilize tke A,-norm (39a) of the individual components ei of e, The relation between the “local’ error (39a) and the ‘global” one is given by (39b). Whether an error should be regarded as acceptable or not, would depend on the size of the time increment h. For this reason, the scaled estimate

may be more ~onven~en~~ ~e~~~~~ng the ~~te~nd Zi in (40) by i’&? one realizes that G+.can be viewed as a root mean square of Gj; i.e. if the trnn~a~io~ error is (and indeed eczufd be) constant, then Ztzjis the estimated absolute value of the constant, Error Iimits could be given on ghj in the same way as in (41). Both overall and partial p-refinement can be handled.

Figure 2 shows a thick circular arch under plane stress. The p-method has been used to compute the eigenpairs based on three different h-refined meshes, [31J The hierarchical variables are activated for the condition qy+’ 2 y,,,ax(~m’l), where y is a threshold value, Table 1 shows for the meshes in Fig.. 2, the 1.5 lowest computed eigenvalues for the error ~~d~~~t~r threshold y = 0.1 and for two cases of the error ~~tirn~tor E = C positive of”” + Eto1=L:CL05 and ctOl= 0.01. The eigenvalue number is EN. The hierarchical functions (up to second order due to limited space) are introduced successively and the number of variables is NV when erErat is reached far the highest eigenvafues.

N.-E.

385

Wiberg et al., Error estimation and adaptivity

Table 1 The 15 lowest eigenvalues of the thick circular arch EN E,‘,,

=

0.01

Mesh 3

Mesh 2

Mesh 1 *

Et”, = 0.05

E,,, = 0.01

*

l1”, =

0.05

Et,, = 0.01

*

0.14676 1.68757 3.11368 7.53954 12.7882

0.14676 1.68757 3.11368 7.53954 12.7882

0.13699 1.65421 2.88571 7.37153 10.7278

0.13699 1.65350 2.76334 7.35959 10.5188

0.13651 1.65345 2.76246 7.35940 10.5188

0.13646 1.64573 2.74327 7.33142 10.3929

0.13502 1.64541 2.71044 7.32292 10.1631

0.13489 1.64534 2.70427 7.32092 10.1508

6 7 8 9 10

24.9755 34.5761 52.2911 61.7774 64.5273

24.9755 34.5761 52.2911 61.7774 64.5273

24.1940 25.6894 46.8722 50.5710 53.5118

24.0621 24.1669 43.7350 49.7954 52.1584

24.0507 24.1656 43.6960 49.7415 52.1102

23.7219 23.9960 41.1999 49.4859 51.5753

22.9994 23.9288 40.8180 49.1113 51.4715

22.9228 23.9117 40.5987 49.0290 51.1787

11 12 13 14 15

91.8795 96.7110 106.754 119.063 120.633

91.8795 96.7110 106.754 119.063 120.633

72.7535 78.8856 81.5250 105.830 106.148

68.3978 78.0354 80.3944 99.1808 102.390

68.2837 77.6855 80.3659 98.9454 101.694

65.1920 77.0520 78.2549 93.3309 99.5498

62.4010 76.5377 76.9337 88.0232 98.9789

62.3098 76.2547 76.4692 87.1718 97.5562

1 2 3 4 5

NV

48

48

140

158

160

464

534

576

* All hierarchical functions are activated.

7.1.2.

h-version, direct integration, stress wave propagation

Figure 3(a) shows a rectangular domain subjected to a half-sinusoidal pulse q(t) giving a stress wave propagation [23]. The mesh is updated from time to time based on an error estimation (23). Only the spatial discretization error is studied, so the time-step chosen is small. Figure 3(b) shows the predicted graded mesh for t E [3.1,3.3] . Figure 3(c) shows the variation of the temporal number of degrees of freedom. Figure 3(d) indicates that the prescribed accuracy requirement has been well achieved. It also shows that, although the number of degrees of freedom used in the second conventional analysis is about the same as the maximum number of degrees of freedom used in the adaptive analysis, the temporal error is larger for almost all stages than that achieved by the adaptive analysis. The average number of degrees of freedom used for the adaptive analysis is about one-third of those used for the first conventional analysis, i.e. on the denser mesh. The CPU time used for the adaptive analysis is about 94% of that for the conventional analysis on the denser mesh. These results indicate the good potential that the presented adaptive procedure possesses. It should be noticed that the efficiency can be improved by specifying a higher upper error limit so that the mesh is not updated too often. In Fig. 3(c), we see that the remeshings have been done quite frequently in order to monitor the moving of steep stress regions. 7.1.3. h-version: L-shaped domain subjected to a triangular pulse Figure 4(a) shows an L-shaped domain subjected to a triangular pulse [23]. Plane stress conditions are assumed. The dynamic response in a time interval of (0.0, 8.0) is computed. A

386

(b)

NEL=681

t E 13.43.31

r NDOF=663.

NDOF=

CC)

g,alo- - - - -

NDOF=940

2094 for 1st conventional

for 2nd conventional

analysis- - - - - - - - - - - - - - -

NDOF for the adaptive ,

4.0

analysis

analysis I

I

0.o

a0

I.0

llma

I-

adaptive

Ji *to= 25

-x (4

P

20

1

--1 -----

15 ifi

--L

analysis

, NDOF E [194.801]

the upper error limit :yz

* F = 8.8 x 10S3

conventional

analysis with NDOF=2094

conventional

analysis with NDOF=940

Fig. 3. (a) A rectangular domain subjected to a half-sinusoidal pulse q(t) = sin(nti4) and a sketch of its stress solution (the upper part). (b) Predicted mesh at time station t E [23]. (c) Temporal number of degrees of freedom for the adaptive analysis. (d) Temporal error I/PII~/v? for FE analysis.

N.-E.

Wiberg et al., Error estimation and adaptivity

387

relative error tolerance of G = 5% is requested with -yl = 0.8 and y2 = 1.4. The relative error defined in (3.21b) is again used to control the computation. The minimum and maximum element size is set to 1.0 and 50.0, respectively. The time-step size used for the Newmark scheme is set to 0.1. Meshes predicted in the adaptive analysis for time, t E [0.3, 0.51 and t E [5.1, 7.51 are plotted in Fig. 5(a) and 5(b). From the predicted mesh in Fig. 5(a), we can see the very strongly graded meshes due to the steep stress location. As time goes on, bending and shear deformations occur, making the solution relatively smoothed in Fig. 5(b). The temporal number of degrees of freedom used in the adaptive computation is given in Fig. 6(a). The temporal relative error is given in Fig. 6(b). The horizontal displacement response at point B and the stress response of the component aI at point C are plotted in Figs. 7(a) and 7(b), respectively. We observe that the responses computed by the adaptive analysis again agree well with those obtained by the conventional analysis made on a quite uniform dense mesh with 4394 elements. These indicate further that the adaptive analysis makes the solution reliable. The CPU time for the adaptive analysis is about 86% of that for the conventional analysis. The adaptive analysis is thus done in a reasonably efficient and economical way.

7.2. Time discretization 7.2.1.

Single degree of freedom:

The equation the form

error estimators

of motion and the initial conditions for an undamped

i+w2x=o;

0 = 1.571 )

50.0

x(0)

= x0 =

1.0 )

single-DOF

system has

X(0) = ii, = 0.0.

(52)

50.0

1

50.0

50.0

i a)

L-shape

domain

Fig. 4. An L-shaped domain subjected to a triangular pulse (plane properties are E = 1.0 X 104, u = 0.3, t = 1.0 and p = 25.0).

b)

triangular

pulse

stress conditions

are assumed;

material

388

N.-E. Wiberg et al., Error estimation and adaptivity

a)

NNO=904 NEL=l737 NDOF=1789 mesh

for

t E [0.3,0.5]

Fig. 5. Meshes predicted for time t E f0.3, OS] and t E (5.1, 7.51 in the adaptive analysis for the L-shaped domain.

As we know the exact solution of (52), we want to assess the accuracy and reliability of the proposed error estimator (35). The effectivity index for the new error estimator lleN I[ is plotted in Fig. 8 for At/T = 0.25 and 0.1, where T is the free vibration period. For the original Newmark method, the results by Zienkiewicz’s original error estimator Ileo 11and by the former error estimator 11 P 11are also given. The original estimator cannot give an accurate measure for local error, the former estimator converges asymptotically to the exact local error for most stages, and the new local error estimator (29), gives a more accurate and reliable estimation. For the linear acceleration method, the error estimator /fee// and Ileo”// is invalid.

time

Fig. 6. L-shaped domain: (a) temporal number of degrees of freedom used in the adaptive and the conventional analysis; (b) temporal error percentage achieved in the adaptive computation.

Li et al. E30) demunstra~e a number of numerical examptes that when included intcl an adaFt~ve tame-stealing procedure, the error estimator gives a simple, refiabfe and accurate procedure. The ~~rres~~nding result for the linear acceleration method is found in Fig. 9. We observe that results are only obtained for 11e”l], because the error estimates I[eoctI/ and /le@Ildo not hold for a linear acceleration scheme.

The shear beam model shown in Fig. 10 may approximately represent an off-shore tower subjected to a harmonic wave force near the surface [303, The adaptive time-stepping procedure is applied to this model to test its performance in a rest application. Without losing generality, we unfy show the results by the &near a~~e~erat~~n method. A relative focat error toolerance of 73= I.0 X 1O-49and three different sets of parameters ‘yeq y2 and Kg,are specified.

N.-E. W&erg et al., Error estimation and adaptivity

390 20 18 16 14 12 10 8 6

-

adaptive analysis

4

*-I-_

conventianal

2

FEM analysis

(NDOF=4554)

0 0.0

2.0

4.0

60

60

a.0

time (x

10-3 )

6 I

*-

-

adaptive analysis

21

- - - - -

conventional

FEM anaiysis

(NDOF=4554) o-2-f 49 -(I-6-

I

0.0

190

1

20

i

I

3,o

4.0

.

I

I

1

5.0

6.0

7.0

a0

time Fig. 7. L-shaped domain:

(a) horizontal

displacement

response at point B; (b) stress response a; at point C.

The estimated relative local error shown in Fig. 11 reveals that the local errors are reasonably controlled in the vicinity of the given error tolerance x = 0.8, y2 = 1.0, Ic, = 5. The time-step sizes used temporally are shown in Fig. 12. If q = 0.8 is chosen, the total number of time-steps for the computation in a time domain [0, 10 s) is 643, in which only seven different time-step sizes are used. From these examples, we see that the numerical solutions have been obtained in such a way that their local errors are effectively controlled in the vicinity of the given error tolerance by a proper time-step size adjustment. 7.2.3. Time ~~~e~~a~i~~,p-version, ad~p~~v~~ro~ed~r~ A number of numerical examples showing adaptive procedures found in [Ml.

for the p-version in time are

N.-E.

Wiberg et ul., Error ~~ti~atio~ and adaptivity

Fig. 8, Effectivity index of the locai error estimator for the original Newmark method for the S-DOF undamped freevibration. (a)bt/T=0.25,(b)dtlT=O.l.(t)-----, originat error estimator /P]t [29]; (2) * . . . . , former error estimator IIZ”ll [19, 221; (3) -, new error estimator jJPI[ [30].

8. Conclusions and adaptivity of spatial and time discretizatian in the semi-discrete finite element analysis for linear elastodynamic problems have been studied. We have analyzed an a priori estimate for the spatial discretization error in the total energy norm, and found that it is reascmable to extend the a posteriori error estimate of post-processing type to elast~d~am~c analysis. It is found that in most situations, the total energy of the error can be approximated Error

estimatiofi

392

Fig. 9. Effectivity S-DUF undamped

N.-E.

Wiberg et al., Error estimation and adaptivity

index of the presented free vibration modet.

40

local error estimator

t

39

//eNl] for the linear acceleration

method in the

M p(t) = lDo0 sin 0.7&

t* 38 qrn 37 t

K’ = i

_1-: 1 1

6854

M=49 3

m=7

2 I .i Fig. 10. A multi-DOF

exampie: a shear beam model.

by the strain energy of the error. Based on the extended error estimate and a mesh generator, we have presented a spatial mesh adaptation procedure which, by an automatic remeshing scheme, updates the spatial mesh when necessary to gain control of the spatial discretization error from time to time. Based on the numerical examples, we conclude that the error estimation and the adaptive procedure are capable of monitoring the moving of steep stress regions by updating the spatiaf mesh according to a given error tolerance, thus providing a reliable finite element solution efficiently+ An @-method for eigenvalue problems based on an error measure in the energy norm is

393

N.-E. Wiberg et al., Error estimation and adaptivity

ii = 1.0 x 10-4, & = 5 y1 = 0.8, y2 = 1.0

--

I

M

Al

4tu

slo

slo

Fig. 11. The relative local error achieved in the adaptive time-stepping. _._-

-1 0.04

M3!5

Ff =

1.0 x 10-J

y2 = l.O,&

= 5

Number of time steps Number of step changes CPU time

(1)

(2)

(3)

yl = 0.8 y1 = 0.5 y1 = 0.2 643

679

7

Fig. 12. Temporal time-step size predicted by the adaptive time-stepping original Newmark method.

1

2

2.03

2.52

procedure

756

2.36

for the shear beam with the

394

N.-E.

Wiberg et al., Error estimation and adaptivity

discussed. It shows that the exactness of the calculated eigenvalue depends heavily on the mesh and element approximation. A local error estimator IleNII- At3 for direct integration with the Newmark scheme, based on a postprocessed solution, has been presented in this paper. The numerical examples indicate that this error estimator is simple, reliable and accurate. An adaptive time-stepping procedure based on this error estimator is described. It can automatically adjust the time-step size so that the local error at each time-step is within a prescribed accuracy. In our test examples, it shows a good performance in error control, efficient computation and convenient implementation. New higher order schemes of hierarchical p-type for first order DE in time are presented which are A-stable of order 2, 4, 6, 8 and L-stable of order 1, 3, 5, 7. Selective p can be used. A residual type of error measure is used as a basis for an adaptive P-type scheme. A considerable improvement in terms of accuracy and effectiveness is found. Based upon the work presented in this paper, the authors believe that there is potential for the following further developments: to couple the spatial mesh adaptation procedure with the adaptive time-stepping procedure so that the discretization error can be controlled in space and time simultaneously; to develop global error estimates in the semi-discrete finite element analysis for elastodynamic problems so that the evolutional part of the discretization error can be further taken into account. Furthermore, as the size of elements in space can vary in different regions, the scheme should also permit different time-stepping (or different p-level) for different spatial regions. From this point of view, it is more attractive to develop adaptive procedures based on space-time finite elements. In this way, time and spatial variables can be treated in a consistent manner. It would therefore be interesting to combine the hierarchical time discretization with hierarchical finite element (space) formulations in order to obtain hierarchical formulations of for example, parabolic, hyperbolic and mixed problems.

References [l] I. BabuSka and W. Gui, Basic principles of feedback and adaptive approaches in the finite element method, Comput. Methods Appl. Mech. Engrg. 55 (1986) 27-42. [Z] A. Peano, A. Pasini, R. Riccioni and L. Sardella, Adaptive approximations in finite element analysis, Comput. & Structures 10 (1979) 333-342. [3] O.C. Zienkiewicz, 3. Gago and D.W. Kelly, The hierarchical concept in finite element analysis, Comput. & Structures 16 (1983) 53-65. [4] N.-E. Wiberg and P. Moller, Formulation and solution of hierarchical finite element equations, Internat. J. Numer. Methods Engrg. 26 (1988) 1213-1233. [5] I. Babuska, O.C. Zienkiewicz, J. Gago and E.R. de A. Oliveira, eds., Accuracy Estimates and Adaptive Refinements in Finite Element Computations (Wiley, New York, 1986). [5] B. Guo and I. BabuSka, The h-p version of finite element method, Part I: The basic approximation results. Part II: General results and applications, Comput. Mech. 1 (1986) 21-41, 203-226. [6] L.F. Zeng and N.-E. Wiberg, Adaptive procedures of h-p version for high accuracy finite element analysis in 2D linear elastic problems, Internat. J. Comput. & Structures 42 (6) (1992) 869-886. [7] B.A. Szabo, The p- and h-p versions of the finite element method in solid mechanics, Comput. Methods Appl. Mech. Engrg. 80 (1990) 185-195. [8] J.T. Oden, L. Demkowicz, W. Rachowicz and T.A. Westmann, Towards a universal h-p adaptive finite element strategy, Part 2. A posteriori error estimation, Comput. Methods Appl. Mech. Engrg. 77 (1989) 113-180.

N.-E.

[91 O.C. Zienkiewicz

Wiberg et al., Error estimation and adaptivity

395

and J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Internat. J. Numer. Methods Engrg. 24 (1987) 333-3.57. [101 O.C. Zienkiewicz and J.Z. Zhu, Superconvergent derivative recovery techniques and a posteriori error estimation in the finite element method, Inst. for Numer. Methods Engrg., University College Swansea, UK, CR/671/91, 1991. [Ill K. Eriksson and C. Johnson, An adaptive finite element method for linear problems, Math. Comp. 50 (1988) 361-383. u21 H. Jin and N.-E. Wiberg, Two dimensional mesh generation, adaptive remeshing and refinement, Internat. J. Numer. Methods Engrg. 29 (1990) 1501-1526. N.-E. Wiberg, Adaptive and hierarchical weighted residual and least-squares time integration of convection[I31 diffusion problems, Comm. Appl. Numer. Methods 4 (1988) 499-507. P41 0. Friberg, P. Moller, D. Makovicka and N.-E. Wiberg, An adaptive procedure for eigenvalue problems using the hierarchical finite element methods, Internat. J. Numer. Methods Engrg. 24 (1987) 319-335. P51 K.-J. Joo and Wilson, An adaptive finite element technique for structural dynamic analysis, Comput. & Structures 30 (1988) 1319-1339. Ml K.C. Park and P.G. Underwood, A variable-step central difference method for structural dynamics analysis Part I. Theoretical aspects, Comput. Methods Appl. Mech. Engrg. 22 (1980) 241-258. 1171P.G. Bergan and E. Mollestad, An automatic time-stepping algorithm for dynamic problems, Comput. Methods Appl. Mech. Engrg. 49 (1985) 299-318. WI P. Moller, High order integration methods for stiff systems, Report 91:4, Dept. of Structural Mechanics, Chalmers University of Technology, Sweden, 1991. [I91 O.C. Zienkiewicz and Y.M. Xie, A simple local error estimate and an adaptive time-stepping procedure for dynamic analysis, Earthquake Engrg. and Structural Dynam. 22 (1991) 871-887. PO1 L.F. Zeng, On adaptive finite element procedures for static and dynamic problems, PhD Thesis, Publ. 91:15, Department of Structural Mechanics, Chalmers University of Technology, Gdteborg, Sweden, 1991. WI N.-E. Wiberg and L.F. Zeng, Adaptive FE-procedures in elastomechanics, Final Proceedings, IV-ICCCBE, Tokyo, Japan, Pub1 91:13, Department of Structural Mechanics, Chalmers University of Technology, Goteborg, Sweden, 1991. PI L.F. Zeng, N.-E. Wiberg, X.D. Li and Y.M. Xie, A posteriori local error estimation and adaptive time-stepping for Newmark time integration in dynamic analysis, Earthquake Engrg. Structural Dynam. 21 (1992) 555-571. [231 L.F. Zeng and N.-E. Wiberg, Spatial mesh adaptation in semidiscrete finite element analysis of linear elastodynamic problems, Comput. Mech. 9 (1992) 315-332. P. 1241 Hansbo, Adaptivity and streamline diffusion procedures in the finite element method, Ph.D. thesis, Publ. 89:2, Department of Structural Mechanics, Chalmers University of Technology, Giiteborg, Sweden, 1989. P51 C.I. Bajer, Adaptive mesh in dynamic problems by the space-time approach, Comput. & Structures 33 (1989) 319-325. L.F. Zeng, N.-E. Wiberg and L. Bernspang, An adaptive finite element procedure for 2D dynamic transient WI analysis using direct integration, Internat. J. Numer. Methods Engrg. 34 (1992) 997-1014. 1271J.Z. Zhu and O.C. Zienkiewicz, Superconvergence recovery technique and a posteriori error estimators, Internat. J. Numer. Methods Engrg. 30 (1990) 1321-1340. PI J.T. Oden and J.N. Reddy, Note on an approximate method for computing consistent conjugate stress in elastic finite elements, Internat. J. Numer. Methods Engrg. 6 (1973) 55-61. [291 O.C. Zienkiewicz, W.L. Wood, N.W. Hine and R.L. Taylor, A unified set of single step algorithms. Part 1: General formulation and applications, Internat. J. Numer. Methods Engrg. 20 (1984) 1529-1552. [301 X.D. Li, L.F. Zeng and N.-E. Wiberg, A simple local error estimator and adaptive time-stepping procedure for direct integration methods in dynamic analysis, Department of Structural Mechanics, Chalmers University _ ot ‘I‘echnology, Goteborg, Sweden, 1991; Internat. J. Numer. Methods Engrg., to appear. [31] N.-E. Wiberg, L.F. Zeng and H. Tagnfors, Numerical methods in structural dynamics, adaptive FEprocedures, Slaveb. 39, No. 9 (1989) 697-704, Bratislava.