Error estimation and adaptivity in explicit nonlinear finite element simulation of quasi-static problems

Computers and Structures 72 (1999) 627±644

www.elsevier.com/locate/compstruc

Error estimation and adaptivity in explicit nonlinear ®nite element simulation of quasi-static problems Kjell M. Mathisen*, Odd S. Hopperstad, Knut M. Okstad, Torodd Berstad Department of Structural Engineering, Norwegian University of Science and Technology, Rich. Birkelandsv. 1a, N±7034, Trondheim, Norway Received 16 October 1998

Abstract A program module for error estimation with application to nonlinear ®nite element (FE) analysis of shell structures is coupled with the adaptive solution procedure in the explicit FE code LS-DYNA. The error estimation module provides estimates of the local and global errors and element-level re®nement indicators. Hence, selective re®nement of the mesh in areas where the local error is relatively large compared with a user-de®ned tolerance is made possible. Furthermore, the relative global error is estimated giving a measure of the overall accuracy of the FE model. Projection-type error estimators based on the L2-norm of the stress vector and the accumulated plastic strain are used to predict the discretization error by comparison of the FE solution with an improved C 0continuous solution obtained by the SPR-method. Three example problems including both material and geometric nonlinearities are provided. The numerical results show that the error estimates capture phenomena such as di€use necking and local buckling, and give meshes with high resolution in areas with large deformations or high stress gradients. # 1999 Elsevier Science Ltd. All rights reserved.

1. Introduction Today, industry problems in engineering are associated with complex geometries, nonlinear constitutive behavior, contacts, large deformations, etc. For such problems only numerical methods, like the FE method, yield solutions of the associated mathematical model. Since numerical methods yield approximate solutions only, it is necessary to control the errors inherent in the method. During the last two decades much research has been devoted to development of error estimators and adaptive FE strategies for solving linear problems. When considering nonlinear problems adaptive sol-

* Corresponding author.

ution strategies become even more appealing, particularly for problems where the optimal mesh changes continuously throughout the analysis, since the CPU time consumption when solving for such problems may become prohibitively large. Another reason for moving towards automatic adaptive solution of nonlinear problems is the ever decreasing cost of computation together with the recent developments of more robust error estimators. The objective of mesh adaptivity is to obtain a FE mesh which is optimal in the sense that the computational costs involved are minimal under the constraint that the error in the solution is below a certain prescribed limit. Since the computational e€ort can be linked to the number of unknowns in the FE mesh the task is to ®nd a mesh with a minimum number of unknowns or nodes for a given error tolerance. In

0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 3 2 8 - 9

628

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

general, adaptive methods rely on error indicators and error estimators which can be computed a priori or a posteriori. Two main types of a posteriori estimates exist today: 1. Residual based estimators, in which the error can be estimated by some norm of the residual of the governing equations, see BabusÏ ka and Rheinboldt [1] and BabusÏ ka et al. [2]. 2. Recovery based estimators, in which the error is estimated by taking the di€erence between a recovered (post-processed higher order) solution and the FE solution, see Zienkiewicz and Zhu [3±5]. Based on the computed errors a new re®ned mesh can be constructed which yields a better approximate solution. To obtain an optimal mesh it is desirable to design the mesh such that the error contributions of the elements are equidistributed over the mesh. Both the residual based and the recovery based error estimators are frequently used estimators for elastic problems. According to Rheinboldt [6] and Stein et al. [7] the introduced error estimators for linear elastic problems may be extended into the geometrically nonlinear regime if they are applied to the linearized problem. When considering explicit nonlinear ®nite element analysis in general, theoretical proofs concerning the convergence that were established in the elastic range are no longer valid. Since we are considering quasi-static problems only, a ®rst step in order to justify the use of recovery-based Zienkiewicz±Zhu error estimators [4] is to neglect the errors due to the discretization in time. In prior studies, an overview of key aspects which have to be considered in an adaptive nonlinear solution procedure has been presented, along with the implementation of an h-adaptive mesh re®nement technique into a general purpose implicit nonlinear FE program for solving geometrically nonlinear shell type problems involving linear elastic materials [8±10]. In the case of history-dependent materials, there exist only a few contributions on adaptive strategies for history-dependent nonlinear problems. Johnson and Hansbo [11] developed an adaptive strategy for small strain elasto-plasticity using the Hencky model. However, due to restrictions in the Hencky model, this method is generally not applicable to standard elastic± plastic and elastic±viscoplastic history-dependent problems. Other contributions involving history-dependent materials include the work by Cheng and Kikuchi [12], Ortiz and Quigley [13], Lee and Bathe [14], Peric et al. [15, 16], Wriggers and Scherf [17], Gallimard et al. [18] and Barthold et al. [19, 20]. In this paper, the objective is to present ecient (although not theoretically justi®ed) techniques for controlling the discretization error in quasi-static simu-

