Analysis of shell structures under transient loading using adaptivity in time and space

Computers and Structures 79 (2001) 2117±2131

www.elsevier.com/locate/compstruc

Analysis of shell structures under transient loading using adaptivity in time and space Jens Neumann, Karl Schweizerhof * Institute for Mechanics, University of Karlsruhe, D-76128 Karlsruhe, Germany

Abstract The adaptive analysis of structures under transient loading leads to the question, which time integration scheme, ®nite elements or ®nite di€erences, is most favorably combined with an adaptive spatial FE discretization. In order to judge this, the properties of di€erent discontinuous Galerkin (DG) and the standard Newmark method are investigated ®rst, also concerning eciency. In particular, the damping and dispersion e€ects are discussed in detail for various types of problems. It must be noted that the type of problem has to be carefully checked in order to apply the most appropriate and ecient time integration scheme. It is shown that e.g. for the wave propagation problems the DG method with linear approximations (DG P1-P1) has to be favored when adaptivity in space is applied. Finally, an adaptive time step modi®cation scheme is presented and applied to various problems. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Transient analysis; Solid mechanics; Error estimation; Adaptive mesh modi®cation; FE-analysis

1. Introduction Complex solid mechanics problems with external time dependent loading and other excitation often require a very high solution e€ort, and numerical schemes have to be used for spatial discretization and time integration. As the loading of the di€erent parts of the structures varies not only within the structure but also with time, adaptive solution procedures seem to be very attractive to achieve a high eciency of the solution process with controlled errors. The standard schemes for error estimation in FE meshes known from statics cannot be applied without considering the kinetic and damping e€ects as well as the time integration e€ects. It is particular important to keep in mind that with kinetic e€ects also many physical phenomena appear which have to be well represented by the solution methods.

*

Corresponding author. E-mail address: [email protected] (K. Schweizerhof).

The contribution is therefore directed towards the capturing of kinetic e€ects, respectively the interaction of the numerical solution schemes with various details of the solution. First, discontinuous Galerkin (DG) time integration methods are compared with standard Newmark type methods for di€erent types of problems using stationary and adaptively modi®ed meshes. Then, the focus is on dispersion e€ects, in particular introduced by the spatial discretization and their representation with the two integration methods discussed before. This is followed by an investigation, if and how mesh adaptation changes with the e€ects observed with stationary meshes. Finally, a simple and ecient time step modi®cation method is presented and discussed using a numerical example. 2. Time integration ± discontinuous Galerkin vs. Newmark type methods In order to take also kinetic aspects into account Cauchy's ®rst equation of motion for solids has to be used to describe the behavior of solid type structures under general loading (Fig. 1):

0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 1 ) 0 0 0 5 9 - 1

2118

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

r ˆ C:. The test space Y and the space of the trial functions X are de®ned as: v 2 Y ; Y ˆ f…H 1 †n ; v ˆ 0 on C1 g, u 2 X ; X ˆ f…H 1 †n ; u…x† ˆ u1 on C1 g. Thus, the ®rst step of semidiscretization leads to a second order ordinary di€erential equation (ODE) system. With the assumption of modal damping D ˆ aM ‡ bK we obtain: _ ‡ Ku…t† ˆ f…t†: M u…t† ‡ Du…t†

…4†

In general this ODE system (4) has to be integrated numerically. 3. Numerical time integration

Fig. 1. Solid under general loading.

div r ‡ qd ˆ qu;

…1†

with r: linear Cauchy-stress tensor and d: acceleration ®eld (body force),  u…x; 0† ˆ u0 ; initial displacements; _ 0† ˆ u_ 0 ; initial velocities u…x; and u…x; t† ˆ u1 …t†; 8x 2 C1 ; r…x; t† n ˆ p…t†; 8x 2 C2 :

boundary conditions

B

Applying the Green±Gauss theorem for the ®rst term of Eq. (2): Z Z Z vdiv r dV ˆ vrn dA rvr dV ; …3† C2

un‡1 ˆ un ‡ Dtu_ n ‡ Dt2 ……0:5 u_ n‡1 ˆ u_ n ‡ Dt……1



Initial and boundary conditions must be given in a standard form, such that the problem is well de®ned. The standard procedure for the solution in solid mechanics ± only brie¯y described in the following ± involves a semidiscretization, thus the discretization in time and space is performed separately (see e.g. Refs. [8,27]) e.g. both with FE methods (see Ref. [9]). Other methods, though promising (see Ref. [12]) have not yet been successful in practice. For the spatial discretization the standard Bubnov± Galerkin formulation is chosen: Z v…div r ‡ qd qu† dV ˆ 0: …2†

B

