Journal of Materials Processing Technology 105 (2000) 110±118
Comparison of implicit and explicit ®nite element methods for dynamic problems J.S. Sun, K.H. Lee, H.P. Lee
*
Department of Mechanical and Production Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore
Received 15 April 1999
Abstract
The ®nite element software ABAQUS offers several algorithms for dynamic analysis. The direct integration methods include the implicit and the explicit methods which can be used for linear and nonlinear problems. The performance of these two methods are compared for several dynamic problems including the impact of an elastic bar and a cylindrical disk on a rigid wall. The advantages of the implicit method for small wavefront wavefront problems and the explicit method for short transient transient problems problems are veri®ed. veri®ed. # 2000 Elsevier Science S.A. All rights reserved. Keywords: Implicit ®nite element method; Explicit ®nite element method; Dynamic problems
1. Introductio Introduction n
The ®nite element method (FEM) has become the most popular method in both research and industrial numerical simulations simulations.. Severa Severall algorithms algorithms,, with different different computacomputational costs, are implemented in the ®nite element codes, ABAQUS [1], which is a commonly used software for ®nite element analysis. Understanding the nature, advantages and disadvantages of these algorithms is very helpful for choosing the right algorithm for a particular problem. Comparison of implicit and explicit methods for ABAQUS in nonlinear problems has been reported by ReBelo et al. [2]. The unconditionally stable implicit method will encoun encounter ter some some dif®cu dif®cultie ltiess when when a compli complica cated ted threethreedimensional model is considered. The reasons are as follows: (i) as the reduction of the time increment continues, the comput computati ationa onall cost cost in the tangen tangentt stiffn stiffness ess matrix matrix is dramatica dramatically lly increased increased and even even causes causes diverge divergence; nce; (ii) local inst instab abili ilitie tiess caus causee forc forcee equi equili libr briu ium m to be dif® dif®cu cult lt to achieve. The explicit techniques are thus introduced to overcome the disadvantages disadvantages of the implicit method [3]. For the explicit method, the CPU cost is approximately proportional to the size of the ®nite element model and does not change as dramatica dramatically lly as the implicit implicit method. method.
The The draw drawba back ck of the the expl explic icit it meth method od is that that it is condi condition tionall ally y stable stable.. The stabil stability ity limit limit for an explic explicit it operator is that the maximum time increment must be less than a critical value of the smallest transition times for a dilatational wave to cross any element in the mesh. Secondly ondly,, the nature nature of the explici explicitt metho method d limits limits it to the analy analysis sis of short short transi transient ent proble problems. ms. If this this metho method d is used used for quasiquasi-sta static tic proble problems, ms, the inertia inertia effec effects ts must must be smal smalll enou enough gh to be negl neglec ecte ted. d. One One way way to assu assure re this is to set the limit on the kinematic energy to be less than 5% of the strain energy. Another limit is that only the ®rst-order, displacement methods elements (four-node quadrilaterals, eight-node bricks, etc.) are available for the present version. For dynamic problems, ABAQUS also offers some other methods such as a modal dynamic algorithm. However, only direct integration methods Ð implicit dynamics and explicit methods Ð are suitable for nonlinear problems. Most of the reported works on the comparison of implicit and explicit methods are on quasi-static nonlinear problems [2,3]. In this paper, attention is paid to comparison of the implicit and the explicit method for linear dynamic problems.
2. Solution Solution procedures procedures *
Corresponding author. Tel.: 65-874-2205; fax: 65-779-1459. E-mail address :
[email protected] (H.P. Lee).
The implic implicit it proced procedure ure uses uses an autom automati aticc increm increment ent strategy based on the success rate of a full Newton iterative
0924-0136/00/$ ± see front matter # 2000 Elsevier Science S.A. All rights reserved. PII: S0924-0136(00)00580-X
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
solution method: Du
i1
1 i i Dui KÀ t Á F À I
(1)
where Kt is the current tangent stiffness matrix, F the applied load vector, I the internal force vector, and Du is the increment of displacement. For an implicit dynamic procedure, the algorithm is de®ned by Hilber et al. [4,5]:
i1 1 aKui1 À aKui Fi1 Mu
u
i1
u
i
u
i
Dt u
i
2
Dt
12
i
À b u bu
i1
(3)
and u
i1
u Dt 1 À g
i
i1
gu
Dt
2 omax
where omax is the element maximum eigenvalue. A conservative estimate of the stable time increment is given by the minimum value for all elements. The above stability limit can be written as Dt
min
L e cd
where L e is the characteristic element dimension and cd the current effective, dilatational wave speed of the material. 3. Example problems
(4) 3.1. Impact of an elastic bar against a rigid wall
with b 14 1 À a2 Y
As mentioned above, the explicit integration operator is conditionally stable, so that the time increments must satisfy
(2)
where M is the mass matrix, K the stiffness matrix, F the vector of applied loads and u the displacement vector:
111
g 12 À aY
1 3
a0
(5)
a À0X05 is chosen by default in ABAQUS as a small damping term to quickly remove the high frequency noise without having a signi®cant effect on the meaningful, lower frequency response. The explicit procedure is based on the implementation of an explicit integration rule along with the use of diagonal element mass matrices. The equation of motion for the body is integrated using an explicit central difference integration rule: ui1 ui Dt i1 u i1 Y u i1a2 u iÀ1a2 12 Dt i1 Dt i ui
Dt
the acceleration, and i and i À 12, where u is the velocity, u i 12 refer to the increment number and mid-increment numbers.
i MÀ1 Á Fi À Ii u
As the ®rst example, the impact of an elastic bar against a rigid wall is presented. This model is used to compare the results by ABAQUS and the reported results by Zhong [6]. The dimensions of the bar are 1 Â 1 Â 10 m3 . Young's modulus is given as E 100 kPa. Possion's ratio is n 0, and the mass density is r 0X001kgmÀ3 . Eighty C3D8 and C3D8R (eight-node bricks) elements are used for the implicit and explicit methods, respectively. An initial velocity of 1 m sÀ1 is prescribed to the bar striking the rigid wall. The elements of the rigid wall are R3D4. The mesh is shown in Fig. 1. The impact time can be simply estimated by
(6)
where M is the mass ``lumped'' matrix, F the applied load vector and I the internal force vector.
2 L c
where L is the length of the bar, and c the wave speed for this kind of material, which is determined by À3 À1 E ar 2X0 Â 10 m s . For the implicit method, the half-step residual tolerance HAFTOL is a very important parameter to control the computational accuracy [1]. A large value of 10 is set initially. In Fig. 2a, the contact reaction forces are compared.
p
Fig. 1. Mesh of the bar.
1 1 2
J . S . S u n e t a l . / J o ur n a l o f M a t er i a l s P r o c e s s i n g T e c h n o l o g y 1 0 5 ( 2 0 0 0 ) 1 1 0 ±1 1 8
Fig. 2. (a) Contact reaction force RF1 with HAFTOL10; (b) the stress component S11 with HAFTOL 10; (c) contact reaction force RF1 with HAFTOL1; (d) the stress component S11 with HAFTOL 1.
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
113
Fig. 3. Mesh of the cylindrical disk.
It only takes 19 steps and less than 1 s CPU time for the explicit method, while it takes 106 s of CPU time for the implicit method. However, under this condition, the accuracy of the results by the explicit method is better than that by the implicit method. Fig. 2b shows the stress component S11. The values of S11 by the implicit method are found to oscillate much more than the corresponding results by the explicit method. From Fig. 2a, the reaction force is found to be about 10 N. In order to secure more accurate results, 10% of the reaction force, 1.0, is set as the revised HAFTOL value. The
analyze and predict failure. Comparing the stresses by these two methods is helpful for failure diagnosis. A three-dimensional example is needed for the above purpose. This example is shown in Fig. 3. A cylindrical disk is formed by two cylinders. The radius of the bottom cylinder is 1 m, and the radius of the top cylinder is 0.25 m. The bottom cylinder is of 0.5 m height, and the height of the top one is 0.1 m. The distance between the centers of these two cylinders is 0.31 m. There are 548 eightnode brick elements and 810 nodes for the mesh. Elements 173 and 409 are chosen for the purpose of stress comparison.
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
113
Fig. 3. Mesh of the cylindrical disk.
It only takes 19 steps and less than 1 s CPU time for the explicit method, while it takes 106 s of CPU time for the implicit method. However, under this condition, the accuracy of the results by the explicit method is better than that by the implicit method. Fig. 2b shows the stress component S11. The values of S11 by the implicit method are found to oscillate much more than the corresponding results by the explicit method. From Fig. 2a, the reaction force is found to be about 10 N. In order to secure more accurate results, 10% of the reaction force, 1.0, is set as the revised HAFTOL value. The reaction force and S11 are shown in Figs. 2c and d. The values of S11 by the implicit methods still oscillate more than those by the explicit method, but both the stresses and reaction force oscillate less than the corresponding results for the previous case with HAFTOL 10. The CPU time for the implicit method signi®cantly increases, from 106 to 524 s. The periods of the FEM deviate from that of the exact solution. As mentioned above, a is set to À0X05 for the implicit method to induce arti®cial damping. The numerical dissipation leads to an amplitude decay of 2 px [5], where x is the algorithmic damping ratio. This is illustrated in Figs. 2a and c.
analyze and predict failure. Comparing the stresses by these two methods is helpful for failure diagnosis. A three-dimensional example is needed for the above purpose. This example is shown in Fig. 3. A cylindrical disk is formed by two cylinders. The radius of the bottom cylinder is 1 m, and the radius of the top cylinder is 0.25 m. The bottom cylinder is of 0.5 m height, and the height of the top one is 0.1 m. The distance between the centers of these two cylinders is 0.31 m. There are 548 eightnode brick elements and 810 nodes for the mesh. Elements 173 and 409 are chosen for the purpose of stress comparison. Nodes 311 and 709 are one of the eight nodes for elements 173 and 409, respectively. The material is ductile steel with a Young's modulus of E 200 GPa and density of r 7833kgmÀ3 . Poisson's ration is taken to be n 0 in order that the stresses can be simply calculated by
E d ui
d u j siY j Y 2 d x j d xi
iY j 1Y 3
3.2. Impact of an elastic cylindrical disk against a rigid wall
In the ®rst example, the results for the displacement by both the implicit and explicit method in direction 1 agree very well with the initial velocity in direction 1, which is used to control the movement of the bar. Moreover, the displacements in the other two directions are very small. Actually, this example cannot show the effects of these two methods on the displacements of the other two directions and subsequent effects on the stress distribution. In engineering applications, stress analysis is usually important to
Fig. 4. The maximum principal stress SP3 for static contact.
(7)
114
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
Fig. 5. (a) The contact reaction force RF2; (b) the displacement component U1 of node 311; (c) the displacement component U2 of node 311; (d) the stress component S11 of element 173; (e) the stress component S22 of element 173; (f) the stress component S33 of element 173.
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
Fig. 5. (Continued ).
115
116
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
Fig. 6. (a) The displacement U3 by explicit; (b) the displacement U3 by implicit.
in which siY j are the six independent components of stress, ui are the three components of displacement, and x j are the three coordinates. The maximum principal stresses (SP3) by the static algorithm are shown in Fig. 4. The original gap between the model and the rigid wall is 0.01 m. A static displacement of 0.02 m is applied on the bottom surface of the large cylinder. An initial velocity of 10 m sÀ1 is prescribed. Fig. 5a shows the agreement of the contact reaction force RF2 for the results by the two methods. The displacements U1, U2 and U3 of node 311 are shown in Figs. 5a±c. U3 shows a large difference, but the absolute value of U3 for node 311 is relatively small because node 311 is on the axisymmetric axis. Due to the fact that the displacement of the rigid body in direction 2 is comparatively large relative to the deformation, U2 shows an insigni®cant change, but Figs. 5d±f show apparent differences of stresses of element 173 by these two methods. Although at some region such as node 331, the values of U3 are quite different for the two methods, the whole displacement pattern of U3 is similar, as shown in Figs. 6a and b. On the other hand, the stress components S11, S22 and S33 of another element 409 which connects to the contact surface, show agreement by the two methods in Figs. 7a±c. Displacements have higher accuracy than stresses, thus the three displacements U1, U2 and U3 of node 709 obviously agree. The results are testi®ed in Figs. 6a and b. For stress distribution, the regularity is the same. Figs. 7d and e show that the pattern of stress distribution is similar for most regions, especially at regions where the stresses are compressive and the absolute values are large. This is a compressive model, the absolute values of tensile stresses being much smaller than the compressive stresses. The differences shown in Figs. 7d and e are all in the region of tensile stresses. For this fast impact problem, the total time is usually very short, the computational cost of the explicit begin about one tenth that of the implicit.
3.3. Slow contact between an elastic cylindrical disk and a rigid wall
The displacement for the bottom surface of the elastic cylindrical disk is set to be 0.02 m. The total time is 1 s, which is much larger than in the fast case (0.002 s). Fig. 8a shows the displacement U2 of node 311, whilst Fig. 8b shows the stress component S22 of element 173: differences are apparent. Fig. 8b shows the stress component S22 of element 409. The results show good agreement. Figs. 8c and d compare the maximum principal stress SP3. The regularities are the same as mentioned in the last section. The perfect agreement between Fig. 8c and Figs. 6a and b show that the results by implicit methods are very accurate, where under these conditions inertia effects can be neglected. The stable time step for explicit is 7 Â 10À6 s for the slow case and 8 X5 Â 10À6 s for the fast case. For the slow case, the time increment 139935 is much more than that of the fast case. Thus the CPU cost increases dramatically from 7 min to 1 h 17 min. The CPU cost for the implicit method is 10 min.
4. Closure
Fast and slow linear contact problems have been analyzed by different methods in ABAQUS. For the fast case, the advantages of the explicit method are apparent within the desirable tolerance. The cost of the explicit method is much less than that of the implicit method. Due to numerical damping, amplitude decay is observed for the implicit method. For the slow case, the solutions are more unstable because high frequency numerical noise becomes more important. The numerical damping induced in the implicit method shows its function to remove the noise and keep the results more accurate. Because the explicit method is conditionally stable, the stable time period is much smaller than that of the
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
117
Fig. 7. (a) The stress component S11 of element 409; (b) the stress component S22 of element 409; (c) the stress component S33 of element 409; (d) the maximum principal stress SP3 of the impacting disk by explicit; (e) the maximum principal stress SP3 of the impacting disk by implicit.
118
J.S. Sun et al. / Journal of Materials Processing Technology 105 (2000) 110±118
Fig. 8. (Continued ).
implicit method. If the whole procedure is very long, it takes many time increments for the explicit method. Thus the advantages of the implicit method are obvious for the slow contact problem.
References
Fig. 8. (a) The displacement U2 of node 311; (b) the stress component S22 of element 173; (c) the stress component S22 of element 409; (d) the maximum principal stress (SP3) of the slow contact disk by explicit; (e) the maximum principal stress (SP3) of the slow contact disk by implicit.
[1] ABAQUS User's Examples and Theory Manual, Version 5.7. Hibbitt, Karlsson & Sorensen Inc., 1998. [2] N. Rebelo, J.C. Nagtegaal, L.M. Taylor, Comparison of implicit and explicit ®nite element methods in the simulation of metal forming processes, in: Chenot, Wood, Zienkiewicz (Eds.), Numerical Methods in Industrial Forming Processes, 1992, pp. 99±108. [3] Hibbitt, Karlsson & Sorensen Inc., Application of implicit and explicit ®nite element techniques to metal forming, J. Mater. Process. Technol. 45 (1994) 649±656. [4] H.M. Hilber, T.J.R. Hughes, Collocation, dissipation and `overshoot' for time integration schemes in structural dynamics, Earthquake Eng. Struct. Dyn. 6 (1978) 99±117. [5] T.J.R. Hughes, The Finite Element Method Ð Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 07632, 1987. [6] Z.H. Zhong, Finite Element Procedures for Contact Ð Impact Problems, Oxford University Press, New York, 1993.