lations by explicit nonlinear ®nite element calculations. It is demonstrated that Zienkiewicz±Zhu type error estimates can be used to provide accurate resolution of the current critical zone being formed. 2. Error estimation 2.1. Preliminaries Among the various sources of error that can be identi®ed in explicit FE simulation of quasi-static problems, the following are of great importance: 1. The discretization error: spatial discretization errors are inevitably connected with the application of the FE method, because of the discrete representation of the computational domain. In that way any FE mesh is associated with this type of error. For the particular case of linear FE analysis of a regular problem the solution converges towards the exact one and the error decreases to zero when the element size decreases. 2. Error in the choice of the constitutive model: this gives the di€erence between the real behavior of the material and the one represented by the constitutive law. In addition there may be errors in the numerical values of the chosen constitutive law parameters. 3. Error due to numerical integration of the constitutive law: this type of error is characteristic for material nonlinear analysis. As material equations are written in rate form they must be integrated in time with an appropriate integration scheme. The choice of the integration scheme performed at every step (explicit, implicit) is fundamental for the accuracy. In this study, only estimation of spatial discretization error is considered. Two major questions should be answered while analyzing a solution obtained with a given FE mesh: 1. What is the precision of the computed solution in global (structural) and local (elemental) sense, i.e. how big is the error? 2. How to obtain a solution within a prescribed accuracy at a minimum computational cost? Error estimation and adaptive mesh re®nement give the answer to the two questions above. The aim is to estimate the error of the FE method, using already computed data. A posteriori estimates employ the solution obtained by the FE analysis in addition to some a priori assumptions about the solution. Without an a posteriori error estimate there is no reliable way of judging the quality of the solution. If the analyst's intuition fails when designing the FE mesh, it is also

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

likely that it fails when evaluating the precision of the computed solution. Therefore the philosophy of an adaptive process is that the computer (code) decides on where to re®ne the mesh and does so automatically. There are two major tasks, closely tied together, that should be realized in order to obtain the answers to the two preceding questions: 1. Establish a posteriori local and global error estimates which accurately and reliably characterize the error in the FE solution according to some speci®ed norm. 2. Design an adaptive procedure for automatic mesh re®nement in order to achieve a solution with a prescribed accuracy in an e€ective and economic manner. The ®rst task is addressed below, while the latter is considered in the next section. 2.2. Error norms The FE solution, in terms of the displacements uh and stresses s h , di€er from the corresponding exact values, u and s . The di€erence de®nes the point-wise errors eu ˆ u ÿ uh

…1†

es ˆ s ÿ s h

…2†

These point-wise de®nitions of errors are generally dif®cult to specify, and various integral measures are more conveniently adopted. One of the most common of such measures is the energy norm … kekE ˆ

O

eTs Cÿ1 es

1 dV

2

…3†

A more direct measure is the so-called L2-norm, which can be associated with the error in any quantity. Thus, for the displacements and stresses, respectively, the errors in L2-norm are … keu kL2 ˆ

O

… kes kL2 ˆ

O

eTu eu dV

eTs es

1 2

…4†

1 dV

2

…5†

The latter expression di€ers from the energy norm only by a weighting matrix Cÿ1. Although the norms (3)±(5) here are de®ned on the whole domain, O, we note that the square of each one

629

can be obtained by summing element contributions, i.e. kek2 ˆ

nel X eˆ1

kek2e

…6†

where keke represents the contribution from element e and is obtained by integrating over the elemental subdomain, Oe, instead of O. Although the value of the energy-norm error does have a physical meaning (indeed, it measures the absolute error in strain energy), it may be easier to interpret relative quantities, especially for L2-norms. For this purpose we adopt the relative percentage error Zˆ

kek
100% kuk

…7†

where kuk denotes the corresponding norm of the exact solution of the problem. The percentage error (7) can be evaluated for the whole domain and for each elemental sub-domain. 2.3. Zienkiewicz±Zhu error estimation Usually, the FE solution is less smooth than the exact one. In other words, it has some de®ciency: . The degree of the interpolation polynomials of the elements may be insucient to represent the exact solution which may be a higher order polynomial or some other function. . The derivatives of the FE solution (strains and stresses) are discontinuous across the element interfaces (see Fig. 1a), as C 0 interpolation functions are usually employed for the FE approximation of displacements, while the derivatives of the exact solution generally are continuous. These two features of the FE solution may constitute a useful indication to estimate the discretization error. In order to construct error estimators based on these indications we must recover some higher order solution. A simple and very ecient a posteriori error estimator for practical engineering analysis is proposed by Zienkiewicz and Zhu [3] and thus frequently referred to as the ZZ-estimator. The basic idea consists in the use of a continuous stress ®eld, s  , which eliminates the stress jump across the element interfaces (see Fig. 1b), and is (at least) one order more accurate than the FE stress ®eld, s h . In the ZZ-estimator the point-wise error (2) is approximated by the di€erence between the higher order solution and the FE solution, viz. es ˆ s  ÿ s h

…8†

630

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 1. Stress ®eld on a four-element patch. (a) Discontinuous FE ®eld; (b) recovered C 0-continuous ®eld.

By substituting the estimate es into Eq. (3), we thus obtain the ZZ-error estimator in global energy norm ke kE ˆ

… O

…s s ÿ s h †T Cÿ1 …s s ÿ s h † dV

1 2

…9†