Mainly either ®nite di€erences (FD) or ®nite elements in time (FET) are used for the time integration (see Refs. [3,6,7,9,11,17,26]). Both types of methods have in common that the total integration interval t 2 ‰0; T Š is divided into `local' intervals In ˆ ‰tn ; tn‡1 Š with a time step size Dt ˆ tn‡1 tn . Thus a piecewise integration is performed. In the following, the major di€erences of both approaches are investigated, in particular with respect to spatial adaptivity. The most popular FD scheme is the basic approach of Newmark [19]:

B

leads to the well known weak form of Eq. (1). The following investigations are restricted to linear behavior with the linear strain tensor : ˆ 1=2…ru ‡ ruT † and linear material laws, e.g. the St.-Venant Kirchho€ model

b† un ‡ b un‡1 †;

c† un ‡ cu_ n‡1 †:

…5† …6†

The parameters b and c control the accuracy, arti®cial damping and the stability of the method. Standard FET have been proposed already since long e.g. by Oden [21], who used the Hamiltonian and the variation as the starting point of his derivations. Zienkiewicz [26] showed in a paper the equivalence between the Newmark scheme and a distinct FET method. All these schemes suggested were based on continuity assumptions at the beginning resp. end of the local time interval. More recent developments, see Refs. [1,9,13] as well as Refs. [16,23] allow discontinuities of the state variables at the internal boundaries with the penalty of an enlargement of the equation systems. In general, various ways to develop the ®nal algebraic system of equations can be chosen. If the ODE (4) is directly taken as starting point, a so-called two-®eld formulation with displacements and velocities as variables seems to have some advantages: u_ ˆ v;

…7†

M_v ‡ Dv ‡ Ku ˆ f:

…8†

A Petrov±Galerkin method is used for the solution of this coupled system; both equations are weighted with di€erent test functions wu and wv and then integration is

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

0

1 K 2

1 K 2

1 Dt K 3 n

1 K 2

1 Dt K 6 n

B B 1K B 2 B B1 B 3Dtn K @

1 Dt K 6 n

1 Dt K 6 n

1 Dt K 3 n

0

1 M 2

Kun B B 0 B ˆB B F1 ‡ Mv @ n

Fig. 2. Time discretization involving jumps in displacements.

with: performed over the time interval Dt. In addition, the jumps in displacements Dun ˆ u‡ un and velocities n Dvn ˆ v‡ v at the beginning of the time interval are n n introduced (see Fig. 2) and the following (variational) form is obtained: Z

Z ‡ wv …M_v ‡ Dv ‡ Ku†ds ‡ wu K…u_ v†ds ‡ w‡ v
Mvn s s Z ‡ ‡ ‡ w‡
Ku ˆ wv f…t†ds ‡ w‡ n u v
Mvn ‡ wu
Kun ; s

8fwu ; wv g 2 fH 1 …0; Dt†  H 1 …0; Dt†g:

…9†

Due to the jump terms other test functions with less restrictions can be used, as continuity is not required. Such schemes are known as DG methods (see Ref. [14]). The test and trial functions for displacements Su and velocities Sv can be chosen as any order polynomials. Usually, low order schemes ± linear or quadratic ± are preferred [11,16,23] due to several aspects discussed in the following. For example if linear functions are chosen for both displacements and velocities, the so-called DG P1-P1 method is obtained: Su ˆ fuh ; wu 2

N [

1

…P …In ††

neq

g;

nˆ1

Sv ˆ fvh ; wv 2

N [

…P1 …In ††neq g:

nˆ1

P1 is the space of ®rst order polynomials. uh ˆ …1

with :

s†u‡ n ‡ sun‡1 ;

s ˆ h=Dtn ;

vh ˆ …1

Dtn ˆ tn‡1

s†v‡ n ‡ svn‡1 ;

tn ;

h 2 …0; Dt†:

The resulting equation system can be written as (without any further transformation):

F2 R F1 ˆ s …1

10

1 Dt K 6 n

u‡ n

1

C CB CB u C CB n‡1 C C CB CB ‡ C 1 1 B C C v M ‡ Dt D n n 2 6 A A@ 1 Dt K 3 n

‡ 13Dtn D

1 M 2

1

2119

‡ 16Dtn D

1 M 2

‡ 13Dtn D

C C C C C A

vn‡1

…10†

s†f…t†Dtn ds;

F2 ˆ

R s

sf…t†Dtn ds:

3.1. Comparing discontinuous Galerkin and Newmark type algorithms After the description of the general schemes of both time integrators our focus is on the speci®c features, the advantages and disadvantages. It is well known that essential properties concerning the quality of time integrators can be found through investigations of the ampli®cation matrix A (A 2 …R†22 ) of a single degree of freedom (SDOF) system. With the eigenvalues k1 , k2 of A general statements concerning the spectral stability, algorithmic damping and the phase error (time dispersion) can be developed. For an undamped linear SDOFsystem with parameters m (mass) and k (sti€ness) the following quantities are obtained: Phase: X ˆ x Dt, Spectral radius: q ˆ max jki j;

i ˆ 1 . . . 2,

Natural frequencies, FET-FD: xIh ˆ

arctan Im…k† Re…k† Dt

;

xII h ˆ

p 2

‡ arctan Re…k† Im…k† Dt

:

…11†

x is the exact natural frequency, xh is the approximated frequency from the time integration scheme. The ®rst equation for xIh (12) is valid in the ®rst quadrant of the complex numbers, while the second equation for xII h stands for the second quadrant. The formulas for the other quadrants can be found analogously. There are two bifurcation points for the DG P2 formulation at X ˆ 3:1 and X ˆ 10:8 (visible in Fig. 3). Beyond X ˆ 10:8 all roots are real, thus the numerical phase remains constant there. It is known that both the DG and the standard Newmark scheme are `A' stable, however there are major di€erences concerning the spectral radii (see Fig. 3) and the algorithmic damping. The spectral radius of the Newmark scheme is independent of the phase, whereas the DG P1-P1 scheme shows a strongly decaying radius for larger phase resp. higher frequencies. The values of the DG P2 scheme, however, do not show a

2120

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 3. Spectral radii q vs. phase for Newmark and DG time integration.

clear dependency, in particular, for the medium range of the phase. Remarkably, the spectral radius remains 1 for the larger phase, thus higher frequency range. This has consequences for the algorithmic damping n and the phase error ep , which are de®ned as: nˆ

ln q ; Dtxh

ep ˆ

X Xh

1;

with : Xh ˆ Dtxh :

…12†