Obviously, the quality and reliability of the ZZ-estimator depend on the accuracy of the recovered solution s  . Zienkiewicz and Zhu [21, 22] show that if the recovered stress ®eld is superconvergent, the ZZ-error estimator will always be asymptotically exact in energy norm. 2.4. Error estimators for nonlinear materials Two error estimators of ZZ-type are considered herein. Each estimator measures the discontinuity across the elements in some quantity that is derived from the approximate FE solution. In the present study, the following criteria have been used in the initial selection of potential error estimators for comparative evaluation. These criteria state that the error estimator must be: 1. Valid for the material models considered. 2. Valid for large displacements (and rotations) but small strains. 3. Easy to compute (i.e. with additional necessary computations kept to a minimum). 4. Based on results available from the analysis code (LS-DYNA [23]). In the literature several possible error estimates have been considered. Tetambe et al. [24] tested several L2normsÐthe L2-norms in stresses, total strains, equivalent total strain and total strain rate, together with the rate energy norm. Further, Peric et al. [15] tested error estimates based on incremental plastic work, plastic dissipation and the energy norm. With LS-DYNA we

are limited to using error estimates based on either rate form quantities, such as strain rate or plastic strain rate, or an estimate based on stresses or accumulated plastic strain. Therefore, we focus on the L2norm errors in stresses and accumulated plastic strain. These are also found to be good error estimates by the authors mentioned above. Adopting the ZZ-approach, the point-wise error in accumulated plastic strain, p (see [25]), may be approximated by ep ˆ p ÿ ph

…10†

We then de®ne the following two error estimators: 1. L2-norm error in stresses kes kL2 ˆ

… O

…s s ÿ s h †T …s s ÿ s h † dV

1 2

…11†

2. L2-norm error in accumulated plastic strain kep kL2

… ˆ



O

h 2

… p ÿ p † dV

1 2

…12†

2.5. Superconvergent patch recovery The ZZ-error estimator was originally used in combination with global least-squares ®t methods and direct nodal averaging for obtaining the higher order ®eld. More recently, a local least-squares ®t approach, popularly denoted superconvergent patch recovery (SPR) has proven to be more ecient and to give more accurate results. Herein, we use the SPR approach in its original form as proposed by Zienkiewicz and Zhu [4], but extended to shell problems.

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 2. Examples of nodal patches: (w) nodal points; (*) the patch node (node de®ning the nodal patch); () result sampling points (superconvergent points).