Their properties for a wide range of frequencies are shown in Fig. 4a resp. b. As already visible from Fig. 3 there is no dependency of the spectral radius from the frequency for the Newmark scheme thus no algorithmic damping over the full range of frequencies. There is a varying dependency of the damping in the low frequency range and no damping of the high frequencies for the DG P2 scheme. A somehow better behavior concerning damping can be found for the DG P1-P1 method. All high frequencies will be damped out. Thus, the phase error shown in Fig. 4b, which causes a dispersion in the time range, will also be damped out, though the phase error in time is growing steadily. It is also remarkable that the phase error of the Newmark scheme remains below the error of the DG P1-P1 method in the range X=2p > 1:2, however, it is not damped out. The phase

error of the DG P2 scheme remains well below the error of the Newmark scheme for the full frequency range and similar to the Newmark scheme there is no damping of this error in the high frequency range. Comparing the accuracy of the DG P1-P1 and the DG P2 method by Taylor series expansion we get for the remainder: DG P1-P1 : s ˆ 361 O…Dt3 †; see Ref: ‰11Š:

DG P2 : s ˆ 721 O…Dt3 †;

Thus, there is only little di€erence in the accuracy. In order to get the arti®cial dissipation in the upper modes also for the DG P2 scheme least-square operators or socalled shock capturing operators [10] could be used. Beyond the above mentioned damping properties there are other aspects, which favor the DG P1-P1 formulation instead of the DG P2 method. First of all, this is the symmetric block-matrix form for the DG P1-P1 scheme with identical main-diagonals, which is obtained for both DG P1-P1 formulations [16,23] discussed in detail in the following. Taking the block structure into account Wiberg/Li obtain the following system of four rows of block matrices after a FE discretization:

Fig. 4. Algorithmic damping resp. phase error vs. phase X=2p for Newmark and DG integration schemes.

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

0

1 I O DtI 6 1 BO I DtI 2 B @O O M O O 12DtD ‡ 13Dt2 K 0 1 un B un C C ˆB @ F A 1  F2

10

1 u‡ DtI n 6 1 CB un‡1 DtI 2 B C 2 M ‡ 16DtD A@ v‡ n 3 vn‡1 M

1 Uˆ

C C A

…13†

with: M ˆ M ‡ 12DtD ‡ 16Dt2 K;

M 2 Rneqneq

It is obvious that the third and the fourth row are fully decoupled from the ®rst two rows, thus only a system of 2  neq equations remains. Furthermore, this large nonsymmetric equation system can be favorably solved by Gauss±Seidel block iteration. For an ecient iterative solution the spectral radius of the iteration matrix should be small, which may not be the case for some problems in general. The iterative solution for the velocities is then written as (see Ref. [16]): Iteration i:

…i†

M vn‡1 ˆ F2

…i 1†

…23M ‡ 16DtD†vn‡1

…14†

…12DtD ‡ 13Dt2 K†v‡…i† n

…15†

Convergence check with:  ‡…i†   ‡…i 1†  vn vn …i 1† v …i† v n‡1 n‡1  ‡…i†  diff ˆ vn v …i†

…16†

if diff 6 tol then ®nish; else continue with next iteration cycle. The advantage of the iterative block solution in the manner shown above is that the LDLT decomposition of the matrix M has to be performed only once. Thus, only forward±backward substitutions are required for the solution in each iteration cycle. After convergence of the solution for the velocities the displacements can be computed without much e€ort from the two remaining sets of equations, row one and two: ‡ 1 u‡ n ˆ un ‡ 6Dt…vn

X

…Bu†T CBu dV ;

vn‡1 †:

An alternative to the scheme of Wiberg/Li is suggested by Ruge [23] based on the variation of the scalar integral I, where H ˆ T ‡ U , kinetic plus strain energy, is the Hamiltonian:

T ˆ

1 2

Z X

pT q 1 p dV ;

p ˆ qv:

The integral I over a arbitrary time interval Dt ˆ t1 t0 R is: I ˆ Dt H dt. The standard requirement concerning I (referring to Hamilton's principle) is, that it remains stationary, once the work of the external forces is included. Thus the variation becomes: Z Z dI1 ˆ dI ‡ duT f dX dt ˆ 0: X

By partial integration with respect to time we obtain:   Z Z oH dX dt dpT u_ dI1 ˆ op Dt X   Z Z oU f dX dt ˆ 0: …17† duT p_ ‡ ou Dt X Usually this formulation leads to a conditionally stable integrator with an order of accuracy O…Dt†, if u, v 2 P 1 …t†; t 2 ‰tn ; tn‡1 Š. To overcome the conditional stability, Ruge [23], following Aharoni/Bar-Yoseph [1], proposed to add the work of the state variables at the beginning of the time step to the integral, thus introducing them as independent variables. Then the varia‡ tion du‡ n and dvn are now also not zero leading to the following form: Z Z ‡ ‡ ‰u u ŠdV q du‡ vn ŠdX dI  ˆ dI1 ‡ q dv‡ n n n n ‰vn X

X

ˆ 0: With linear test functions u…t† and v…t† ± a P1-P1 scheme ± the following non-symmetric set of equations is developed: ! ! 2 1 Dt M ‡ 23D ‡ 185 DtK D Dt9 K un‡1 Dt 6 1 M Dt

n‡1

un‡1 ˆ un ‡ 12Dt…v‡ n ‡ vn‡1 †;

Z

Dt

F1 ˆ 53…F1 ‡ Mvn † 13F2 23DtKun ; F1 ˆ F1 ‡ F2 ‡ Mvn DtKun :

M v‡…i† ˆ F1 n

1 2

2121

ˆ

2

‡ 23D ‡ 187 DtK 1 M ‡ 23D Dt 1 M ‡ 13D Dt

M ‡ Dt6 D Dt18 K !  2 DtK M un 9 1 DtK 9

O

vn

‡

vn‡1   f1 f2

Written in short form as : Az1 ˆ Bz0 ‡ r:

…18† …19†

This system can be transformed to a similar structure with identical block diagonals as for the rows 3 and 4 of the Wiberg±Li scheme. With speci®cally chosen hypermatrices the following system of equations is found: !   1 1 0 1 T1 ˆ ; T2 ˆ 2 ; 1 Dt3 1 1 1 Dt with Eq: …19† and z1 ˆ T2 y1 T1 AT2 y1 ˆ T1 Bz0 ‡ T1 r; ! A y1 ˆ B z0 ‡ r

…20† …21†

2122

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 5. Spectral radius qit vs. sti€ness k for the DG P1-P1 schemes proposed by Wiberg/Li [16] and Ruge [23] combined with relaxation; SDOF system: m ˆ 0:5 ˆ constant, Dt ˆ 0:1 ˆ constant.



with A ˆ

M Dt

 ‡ 23D ‡ 185 DtK 2M ‡ 56D ‡ 29DtK Dt : D ‡ 29DtK A11 3

…22†

Then the same iteration procedure with a block Gauss± Seidel scheme as described for the Wiberg/Li equation system can be used. However, the major di€erence between both equation systems is that Ruge's scheme contains the variables un‡1 and vn‡1 , whereas Wiberg/Li solve for the velocities v‡ n and vn‡1 directly. Thus, a direct comparison is dicult and the judgment is preferably based on the eciency of the solution. The latter is mainly in¯uenced by the condition number of the system matrices, which can be decomposed as follows: A resp: A ˆ

L‡S

U;

…23†

S . . . diagonal matrix of A (13) resp. A (22); L, U . . . lower resp. upper triangular matrix of matrix A resp. A . A comparison shows that the diagonal of matrix A (Eq. (13)) of the Wiberg scheme is far more dominant than the diagonal of matrix A (Eq. (22)). In Fig. 5 the spectral radius qit of the respective iteration matrix Tit ˆ …S L† 1 U is presented for both integration schemes applied to a SDOF with varying sti€ness. It is obvious that the spectral radius is much smaller for Wiberg/Li's scheme resulting in a better convergence. Furthermore, an attempt was made to improve the spectral radius by a combination with a relaxation scheme. Then the modi®ed iteration matrix Tit becomes:

Tit ˆ …S

gL† 1 ‰…1

g†S ‡ gUŠ:

…24†

In the example (SDOF) the relaxation parameter g was varied between 0.95, 1.00 and 1.1. For the chosen system parameters a considerable improvement in the spectral radius was found for g ˆ 1:1. However, in real systems with a widespread band of frequencies no visible advantage using such a relaxation scheme could be achieved. Finally, the wave propagation in a slab (see Fig. 6) is considered to compare the performance of the solution schemes without relaxation using a constant time step size and a constant mesh. From Fig. 7 it can be concluded that the e€ort required for Ruge's scheme is considerably higher at least in the beginning of the process. Similar observations have been made for other problems. 4. Comparing the performance in time of DG and Newmark time integration for wave propagation and vibration problems In order to show the di€erence between a DG and the standard Newmark method two completely di€erent types of problems are investigated; an example of an instationary problem and one for stationary motion. A classical example for an instationary problem is the wave propagation in solids, while the vibration of a cantilever plate represents a stationary vibration prob-

Fig. 6. Wave propagation problem; system and load function vs. time; material parameter: E ˆ 1  104 , q ˆ 10, t ˆ 1, v ˆ 0:0.

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 7. Time history of number of iterations Nit for DG P1-P1 methods; example ± wave propagation.

lem. For both cases di€erent parts of the frequency spectrum are dominantly excited. In the wave propagation problem (see Fig. 8) the response is in the complete frequency range, whereas in the stationary problem predominantly the response is in the lower frequencies. In the left diagram in Fig. 8 the time history for the total energy norm of the solution jjujj is presented. The dissipative property of the DG P1P1 is obvious. About 25% of the energy is dissipated after 80 time units, whereas the energy is perfectly conserved for the Newmark scheme. The comparison of the velocity component vx shows the positive aspects of the DG P1-P1 scheme on the other hand. As a result of the algorithmic damping, the higher (less accurate) frequencies are ®ltered out and a smooth response is achieved. It is also obvious that both schemes cannot

2123

maintain the correct velocity amplitudes. Although the Newmark scheme conserves the energy, the velocity amplitude shows a strong variation due to the inherent spatial dispersion of the FEM-mesh. The conclusion is rather clear: Both the DG P1-P1 and Newmark schemes are not applicable for long time integration in such a case. However, for short duration problems the DG P1P1 scheme performs rather well. The second case, a stationary motion of a clamped cantilever plate with a line load at the free end (see Fig. 9) is investigated next. In the diagrams in Fig. 10 the total energy norm of the solution and the vertical velocity vz at the free end are presented. As expected there is hardly any di€erence between the results obtained with both schemes. Obviously higher frequencies are not excited and there is no accuracy advantage for any method. However, the solution e€ort needed for the Newmark scheme is far below the e€ort required for any DG scheme for such a problem and it should be preferred then.

5. Dispersion e€ects depending on the spatial FE discretization 5.1. General The time integration schemes discussed above showed some phase error (Eq. (12)), thus some dispersion changing in time was already found. A similar e€ect is introduced by the spatial discretization. This arti®cial

Fig. 8. Wave propagation problem ± time history plot of (a) total energy norm jjujjtotal;X and (b) velocity vx at x ˆ 140; comparison of Newmark and DG P1-P1 schemes.

Fig. 9. Clamped plate with line load; example for stationary vibration.

2124

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 10. Clamped plate, stationary motion problem ± time history plot of (a) total energy norm jjujjtotal;X and (b) velocity vz at the end of plate; comparison of Newmark and DG P1-P1 time integration.

Fig. 11. Propagation velocity vs. inverse of wavelength. (a) Flexural wave in plates after Mindlin [18]. (b) Longitudinal waves in 1D and 3D continuum after Gra€ [5], c ˆ wave number.

dispersion, however, has to be separated from any dispersion which is inherent in the problem described by the corresponding di€erential equation. Such phenomena can be investigated rather well in wave propagation problems. There, every wavelength ki has its own propagation velocity c…ki †. In Fig. 11a the inherent dispersion for the wave propagation of a ¯exural mode in plates is shown depending on the kinematic model chosen [18]. Every model (Mindlin plate, Kirchho€ plate and 3D continuum) has its own dispersion characteristics which are discussed brie¯y in the following with the de®nitions: h ± thickness of the plate, k ± wavelength of ¯exural wave, c ± propagation velocity of for ¯exural wave, cs ± propap gation velocity of the shear wave; cs ˆ E=2…1 ‡ m†q. From Fig. 11a it is obvious that the Kirchho€ model is able to represent the bending waves only for large wavelengths resp. low wave numbers, and the kinematic model of Mindlin corresponds very well to the exact solution of the general 3D model for the full range of wave numbers. Comparing the longitudinal resp. membrane waves for the 3D continuum and 1D continuum e.g. a bar with circular cross-section, a signi®cant dispersion is found for the exact of the 3D case psolution  with the de®nitions: c0 ˆ E=q ± propagation velocity

of 1D continuum, c ± propagation velocity of longitudinal wave; k ± wavelength of longitudinal wave. However, it is also clear that the simplest reduced continuum model, the 1D case is able to represent the exact solution of the 3D case …c ˆ c…k†† only for large wavelengths. For the longitudinal waves in Fig. 11b the exact dispersion in 3D converges to a constant value;  ˆ 0:57. here: limk!0 c…k† Now the non-dispersive characteristics of the 1D wave propagation problem are used to investigate the arti®cial dispersion due to the spatial discretization. The well known linear homogeneous di€erential equation for the 1D wave propagation problem in a bar with a constant cross-section is ± assuming for completeness some boundary and initial conditions: q

o2 u ot2

E

o2 u ˆ 0: ox2

…25†

For the solution the standard assumption is taken: u…x; t† ˆ A exp ix…x=c t†. Without loss of generality the following can be used as test and trial functions for discretization in space: test space : v…x† 2 H01 …X†;

trial space : u…x† 2 H 1 …X†:

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

ch ˆ

Fig. 12. Discretized bar; unknown quantities at nodes.

k h xh l : ˆ a Th

2125

…28†

Thus, it can be concluded that even for the most simple problem without dispersion the corresponding FE discretization introduces arti®cial dispersion. It is obvious that a strict distinction between arti®cial and physical dispersion is not possible for 2D and 3D continua. 5.2. Uniform discretization in space

With the standard Bubnov±Galerkin method, weighting with the test functions and using Gauss±Green's theorem the ODE for node j (see Fig. 12) in a FE mesh with uniform element length l is obtained: c20 … uj‡1 ‡ 2uj l2 j ± index for nodes:

 uj‡1 ‡ 4uj ‡ uj

1

‡6

uj 1 † ˆ 0; …26†

Using the same assumption as for the exact solution: uj ˆ C exp i…xh t aj†, with: C ± constant to ful®ll the initial conditions, xh ± circular frequency, k ± wavelength, a ± normalized wave number; a ˆ 2p l=k, the following relation between the normalized wave number a and the circular frequency xh for the discretized solution is found, see Ref. [20]: xh ˆ

c0 l

r 1 cos a : 6 2 ‡ cos a

…27†

Contrary to the analytical exact solution the approximate solution has a so-called cut-o€ frequency xco , the maximum frequency which can be represented in the discretized system (see Fig. 13a). Thus, using a ®xed l for an element, one can only ®nd a maximum frequency xco corresponding to a constant ratio l=k. As a consequence of this non-linear wavelength±frequency relation a nonconstant propagation velocity ch , which also depends on the wavelength, is obtained:

In order to show the combined e€ect of space and time discretization, a wave propagation example, the 1D continuum problem as shown in Fig. 6, is investigated. The mesh is kept uniform and unmodi®ed during the analysis. Three types of time integration are compared: First, the fully exact solution in time and space is given by an analytical solution. Second, the structure is discretized in space by FE and followed by modal decomposition. Then each modal response is integrated separately and summed up exactly as modal solution. As third action the FE discretization in space is combined with either FET (DG P1-P1) resp. FD (Newmark) in time. In the following diagrams the stress state r…x† of the whole space domain is plotted at several times. As expected (see Fig. 14) the modal solution is somehow in front of the exact solution because of the arti®cial dispersion inherent in the discretized continuum. The short waves have a higher propagation velocity than the long waves causing the waviness in front of the sinusoidal half-wave. With progressing time this arti®cial dispersion increases. The Newmark solution shows some arti®cial waves traveling behind the exact and the modal solution. This is mainly due to the phase error (see Fig. 4b) which accumulates in time. However, it may be noted that taking a much smaller time step, a similar solution as for the modal method could be obtained. For the DG P1-P1 method (see Fig. 15) the modal solution is approximated rather well mainly due to the order of accuracy O…Dt3 † of the integration

Fig. 13. Circular frequency resp. propagation velocity vs. ratio of element length/wavelength of discretized 1D continuum. (a) Circular frequency with cut-o€. (b) Propagation velocity with cut-o€.

2126

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 14. 1D wave propagation ± stress state of complete space domain at di€erent time states; uniform mesh; behaviour of Newmark scheme compared to analytical (d-almb) and modal solution. (a) t ˆ 1:0, (b) t ˆ 4:0, (c) t ˆ 9:0 and (d) t ˆ 17:2.

Fig. 15. 1D wave propagation ± stress state of complete space domain at di€erent times; uniform mesh; behavior of DG P1-P1 scheme compared to analytical (d-almb) and modal solution. (a) t ˆ 1:0, (b) t ˆ 4:0, (c) t ˆ 9:0 and (d) t ˆ 17:2.

scheme. The e€ect of the phase error seems to be small for this type of excitation which is also due to the arti®cial damping of the higher frequencies. The latter also

leads to some decay of the major amplitude and to reduced amplitudes of the small waves running in front of the sinusoidal large wave.

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

2127

Fig. 16. 1D wave propagation ± stress state of complete space domain; comparing DG P1-P1 and Newmark time integration with exact solution; adaptive discretization in space. (a) t ˆ 1:0, (b) t ˆ 4:0, (c) t ˆ 9:0 and (d) t ˆ 17:2.

5.3. Adaptive discretization in space If instationary meshes are used, e.g. adaptive mesh generation takes place, then the situation changes. In the following for the estimation of the spatial discretization error the well-known estimator of Zienkiewicz=Zhu [29] for the strain energy norm of error and an additional indicator for the kinetic energy part, proposed by Riccius and Schweizerhof [22] are applied. In order to judge the e€ect of the time integration schemes the time step is kept constant. Thus, the exact analytical solution and fully numerical solutions are compared. Due to the graded meshes, generated in the adaptation process, the maximum amplitude in the response rx presented in Fig. 16 at various time states is considerably reduced for both numerical schemes. Though the arti®cial dispersion cannot be controlled perfectly, very good results can be achieved with the DG P1-P1 method. Only a small increase of the width of the wave and no longer additional small waves are obtained. Furthermore, there is only a very small phase error visible after a while, e.g. at time t ˆ 17:2. However, the solution errors with the Newmark method increase considerably with time. In particular, the arti®cial wave following the physical wavefront grows. The reason is the mesh re®nement in space also for the region of these smaller waves. An amplitude decay is visible for both methods. For the energy conserving Newmark scheme this is caused by the energy transfer from the long wave to the smaller

arti®cial waves. The latter is also the reason for the energy decay ± by far less pronounced ± in the DG P1-P1 method, as it is not possible to damp out the higher modes without a€ecting the lower modes. Thus, the arti®cial damping of the DG P1-P1 scheme will damp out all modes in time as t ! 1. Typical meshes created in the adaptive process are shown in Fig. 17. Spurious re¯ections: Besides the discussed e€ect of arti®cial dispersion the adaptive meshing could also create so-called spurious re¯ections. Bazant/Celep [2] investigated this phenomena on stationary graded meshes with the conclusions that only short waves will be re¯ected at the element interface in the case l2 P 2l1 (see Fig. 18) and l1 P k=10, if a graded mesh is created. However, such an e€ect was not observed during all investigations with instationary (adaptive) meshes so far. 6. Adaptive time stepping Though in general it is possible to integrate the complete time dependent problem with e.g. time-space ®nite elements with di€erent time step size for each spatial element, see Ref. [9] or at least parts of structures with di€erent time step size and/or di€erent time integration schemes (see e.g. Refs. [8,12]) it is favorable from a practical point of view to obtain the total information of the status of a structure at a given time. Thus integration schemes with a uniform time step size Dt for the

2128

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 17. Adaptive meshes at various times created in combination with two di€erent time integration schemes; colored: distribution of the total energy. (a) meshes with Newmark method, (b) meshes with DG P1-P1 method.

Fig. 18. Spurious re¯ections in graded meshes.

complete structure are discussed and our focus is on their properties within a typical interval In ˆ ‰tn ; tn‡1 Š with Dt ˆ tn‡1 tn . The error within an interval can be viewed as local error and its accumulation over the complete process time as global error. The latter is of major interest for the ®nal result and the question is, how to obtain a proper information about this error in order to adjust the time step size. Currently, some global schemes [4,24] based on duality considerations are developed. However, though they give some hints about the global error, their use for real-size FE problems seems to be limited. Thus, our focus is on the local error estimates and the use of such an estimate for a global estimate and time step adaptation. The error in displacements and velocities at time t ˆ tn‡1 is given as en‡1 ˆ uex;n‡1

uh;n‡1 ;

e_ n‡1 ˆ u_ ex;n‡1

u_ h;n‡1 :

…29†

For the DG P1-P1 method the velocity error is given by the separately interpolated velocity-®eld v. In order to couple the time integration error with the spatial dis-

cretization error the total energy-norm of error can be computed (see Ref. [22]): Z ke…Dt†ktotal;X;tn‡1 ˆ …Len‡1 †T
CLen‡1 dX B

‡

Z B

q < e_ n‡1 ; e_ n‡1 > dX

1=2 : …30†

For the Newmark-scheme a Taylor series expansion can be performed to obtain some information about the local error within a time step, assuming the quantities in the previous time steps are computed exactly: Dt2 Dt3 * Dt4 <  un ‡ un ‡ un ‡ O…Dt5 †; 2 6 24 Dt2 * Dt3 < un ‡ un ‡ O…Dt4 †:  u_ n ‡ Dt un ‡ 2 6 …31†

uex;n‡1  un ‡ Dtu_ n ‡ u_ ex;n‡1

To eliminate terms of order 3 and 4 in Eq. (31) the central di€erence scheme;

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131 *

With a required energy norm of the error ev;T at the end of the process the constant c can be computed:

un ˆ un‡12Dtun 1 ‡ O…Dt2 †; < un

ˆ

 un ‡ un un‡1 2 Dt2

1

2

‡ O…Dt †;

is used twice, similar to Li [15]. Then the displacements and velocities depend only on displacements, velocities and accelerations of the last two time steps plus higher order remainders: 2

uex;n‡1  un ‡ Dtu_ n ‡ Dt  1 u u_ ex;n‡1  u_ n ‡ Dt 12 n

 1

1  u 24 n 1

‡

2  u 3 n

‡

‡

5  u 12 n

5  u 12 n‡1

1  u 8 n‡1



5

‡ O…Dt †; ‡  4 ‡ O…Dt †:

With the Newmark interpolation for the displacements uh;n‡1 and velocities u_ h;n‡1 , see Eqs. (5) and (6) the local error estimate for the Newmark scheme can be given: en‡1 ˆ

e_ n‡1 ˆ

Dt2 …un 24

Dt …un 12

1

‡ …2

24b†un ‡ …24b

3†un‡1 † ‡ O…Dt5 †; …32†

1

‡ …4

12c†un ‡ …12c

5†un‡1 † ‡ O…Dt4 †: …33†

For the DG P1-P1 scheme a very simple local error indicator can be constructed (see Refs. [16,23]) based on the discontinuity of displacements and velocities: en‡1 ˆ u‡ n‡1

un‡1 ;

2129

e_ n‡1 ˆ v‡ n‡1

vn‡1 :

…34†

However, the latter cannot be used for the global error estimation. As a short remark for the DG P2 method, we have to note that there is no displacement discontinuity in time, because u‡ n‡1 ˆ un‡1 . That is another aspect, why the DG P1-P1 can be considered as superior concerning eciency when compared to DG P2. The global error of the DG schemes discussed is certainly also dependent on their properties concerning the high frequency damping, thus global error estimates based on the rather simple local error estimates are dicult and an error estimation function as proposed in Ref. [24] seems appropriate for the DG methods. For the Newmark scheme a rather simple formula for the global time integration error ke…Dt†ktotal;X;T is based on an extrapolation proposed by Zienkiewicz/Xie [28] and Warburton [25]: ke…Dt†ktotal;X;T ˆ NT ke…Dt†ktotal;X;tn‡1 ;

cˆ ˆ

ev;T ke…Dt†ktotal;X;T

!1=2

ev;T Dtold T ke…Dt†ktotal;X;tn‡1

!1=2 :

…37†

For the limitation of the estimated error at the end of the process after Eq. (35) an upper limit is given for accuracy and a lower limit for eciency reasons: emin;T 6 ke…Dt†ktotal;X;T 6 emax;T . The time step size is modi®ed, if any of the tolerances is violated. Then the question remains, how often the time step size has to be adjusted. As the matrices remain constant with constant time step size, then it is advantageous to keep the same time step size, if it is smaller than the required one, as the factorization of the matrices can be avoided then. It is suggested, however, that a time step size adjustment is also performed after each mesh modi®cation. The same strategy, global error estimation and time step control, was applied as a simple test with the DG P1-P1 scheme. As numerical example the clamped plate system (see Fig. 9) is taken to show the e€ects of both error estimators resp. indicators (Eqs. (32)±(34)). For both the Newmark and the DG P1-P1 schemes ®rst a constant and second an adaptively modi®ed time step is used. For comparison also the solution for the velocity vz at the end of the cantilever using modal decomposition of the FE model was computed. In order to get a physically motivated change of the time step size, the modal damping parameters a and b were chosen as: a ˆ 0:4, b ˆ 0:0014. The default resp. starting time step size for all computations was set to Dt ˆ 0:4; the integration interval considered is t 2 ‰0; 5Š. In Fig. 19 the evolution of the time step size is presented. It is expected, that the DG P1-P1 method would lead to a larger time step size than the Newmark method due to the higher order accuracy. However, as found

…35†

with NT ˆ T =Dtn as remaining number of time steps until the end of the process. Time step control by local error control: If a constant time step size is assumed, e.g. Dtnew ˆ cDtold , and the order of accuracy of the Newmark scheme is taken into account, we obtain: ke…cDt†ktotal;X;T ˆ c2 ke…Dt†ktotal;X;T :

…36†

Fig. 19. Evolution of the time step size in time for Newmark and DG P1-P1 methods.

2130

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131

Fig. 20. Velocity history at the end of cantilever with and without time step adaptation using DG P1-P1 and Newmark scheme compared to the modal solution. (a) Newmark scheme, (b) DG P1-P1 scheme.

also in other numerical examples, the Newmark scheme leads to a larger time step size. The main reason is that the jumps …Du; Dv) give only a restricted information about the error at a discrete time state, whereas the local error estimator for the Newmark scheme describes the error evolution in one time step. The time history plots of the velocity vz at the end of the clamped plate in Fig. 20a show the solution with time step adaptation for the Newmark method which approaches the modal solution very well. The major improvement of the adaptive solution compared to using a constant time step size is in the beginning of the process, when a small time step is needed. In the time range t 2 ‰2:5; 5:0Š the solution with constant time step size will also approximate the modal solution, which is caused by the physical damping. For the DG method only very small di€erences between the solutions with constant and variable time step were obtained, also indicating the higher accuracy properties for the fairly large time steps.