The SPR procedure utilizes the concept of superconvergent points. These are certain locations within each element where the computed FE solution has higher accuracy than elsewhere. The rate of convergence is also higher in these points than at other locations. The ®rst demonstration of superconvergence is due to Barlow [26]. It is also discussed by Hinton and Campbell [27] and later by Barlow [28] for distorted elements. For quadrilateral ®nite elements, the superconvergent points for stresses and strains coincide with the reduced Gauss quadrature points, also referred to as Barlow points. Linear elements have thus one superconvergent point each located at the centroid (see Fig. 2). The Barlow points are generally the best sampling points for stresses whereas the element nodes are the worst, as the derivatives of the shape functions tend to behave poorly near the border of the interpolation region. Let vh denote some piece-wise continuous ®eld, i.e. h v ( C 0(O) but vhvOe $ C 0(Oe), which is computed by the analysis code and that we want to recover an improved version of. The corresponding recovered ®eld, which we denote v, is expressed as a polynominal v ˆ Pa

…13†

where P is a matrix of polynormal terms and a is a vector of unknown coecients. For four-noded quadrilateral elements, we use the bilinear basis 3 2 P1 7 6 .. 7; Pi ˆ ‰1,x,y,xyŠ Pˆ6 …14† 5 4 . Pn

1 Boundary nodes are normally not connected to more than two elements. With ®rst-order elements containing one superconvergent point each, there must be at least four elements in the patch in order to be able to recover a bilinear ®eld over the patch. For interior three-element patches (Fig. 2b), we omit the xy-term from Eq. (14) to avoid rank singularity of the resulting equation system.

631

with n denoting the dimension of the vector ®eld v. According to Zienkiewicz and Zhu [4], slightly better results may be obtained for quadrilateral elements when the same polynomial terms as the basis functions used to interpolate the primary variables of the FE problem are used instead of only the completely linear basis [1, x, y ], although the latter choice has the advantage of being independent of the local coordinate system used. However, by carefully selecting a local coordinate system for each patch, we avoid problems with singular or poorly conditioned equation systems which may occur when the bilinear basis is used. This issue is addressed further by Kvamsdal and Okstad [29], and the procedure proposed therein for the local patch coordinate system is also used in this study. The polynomial (13) is de®ned over small patches of elements surrounding each nodal point. The coecients a are then determined from a least-squares ®t of the ®eld v to the values of vh at the superconvergent points within the elements in the patch, i.e. we minimize the following functional with respect to a F…a† ˆ

nsp X …vhi ÿ v …xi ††T …vhi ÿ v …xi ††

…15†

iˆ1

where xi is the coordinates of the ith superconvergent point, vhi is the computed value of vh at that point, and nsp denotes the total number of superconvergent points in the patch. The minimization of Eq. (15) leads to a small linear system of equations which is assembled and solved on each patch of elements. Typical examples of such nodal patches are shown in Fig. 2. Note that patches are established for interior nodes only, as nodes on the exterior boundary rarely are connected to enough elements.1 Neither are patches established for nodes along shell intersections, since some stress components should not be continuous across such intersections, e.g. when two perpendicular shell surfaces meet, the normal stress in one shell surface is transferred into a transverse shear stress in the other surface. 2.6. Global recovered ®eld When result recovery is used for error estimation, we are normally interested in recovered values at element-interior points (integration points) rather than nodal values. Since a speci®c element belongs to more than one patch, the patch recovery does not provide a unique result value at such points. In order to construct a global recovered ®eld, Blacker and Belytschko [30] proposed conjoining the polynomial expansions (13) for all patches containing the actual element using the nodal shape functions as weighting functions. Adopting this approach, the recovered ®eld within a bilinear element e is evaluated through

632

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 3. Nodal patches associated with a node (.) on the exterior boundary. (a) Node connected to exterior element boundaries only. (b) Node connected to one interior element boundary. (c) Node connected to two interior element boundaries.

ve …x† ˆ

4 X aˆ1

N a …x†va …x†

…16†

where Na(x) denotes the bilinear shape function associated with element node a and va (x) is a local recovered ®eld of the form (13). If the element node a is an interior node, va (x) is evaluated on the patch of elements surrounding that node. For nodes lying on the exterior boundary, va (x) is instead evaluated on the patch (or patches) associated with the other node(s) that are connected to node a through an interior element boundary, see Fig. 3. If more than one patch in this manner is connected to a boundary node (Fig. 3c), the corresponding values for va computed on each patch are averaged. 3. Adaptivity within LS-DYNA 3.1. Adaptive solution procedure Fig. 4 gives an algorithmic overview of the built-in adaptive solution procedure in LS-DYNA [23]. Here, lt is the load parameter at time t, rt denotes the primary unknowns of the FE problem (e.g. nodal displacements) at time t, whereas st represents the other solution-dependent state variablesÐstresses, strains, etc. The total duration of a simulation is denoted ttotal whereas the time step is denoted dt. In an adaptive simulation the user speci®es a birth time, tbirth, at which the adaptive remeshing begins, the death time, tdeath, at which the adaptive remeshing ends, and a time interval, Dt, between each check for adaptive re®nement. In addition, the maximum number of re®nement levels, nmax, the minimum element size, hmin, and of course the prescribed error tolerance, Zmax, are speci®ed. Explicit FE analyses are usually characterized by very small (but many) time steps, dt. However, the computational cost for each step is mod-

est when compared with an implicit analysis. It is therefore not desirable to check for adaptive remeshing at every time step, and typically Dt >> dt. The time step dt is continuously adjusted by LS-DYNA in order to maintain a stable solution. LS-DYNA o€ers two types of adaptivityÐa onepass strategy and a two-pass strategy. In the latter variant, the solution procedure jumps back to the time instant tk (see Fig. 4) at which the previous error control was performed, after the adaptive step is completed. Note that this backward jump is performed also if the estimated error is within the tolerance and not only when the mesh has been re®ned. Thus, the two-pass strategy ends up calculating each time interval Dt between tbirth and tdeath twice and is therefore considerably more expensive than the one-pass strategy. 3.2. Adaptive re®nement strategy The error estimation including calculation of element-wise re®nement indicators is herein performed by a separate program module [31, 32] in which the procedures described in Section 2 have been implemented. LS-DYNA has been tailored by adding a system call to the error estimation module. The adaptivity is based on the so-called h-method, where the size of the elements in the re®ned mesh is de®ned through the expression hnew ˆ hold e e =fe

…17†

where h old denotes the characteristic size of element e e in the existing (old) mesh, and h new is the desired size e of the new elements in the region covered by that old element. The mesh re®nement parameter, fe, is computed from the estimated element errors based on some mesh optimality criterion. The optimality criterion aims to minimize the overall

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 4. Schematic overview of the adaptive solution procedure in LS-DYNA.

633

634

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

computational cost of a simulation for a given accuracy of the solution results. Herein, we adopt the strategy presented in [33, 34], where one tries to minimize the number of elements with the constraint that the total error should be equal to the prescribed tolerance.2 The resulting expression for the re®nement indicator reads fe ˆ

Ze1=…q‡1† Z1=q max

nel X eˆ1

!1=2q Z2=…q‡1† e

…18†

where Ze represents the contribution from element e to  the estimated p el 2 relative global relative error, Z , i.e. Z = …Sneˆ1 Ze †, and q is the theoretical convergence rate of the elements. Thus, for linear elements q=1. The mesh is re®ned by pure sub-division of the existing elements into four new elements introducing hanging nodes between the re®ned and un-re®ned elements. Details on the mesh re®nement procedure may be found in [35]. Note that the mesh is re®ned only one level per adaptive step, regardless of the value on fe. The only two values of fe that we are able to satisfy exactly is therefore 1.0 (no sub-division) and 2.0 (subdivision). Herein, we choose an in-between solution and subdivide all elements for which fe r 1.5. Fusion of elements (when fe < 1) is not available in the current version of LS-DYNA [23].

4. Numerical studies To demonstrate the use of error estimation in adaptive FE plasticity analysis we consider three elastic± plastic example problemsÐa tensile test, stretch bending of a square hollow beam, and axial crushing of a thin-walled square pro®le. All simulations are performed with LS-DYNA using the Belytschko±Lin± Tsay shell element [36] with single-point in-plane integration, ®ve integration points through the thickness and sti€ness-based hourglass control. The material is modeled as an isotropic elastic±plastic material, adopting the von Mises yield criterion, an isotropic strain hardening rule, and the associated ¯ow rule in the plastic domain. The strain hardening of the material is de®ned by 2

It should be noted that this criterion is `optimal' strictly for static or implicit computations only. For explicit computations the overall solution time is governed by the number of elements and the size of the smallest element in the mesh, since the latter controls the time step size through the CFL condition. An optimality criterion for explicit computations should therefore also try to maximize the size of the smallest elements.

Y ˆ Y0 ‡ Q1 …1 ÿ exp…ÿC1 p††

…19†

‡ Q2 …1 ÿ exp…ÿC2 p††

where Y is the von Mises equivalent stress, p denotes the accumulated plastic strain, and Y0, Qk and Ck are material constants. For each test we carry out one reference simulation using a ®ne uniform mesh, and three adaptive simulations with characteristics as given in Table 1. The computational cost of the adaptive simulations with respect to the reference simulation is given in Table 2 for the three problems. Tables 3 and 4 summarize the adaptivity control parameters and the material parameters, respectively, that are used for each example problem. The prescribed tolerances Zp,max and Zs,max that are used to guide the adaptive simulations, are here selected based on the estimated relative global errors of the corresponding reference simulation. LS-DYNA is an explicit FE code using the conditionally stable central di€erence method for time integration of the discrete equations of motions. In general, stability requires time steps in the order of micro seconds, and simulation of quasi-static problems is practically dicult. One possible way to overcome this problem is to increase the loading rate in the simulation provided that inertia and rate e€ects still are negligible in the response. Another method is to scale up the mass to increase the stable time step. The ®rst method is adopted for the stretch bending simulation whereas mass scaling is used in the simulations of the tensile test and the axial crushing test. The mass scaling is here performed by scaling the density by a factor which is kept constant throughout the entire simulation.

Table 1 Adaptive simulations Simulation no.

Error estimator

Adaptive strategy

1 2 3

Zp Zp Zs

One-pass Two-pass Two-pass

Table 2 Normalized CPU time of the simulations Problem

Reference

No. 1

No. 2

No. 3

Tensile test Stretch bending test Axial crushing test

1.0 1.0 1.0

0.89 0.96 0.66

2.04 1.61 0.85

0.87 3.27 1.13

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

635

Table 3 Control parameters for adaptivity Problem

ttotal (s)

tbirth (s)

tdeath (s)

Dt (s)

nmax

hmin (mm)

Zp,max (%)

Zs,max (%)

Tensile test Stretch bending test Axial crushing test

35 0.008 20

0 0 0

1 1 1

0.5 0.0001 0.4

4 4 4

0 0 0

10 10 20

15 10 25

Table 4 Material constants used in the simulations Problem

E (MPa)

v

Y0 (MPa)

Q1 (MPa)

C1

Q2 (MPa)

C2

Tensile test Stretch bending test Axial crushing test

72,000 72,000 70,000

0.33 0.33 0.33

300 300 75

300 300 75.6

8 8 25

0 0 200

0 0 1.9

4.1. Tensile test Fig. 5 shows the geometry of the test specimen. In the simulations, the relative velocity of the bolts is prescribed, increasing the velocity smoothly from zero to a constant level in order to avoid spuriously high frequency oscillations in the solution. Fig. 6 depicts the FE mesh used in the reference simulation, and the initial mesh of the adaptive simulations. The two meshes contain 1316 and 256 elements, respectively. Fig. 7 shows the deformed geometry at four di€erent displacement levels for the simulation using the reference mesh (Fig. 6a) and for the three adaptive simulations. The force±displacement curves for the four simulations are compared in Fig. 8. The estimated relative global errors in the reference solution are plotted against the displacement level in Fig. 9. Both error estimates are less than 10% in the plastic domain prior to necking, but Zs is signi®cantly

higher than Zp. The relative global error increases (Zp more strongly than Zs) as strain localization occurs in the specimen's gage section, hence indicating that even local errors in the solution may be detected by studying the evolution of the relative global error. The peak in Zp, occurring when the specimen is mainly in the elastic domain, is due to the normalization procedure. The relative global error is the sum over all elements of the ratio of the estimated error to the norm of the FE (or improved) solution. At the start of the deformation process the accumulated plastic strain is small or even zero, and the result may be large values of the relative error, even if the absolute error is small. Fig. 8 shows that all simulations give the same force±displacement curve up to necking. In the descending part of the force±displacement curve, the reference solution and adaptive simulation no. 2 compare

Fig. 5. The tensile test: geometry of the test specimen (all units in mm).

636

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 6. The tensile test: (a) ®ne mesh used in the reference simulation; (b) initial mesh of the adaptive simulations.

well, while adaptive simulations nos 1 and 3 predict higher load levels. The element distribution (Fig. 7) and the total number of elements (Fig. 10) during the simulation di€er signi®cantly between the di€erent schemes. In all cases, the element mesh is re®ned in the area around the loading point, where the stress gradients are large and the stress state is multi-axial. Both error estimates capture the necking of the specimen and resolve the high gradients in the localization zone when two-pass adaptivity (nos 2 and 3) is used; this is not the case for onepass adaptivity (no. 1). For the adopted values of Zs,max and Zp,max, the error estimator Zp gives the highest mesh density and CPU time, see Table 2.

Fig. 7. The tensile test: deformed geometry at displacement levels 12.5, 16.5, 20.5 and 24.5 mm. (a) Reference solution; (b) adaptive simulation no. 1; (c) adaptive simulation no. 2; (d) adaptive simulation no. 3.

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 8. The tensile test: force±displacement curves.

4.2. Stretch bending test Fig. 11 shows a simpli®ed model of the experimental stretch bending set-up presented by Welo et al. [37]. The model consists of a rigid bending die with radius 50 mm, and an extruded aluminum square hollow section with nominal wall thickness 2 mm, width 20 mm and length 440 mm. The model is symmetric about the plane x=220 (mm) and about the plane y=0. The pro®le is supported along the y-axis where the displacements in the x-, y- and z-directions and the rotations about the x- and z-axes are ®xed, whereas the rotation

Fig. 9. The tensile test: estimated relative global errors.

637

Fig. 10. The tensile test: number of elements in the mesh during the adaptive simulations.

about the y-axis is free. The pro®le is ¯exible with elastic±plastic material properties between the plane x=135 (mm) and the plane x=220 (mm). Elsewhere it is rigid. The pro®le is loaded by gradually moving the bending die downwards. The maximum displacement in the z-direction is 60 mm. Fig. 12 shows the FE mesh used in the reference simulation and the initial mesh of the adaptive simulations. Symmetry is taken into account in the model and only one quarter of the bending die and the extruded pro®le is discretized. The reference model consists of 1420 shell elements, whereas the initial mesh of the adaptive simulations consists of 384 shell elements. The default version of the surface-to-surface contact algorithm in LS-DYNA is used in the simulations. Friction between the bending die and the pro®le is neglected. As for the tensile test, one reference simulation and three di€erent adaptive simulations with varying error estimator and adaptive strategy are conducted. The results are presented in Figs. 13±16, respectively, as ¯ange sagging against die displacement, deformed mesh at 60 mm die displacement, estimated relative global errors against die displacement (reference simulation only) and ®nally the total number of elements against die displacement. The ¯ange sagging is de®ned (see Fig. 11) as the displacement in the z-direction of point B relative to point A after stretch bending. Fig. 13 shows that the sagging developments predicted in the reference solution and the adaptive solutions are in good agreement. Fig. 15 shows that the relative global errors are less than 10% in the major part of the forming operation, and in this case the two error estimates give relatively consistent values of the relative global error.

638

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 11. The stretch bending test: geometry, boundary and loading conditions (all units in mm).

The element distribution (Fig. 14) and the total number of elements during the simulation (Fig. 16) are almost the same in the one-pass and two-pass adaptive simulations using the error estimator Zp, i.e. adaptive simulations nos 1 and 2. The required number of elements is signi®cantly higher in adaptive simulation no. 3 in which Zs is used, and this is also re¯ected by a signi®cant increase in the computational cost as shown in Table 2. 4.3. Axial crushing test Fig. 17 shows the specimen geometry, support and loading conditions in the axial crushing test of thinwalled square aluminum extrusions [38, 39]. The FE meshes used in the reference and adaptive simulations are depicted in Fig. 18. The specimen is clamped at the lower end, free at the top, and loaded quasi-statically

Fig. 12. The stretch bending test: (a) ®ne mesh used in reference simulation; (b) initial mesh of the adaptive simulations.

in compression. The load and the relative displacement between the base and the crosshead are measured. Owing to symmetry, only one quarter of the test setup is modeled. The specimen is modeled using shell elements. A rigid body modeled by brick elements is used to impose the load at the free end of the pro®le. The rigid body is given a constant velocity during the deformation. The specimen is modeled using, respectively, 240 and 2400 shell elements in the coarse and ®ne model. Initial imperfections are introduced near the top section of the model to trigger a symmetric buckling mode, see Fig. 18. The results from the reference solution are given in Figs. 19 and 21, where the absorbed energy (computed as the work of the external load) and the relative

Fig. 13. The stretch bending test: ¯ange sagging against die displacement.

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

639

Fig. 14. The stretch bending test: deformed geometry at 60 mm die displacement. (a) Reference solution; (b) adaptive simulation no. 1; (c) adaptive simulation no. 2; (d) adaptive simulation no. 3.

640

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 15. The stretch bending test: estimated relative global errors for the reference simulation.

global errors are plotted against the displacement of the rigid body. In this case, the relative global errors di€er signi®cantly: the error estimators Zp and Zs give relative global errors of approximately 12% and 20± 25%, respectively.

Fig. 16. The stretch bending test: number of elements in the mesh during the adaptive simulations.

The results from the three adaptive simulations are presented in Figs. 19, 20 and 22, in terms of absorbed energy against displacement, deformed mesh at 200 mm displacement, and number of elements against displacement, respectively.

Fig. 17. The axial crushing test: geometry, boundary and loading conditions (all units in mm).

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

641

The energy absorption and the deformation mode are similar in the reference solution and the adaptive solutions. The required numbers of elements are signi®cantly higher in adaptive simulations nos 1 and 3 than in no. 2. Hence, in this case the CPU time (see Table 2) increases by about 30% only, when going from one pass to two-pass adaptivity, using the error estimator based on plastic strain. 5. Concluding remarks

Fig. 18. The axial crushing test: (a) ®ne mesh used in the reference simulation; (b) initial mesh of the adaptive simulations.

Fig. 19. The axial crushing test: absorbed energy versus displacement.

A program module for error estimation of nonlinear FE shell analysis was coupled with the adaptive solution procedure in the explicit FE code LS-DYNA. Automatic re®nement of the mesh in areas where the element level error exceeds a user-de®ned tolerance is thus made possible. The discretization error is predicted by comparing the FE solution with an improved C 0-continuous solution obtained by the superconvergent patch recovery (SPR) method. Error estimators based on the L2-norm of the stress vector and the accumulated plastic strain are adopted. The adaptive procedure in LS-DYNA is based on mesh enrichment where a single `parent' quadrilateral element is split into four equally sized `sibling' elements. The data transfer between the old deformed mesh and the new mesh is trivially solved by letting the sibling elements inherit the state of the parent element. One-pass and two-pass adaptivity are available. One-pass adaptivity means that the analysis is restarted with the enriched mesh from the time step of the current error estimation. In two-pass adaptivity, the analysis is always restarted from the time step of the previous error estimation, even if the mesh remains unaltered. In practice, this means that the simulation is run through twice. Example problems involving material and geometric nonlinearities and quasi-static loading conditions were solved with LS-DYNA. It was found that both the stress- and strain-based error estimators produced meshes with high resolution in areas with high stress gradients or large deformations. It is however dicult to compare the accuracy and e€ectivity of the two error estimators since they predict di€erent relative global errors. In addition, the mutual relationship between relative global errors may change during the simulation. One-pass adaptivity is obviously more e€ective than two-pass adaptivity, but its robustness is less since abrupt changes in the response (i.e. due to geometric or material instabilities) between two error checks may be missed. However, a modi®ed adaptive procedure in which the analysis is restarted from the time step of the previous error check only when the current error check leads to mesh enrichment is readily implemented. Such a procedure will always give the

642

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

Fig. 20. The axial crushing test: deformed geometry at 200 mm displacement. (a) Reference solution; (b) adaptive simulation no. 1; (c) adaptive simulation no. 2; (d) adaptive simulation no. 3.

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

643

Acknowledgements This work was supported ®nancially by Hydro Raufoss Automotive Research Centre. The support is highly acknowledged.

References

Fig. 21. The axial crushing test: estimated relative global errors for the reference simulation.

same accuracy and robustness as two-pass adaptivity and in many cases be almost as e€ective as one-pass adaptivity. The main problem inherent in the adaptive simulation is the user-de®ned element-level error tolerance. It is not easy to specify a reasonable value for this tolerance without any a priori knowledge of the problem at hand. Such knowledge may be achieved by running a nonadaptive reference solution using a relatively coarse mesh. However, in the problems considered in this paper, a relative global error between 10 and 20% seems to be a good choice.

Fig. 22. The axial crushing test: number of elements in the mesh during the adaptive simulations.

[1] BabusÏ ka I, Rheinboldt WC. A-posteriori error estimates for the ®nite element method. International Journal for Numerical Methods in Engineering 1978;12:1597±615. [2] BabusÏ ka I, Strouboulis T, Upadhyay CS. A model study of the quality of a posteriori error estimators for linear elliptic problems. Error in the interior of patchwise uniform grids of triangles. Computer Methods in Applied Mechanics and Engineering 1994;114:307±78. [3] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedures for practical engineering analysis. International Journal for Numerical Methods in Engineering 1987;24:337±57. [4] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique. International Journal for Numerical Methods in Engineering 1992;33:1331±64. [5] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 2: error estimates and adaptivity. International Journal for Numerical Methods in Engineering 1992;33:1365±82. [6] Rheinboldt WC. Error estimates for nonlinear ®nite element computations. Computers and Structures 1985;20:91±8. [7] Stein E, Carstensen C, Seifert B, Ohnimus S. Adaptive ®nite element analysis of geometrically nonlinear plates and shells, especially buckling. International Journal for Numerical Methods in Engineering 1994;37:2631±55. [8] Okstad KM, Mathisen KM. Towards automatic adaptive geometrically nonlinear shell analysis. Part I: Implementation of an h-adaptive mesh re®nement procedure. International Journal for Numerical Methods in Engineering 1994;37:2657±78. [9] Okstad KM. Adaptive methods for nonlinear ®nite element analysis of shell structures. Dr Ing. dissertation, Department of Structural Engineering, The Norwegian Institute of Technology, Trondheim 1994. [10] Mathisen KM, Okstad KM. Interactive±adaptive geometrically nonlinear analysis of shell structures. Engineering with Computers 1996;12:63±83. [11] Johnson C, Hansbo P. Adaptive ®nite element methods in computational mechanics. Computer Methods in Applied Mechanics and Engineering 1992;101:143±81. [12] Cheng JH, Kikuchi N. A mesh rezoning technique for ®nite element simulations of metal forming processes. International Journal for Numerical Methods in Engineering 1986;23:219±28. [13] Ortiz M, Quigley JJ. Adaptive mesh re®nement in strain localization problems. Computer Methods in Applied Mechanics and Engineering 1991;90:781±804.

644

K.M. Mathisen et al. / Computers and Structures 72 (1999) 627±644

[14] Lee NS, Bathe K-J. Error indicators and adaptive remeshing in large deformation ®nite element analysis. Finite Elements in Analysis and Design 1994;16:99±139. [15] Peric D, Yu J, Owen DRJ. On error estimates and adaptivity in elastoplastic solids: applications to the numerical simulation of strain localization in classical and Cosserat continua. International Journal for Numerical Methods in Engineering 1994;37:1351±79. [16] Peric D, Hochard Ch, Dutko M, Owen DRJ. Transfer operations for evolving meshes in small strain elastoplasticity. Computer Methods in Applied Mechanics and Engineering 1996;137:331±44. [17] Wriggers P, Scherf O. Adaptive ®nite element methods for contact problems in plasticity. In: Proceedings of the 4th International Conference on Computational Plasticity, Barcelona, 1995. p. 787±807. [18] Gallimard L, LadeveÁze P, Pelle JP. Error estimation and adaptivity in elasto-plasticity. International Journal for Numerical Methods in Engineering 1996;39:189±217. [19] Barthold F-J, Stein E. Error estimation and mesh adaptivity for elasto-plastic deformations. In: Proceedings of the 5th International Conference on Computational Plasticity, Barcelona, 1997. p. 597±603. [20] Barthold F-J, Schmidt M, Stein E. Error indicators and mesh re®nements for ®nite-element-computations of elastoplastic deformations. In: Proceedings of 10th Nordic Seminar on Computational Mechanics, Tallinn, Estonia, 1997. p. 12±37 (abstracts). [21] Zhu JZ, Zienkiewicz OC. Superconvergence recovery technique and a posteriori error estimators. International Journal for Numerical Methods in Engineering 1990;30:1321±39. [22] Zienkiewicz OC, Zhu JZ. Adaptivity and mesh generation. International Journal for Numerical Methods in Engineering 1991;32:783±810. [23] Livermore Software Technology Corporation LS-DYNA Keywords User's Manuals, Version 940 1997. [24] Tetambe RP, Yunus SM, Rajakumar C, Saigal S. Examination of ¯ux projection-type error estimators in nonlinear ®nite element analysis. Computers and Structures 1995;54:641±53. [25] Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge: Cambridge University Press, 1990. [26] Barlow J. Optimal stress locations in ®nite element models. International Journal for Numerical Methods in Engineering 1976;10:243±51. [27] Hinton E, Campbell JS. Local and global smoothing of discontinuous ®nite element functions using a least

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

squares method. International Journal for Numerical Methods in Engineering 1974;8:461±80. Barlow J. More on optimal stress points, reduced integration, element distortions and error estimation. International Journal for Numerical Methods in Engineering 1989;28:1487±504. Kvamsdal T, Okstad KM. Error estimation based on superconvergent patch recovery using statically admissible stress ®elds. International Journal for Numerical Methods in Engineering 1998;42:443±72. Blacker T, Belytschko T. Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements. International Journal for Numerical Methods in Engineering 1994;37:517±36. Okstad KM, Kvamsdal T. AFEMÐadaptive ®nite element module, maintenance manual. Technical report, SINTEF Applied Mathematics, Department of Numerical Simulation, Trondheim 1997. Hopperstad OS, Mathisen KM, Okstad KM, Berstad T, Tiller I. Error estimation in plastic forming simulations. Technical report R-21-97, Norwegian University of Science and Technology, Trondheim 1997. LadeveÁze P, Pelle J-P, Rougeot Ph. Error estimation and mesh optimization for classical ®nite elements. Engineering Computations 1991;8:69±80. Coorevits P, LadeveÁze P, Pelle J-P. An automatic procedure with a control of accuracy for ®nite element analysis in 2D elasticity. Computer Methods in Applied Mechanics and Engineering 1995;121:91±120. Belytschko T, Wong BL, Plaskacz EJ. Fission±fusion adaptivity in ®nite elements for nonlinear dynamics of shells. Computers and Structures 1989;33:1307±23. Belytschko T, Lin JI, Tsay CS. Explicit algorithms for the nonlinear dynamics of shells. Computer Methods in Applied Mechanics and Engineering 1984;42:225±51. Welo T, Paulsen F, Moen KE. On the use of numerical simulation for predicting bendability and geometric tolerances of AA7xxx extrusions formed in stretch bending. In: Proceedings of Formability '94, Ostrava, Czech Republic, 1994. p. 24±7. Langseth M, Hopperstad OS. Static and dynamic axial crushing of square thin-walled aluminium extrusions. International Journal for Impact Engineering 1996;18:949±68. Berstad T. Material modelling of aluminium for crashworthiness analysis. Dr Ing. dissertation, Norwegian University of Science and Technology, Trondheim 1996.