7. Conclusions The basic properties of two numerical time integration schemes, DG and Newmark method were presented and compared. It was also shown, that the DG P1-P1 scheme has some advantages compared to the DG P2 scheme within the DG family of methods. Then, two representative problems for wave propagation and for stationary vibration were chosen to discuss the properties of the numerical integration schemes. In particular two di€erent proposals for the DG P1-P1 method by Wiberg/Li and by Ruge were compared concerning ef®ciency. Further, the arti®cial dispersion introduced by spatial discretization was investigated on stationary and adaptively modi®ed FE-meshes. It was found that the DG P1-P1 scheme leads to rather good solutions in adaptively modi®ed meshes for short duration problems, where the numerical damping does not yet a€ect the global energy too much. Finally, two local error esti-

mators resp. indicators for the time integration error were discussed brie¯y using a numerical example. References [1] Aharoni D, Bar-Yoseph P. Mixed ®nite element formulations in the time domain for solution of dynamic problems. Comput Mech 1992;9:359±74. [2] Celep Z, Bazant ZP. Spurious re¯ection of elastic waves due to gradually changing ®nite element size. Int J Numer Meth Engng 1983;19:631±46. [3] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized a method. J Appl Mech 1993;60:371±5. [4] Estep D. A posteriori error bounds and global error control for approximations of ordinary di€erential equations. SIAM J Numer Anal 1995;32(1):1±48. [5] Gra€ KF. Wave motion in elastic solids. Ohio State University Press; 1975. [6] Hilber HM, Hughes TJR. Collocation, dissipation and `overshoot' for time integration schemes in structural dynamics. Earthquake Engng Struct Dyn 1978;6:99±117. [7] Ho€ C, Pahl PJ. Development of an implicit method with numerical dissipation from a generalized single-step algorithm for structural dynamics. Comput Meth Appl Mech Engng 1988;67:367±85. [8] Hughes TJR. The ®nite element method ± linear static and dynamic ®nite element analysis. Prentice-Hall: International Editions; 1987. [9] Hughes TJR, Hulbert GM. Space-time ®nite element methods for elastodynamics: formulations and error estimates. Comput Meth Appl Mech Engng 1988;62:339±63. [10] Hulbert GM. Discontinuity-capturing operators for elastodynamics. Comput Meth Appl Mech Engng 1992;96: 409±26. [11] Hulbert GM. Time ®nite element methods for structural dynamics. Int J Numer Meth Engng 1992;33:307±31. [12] Hulbert GM, Hughes TJR. Numerical evaluation and comparison of subcycling algorithms for structural dynamics. Technical Report, Stanford University, 1988. p. 1± 65. [13] Johnson C. Discontinuous Galerkin ®nite element methods for second order hyperbolic problems. Comput Meth Appl Mech Engng 1993;107:117±29.

J. Neumann, K. Schweizerhof / Computers and Structures 79 (2001) 2117±2131 [14] Johnson C, Pitk aranta J. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math Comput 1986;46(173):1±26. [15] Li X. Adaptive ®nite element procedures in structural mechanics. Doktorarbeit, Chalmers University of Technology, Gothenborg, 1996. [16] Li XD, Wiberg NE. Structural dynamic analysis by a timediscontinuous Galerkin ®nite element method. Int J Numer Meth Engng 1996;39:2131±52. [17] Borri M, Mello F, Atluri SN. Variational approaches for dynamics and time-®nite-elements: numerical studies. Comput Mech 1990;7:49±76. [18] Mindlin RD. In¯uence of rotatory inertia and shear in ¯exural motions of isotropic elastic plates. J Appl Mech 1951;18:31±8. [19] Newmark NM. A numerical method for structural dynamics. J Engng Mech Div ASCE 1959;85:67±94. [20] Nilsson L. Impact loading on concrete structures. PhD Thesis, Department of Structural Mechanics, Chalmers, Gothenburg, 1979. [21] Oden JT. A general theory of ®nite elements II Applications. Int J Numer Meth Engng 1969;1:247±59. [22] Riccius J, Schweizerhof K. Aspects of hierarchical hadaptive dynamic analyses. In: Topping BHV editor, Third

[23] [24]

[25] [26] [27] [28] [29]

2131

International Conference on Computational Structures Technology. Civil-Comp Press; 1996. Ruge P. Hybrid time-®nite-elements with time-step-adaption by discontinuity control. Comput Mech 1996;17:392± 7. Cirak F, Schleupen A, Ramm E. Error estimators for local and global variables using duality principles. Fourth World Congress on Computational Mechanics. Buenos Aires: IACM; 1998. Warburton GB. Formulae for errors for initial displacement and velocity problems using the Newmark method. Earthquake Engng Struct Dyn 1989;18:565±73. Zienkiewicz OC. A new look at the Newmark, Houbolt and other time stepping formulas. A weighted residual approach. Earthquake Engng Struct Dyn 1977;5:413±8. Zienkiewicz OC, Taylor RL. The ®nite element method ± vol 1. McGraw-Hill: International Editions; 1989. Zienkiewicz OC, Xie YM. A simple local error estimate and an adaptive time-stepping procedure for dynamic analysis. Earthquake Engng Struct Dyn 1991;22:871±87. Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates, part 1. The recovery technique. Int J Numer Meth Engng 1992;33:1131± 364.