INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 44, 155–176 (1999) AN IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS H.-L. CAO∗ AND M. POTIER-FERRY Laboratoire de Physique et MÃecanique des MatÃeriaux, UMR CNRS 7554, UniversitÃe de Metz, Ile du Saulcy, 57045 Metz Cedex, France ABSTRACT In this work, some techniques are proposed to improve the usual Newton–Raphson Method (NRM) used in the numerical analysis of large strain viscoplastic problems. These techniques, based on the ÿrst-order perturbation technique, allow to deÿne an adaptive step strategy and to improve the trial solution (guessed solution) for the ÿrst iteration at each step. To assess the eciency of the proposed techniques, a number of numerical examples are presented: necking of a circular bar, plane strain bending, and an axisymmetric hydrostatic bulging process. Compared with the classical Newton–Raphson method, the proposed scheme requires less matrix inversions, less time steps and less CPU time. For a typical axisymmetric hydrostatic bulging process (606 degrees of freedom), the proposed method needs 5 times less matrix inversions, 1·7 times less time steps, and 4 times less CPU time. Copyright ? 1999 John Wiley & Sons, Ltd. KEY WORDS: Newton–Raphson method; perturbation technique; viscoplasticity; adaptive step 1. INTRODUCTION In the numerical simulation of metal forming processes, the viscoplastic constitutive laws have been widely employed in analysing large strain processes by the ÿnite element method.1–5 Because of material and geometrical non-linearities, some iterative methods such as the Newton–Raphson Method (NRM) are generally used. This kind of methods has two major diculties: the diculty to choose an optimal step size and the diculty to estimate a good trial solution for the ÿrst iteration at each step. For the optimal step size, some ÿnite element codes (ABAQUS, etc) adopt a so-called automatic step size, which is limited according to some criteria (the maximum admissible strain increment, for example). This technique permits the avoidance of very large steps when non-linearities come essentially from geometric variations or from plastic strain increments. It does not, however, allow to choose an adequate step size when the non-linearity is due to the variation of strain rate. For the trial solution, Chung and Wagoner4 and many other authors1–3 have examined the choice of dierent trial solutions. The major idea is to utilize the solution of the previous step or the solution extrapolated from that of several previous steps. As these techniques are based only on ∗ Correspondence to: H.-L. Cao, LPMM, URA CNRS 1215, UniversitÃe de Metz, Ile de Saulcy, 57045 Metz Cedex, France. E-mail: cao@lpmm.univ-metz.fr CCC 0029–5981/99/020155–22$17.50 Copyright ? 1999 John Wiley & Sons, Ltd. Received 17 October 1997 Revised 20 April 1998 156 H.-L. CAO AND M. POTIER-FERRY the solutions of the previous steps, it is dicult to predict a good trial solution particularly when the strain rate varies signicantly. The perturbation technique, coupled with the discretization method, has been applied to numerical analysis of geometrically non-linear problems.6 – 8 Some signicant advances on the perturbation technique, called incremental perturbation method or Asymptotic Numerical Method (ANM), have been achieved by Yokoo et al.,9 Damil and Potier-Ferry,10 as well as by Cochelin et al.11 The same method has been extended to the large strain viscoplastic problems by Potier-Ferry et al.12 Our objective in this paper is to use the same perturbation technique to improve the classical Newton–Raphson Method. The improvements concern particularly the two principal diculties cited above. Recently, a study carried out by Mole et al.13 employed a rst-order perturbation technique to solve large strain viscoplastic problems, although the terminology ‘perturbation technique‘ has not been mentioned in their paper. Their idea is to estimate a solution based on the acceleration (the rst time derivative of velocity) at each step without any iteration. Unfortunately, the given numerical examples were restricted to strain rate sensitive material. The investigation in the present work follows the same principle as Mole et al.13 but diers from theirs in four aspects. Firstly, the proposed improving techniques are deduced in the framework of the perturbation technique and can be extended straightforward to higher orders. Secondly, the proposed computational procedure deals not only with strain rate sensitive materials but also applicable to work hardening materials. Thirdly, the predicted solution, based on the rst time derivative of the velocity, is just taken as the trial solution for the rst iteration at each step. This ensures some accuracy and numerical stability. Finally, on the basis of both strain increment and strain rate increment, an adaptive step strategy is considered to optimize the computational procedure. In addition, the proposed improvements were readily inserted into an existing viscoplastic nite element code. The paper is arranged in the following way. In Section 2, we will give the mechanical model used in this study. In Section 3, a brief review on the classical NRM will be given. In Section 4, we will present the perturbation technique in the framework of large strain viscoplasticity. In Section 5, some techniques will be provided to dene an adaptive step size and to improve the trial solution for the rst iteration at each step. A new computational procedure will be detailed at the end of this section. To assess the eciency of the proposed techniques and check their validity, a number of numerical examples are examined in Section 6 such as necking of a circular bar, plane strain bending, and axisymmetric hydrostatic bulging processes. The computational results obtained by the proposed scheme are compared with that of the classical NRM using equal or adaptive steps, with the numerical solution of the ANM, also with a published computational solution14 as well as with the available experimental results.15 Finally, in Appendix I, an usual computational procedure used in the classical NRM will be given to facilitate the comparison with respect to the proposed method. For the same reason, The ANM will be briey described in Appendix II. Appendix III gives in detail the description to deduce the rst time derivative formulation in large strain viscoplasticity. 2. THE MECHANICAL MODEL TO BE SOLVED This section describes the governing equations of large strain viscoplastic problems and some notations used in the paper. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 157 2.1. Equilibrium equations and boundary conditions By neglecting the inertial and volumic forces, the equilibrium of a solid body occupying a deformed region with boundary @ can be expressed as div b = 0 in ; V = Vd on @ v ; b · n = td on @ (1) where @ v ∪ @ = @ and @ v ∩ @ = ∅. In (1), V and b denote the velocity and stress elds over . V d represents the prescribed velocity over the kinematic boundary @ v and t d the prescribed traction over @ with as a load parameter, which depends on time in a process. For the present work, is taken as time without loss of generality. Finally, n gives the external normal vector to the surface @ . 2.2. The viscoplastic law In the present work, the elastic strain is neglected and a viscoplastic law is assumed: b = b 0 − pI with b 0 = cD 0 tr D = 0 (2) in which b denotes the stress tensor, p the hydrostatic pressure, c the material consistency coecient, and D the strain rate tensor which is the symmetric part of the velocity gradient tensor. The prime 0 and the superscript t imply a deviatoric tensor and a transposed tensor, respectively. The symbol tr(·) gives tr D = Dii , which denes the volume strain rate. The evolution of c is generally governed by r Z t 2 0 0 n m−1 2 D : D = 0 + with D = D ds (3) c = 3 k( + ) D 3 t0 Here k, , n and m are material constants, D is the equivalent strain rate, is the cumulated plastic strain and 0 denes the initial value of the plastic equivalent strain at the beginning of time step t0 : The symbol (:) implies the contracted product of two tensors D : D = Dij Dij . As in most numerical codes, the material incompressibility is treated by using a penality method. In this case, (2)3 can be replaced by tr D = − p g (4) with g a very large constant. To avoid the division by a very small value when the strain rate D is near to zero, equations (3) are regularized as r Z t 2 0 0 n m−1 2 DA ds (5) D : D + Á 2 ; = 0 + c = 3 k( + ) DA ; DA = 3 t0 in which Á is a small numerical parameter. If Á is zero, equations (5) reduce to (3). In most viscoplastic nite element codes, a tolerance is introduced to prevent the division by a very small value. These two techniques have not much essential dierence. From now on, equations (4) and (5) will be employed in place of (2)3 and (3). Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 158 H.-L. CAO AND M. POTIER-FERRY 2.3. Principle of virtual velocities In practice, a spatial discretization technique such as the nite element method is commonly used to obtain approximate solutions to the above equations. To do that, it is necessary to reformulate equations (1) and (2) into a weak form—the principle of virtual velocities. They are thus rewritten as Z Z (b 0 : D 0 − p tr D) d = t d · V d @ Z p p d = 0 (6) tr D + g for all weight functions p and V which satisfy the homogenous boundary conditions on @ v . In (6), represents boundary surfaces. The unknowns are the velocity V and the pressure p elds. The integration is carried out on the deformed conguration estimated at the end of the time step. At a given time step, the Newton–Raphson method uses the linearized form of (6) to iterate until its residual norm is admissible at the end of the time step. In order to compute the time derivatives of the velocity and the pressure at each time step, the above equations should be expressed on a known reference conguration (the conguration at the beginning of the time step for instance), designated as 0 with boundaries @ 0 , in the form Z Z Z p p d = 0 (7) T : L d = t d · V d ; J tr D + g 0 @ 0 0 with T being the rst Piola–Kirchho stress tensor, and t d being the traction vector, T = J bF−t or T = J (b 0 − pI)F−t ; t d = J bF−t N (8) Here, F is the deformation gradient tensor, L the velocity gradient tensor with respect to the reference conguration, J the determinant of F and N indicates the external normal vector to the surface @ 0 . In the cases where the stress is prescribed b = b d instead of t d on the surface, t d varies as a function of F (in an hydrostatic bulging process, for instance). As in most viscoplastic nite element codes, the process here is described using an updated Lagrangian scheme. The co-ordinate vector x of the deformed conguration is updated at the end of each time step as Z t V ds (9) x = x0 + t0 where x0 is the co-ordinate vector at the beginning of time step t0 : This new updated conguration is taken as the reference conguration for the subsequent time step. 3. REVIEW ON THE CLASSICAL NEWTON–RAPHSON METHOD The Newton–Raphson method is well known. A general review on the method can be found in Reference 16. Here, we make precise only the NRM for large strain viscoplastic problems using the updated Lagrangian scheme, which allows us to understand more clearly the proposed improving techniques and the facilities of its insertion into an existing viscoplastic nite element code. An usual computational procedure used in the classical NRM is given in Appendix I. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 159 Generally, the classical NRM uses the linearized form of (6) at each iteration.1 This form can be expressed as Z [(M c : dD 0 ) : D 0 − dp tr D] d = R ext − Rint Z dp p d = Rp tr dD + (10) g with nDA t 2c D0 ⊗ D0 (m − 1·0) + M = cI I + ( + ) 3DA2 Z 1 dD = (dl + dlt ); Rint = (b 0 : D 0 − p tr D) d 2 Z Z p p d tr D + t d · V d ; Rp = − R ext = g @ c where M c is the eective viscoplastic modulus tensor (symmetric). dl, dD 0 , and dp are the corrections between two successive iterations to the eld velocity gradient l, to the eld D 0 and to the eld p, respectively. t represents the time increment. The integration and its integral function in (10) are evaluated by using the newly obtained solution of the preceding iteration. Finally, the 0 . symbols and ⊗ imply: (I I)ijkl = ik jl and (D 0 ⊗ D 0 )ijkl = Dij0 Dkl By using the nite element discretization, we obtain a set of equations from (10), schematically as [K c ]{dU } = {R} (11) where {dU } includes the correction of velocity {dV } and the correction of pressure {dp}. [K c ] stands for the tangent stiness matrix (symmetric) and {R} denotes the residual force array. This is the system to be solved at each iteration in the classical NRM. 4. PERTURBATION TECHNIQUE Let us represent the weak form (7) by the symbol (U; ) to simplify the discussion. Here U stands for both the unknown velocity and the unknown pressure, and denotes the load parameter. Since U and are functions of time, (U; ) can also be considered as a function of time. Supposing that at time t0 we have U(0) and (0) which satisfy the weak form (7): (U(0) ; (0) ) = 0 (12) we can take the Taylor expansion of U, and (U; ) in the vicinity of t0 for any time varying from t0 to t: U() = U(0) + U; t |t0 ( − t0 ) + () = (0) + ; t |t0 ( − t0 ) + (U(); ()) = 0 Copyright ? 1999 John Wiley & Sons, Ltd. 1 2! U; tt |t0 ( 1 2! ; tt |t0 ( − t0 ) 2 + · · · − t0 ) 2 + · · · (13) Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 160 with H.-L. CAO AND M. POTIER-FERRY (U(); ()) = (U(0) ; (0) ) + ; t |t0 ( − t0 ) + 1 2! ; tt |t0 ( − t0 ) 2 + · · · where the symbol (•); t |t0 denotes the rst time derivative taken at time t0 and (•); tt |t0 is the second time derivative. From now on, we will use the symbol (•)(1) to designate (•); t |t0 and (•)(2) for (•); tt |t0 : As relation (13)3 should be valid whatever is the value , we obtain, by taking into account (12): (1) = 0; (2) = 0; : : : (14) Futhermore, (14) can be reformulated as (1) = ; U |t0 U(1) + ; |t0 (1) = 0 (2) = ; U |t0 U(2) + ; |t0 (2) + ; UU |t0 U(1) U(1) + 2; U |t0 U(1) (1) + ; |t0 (1) (1) = 0 ::: (15) From (15)2 , we observe that the second time derivatives U(2) and (2) have the same operator as the rst time derivatives U(1) and (1) : The remaining terms involve only terms dependent on lower time derivatives (the rst time derivatives in (15)2 , for instance). For the higher time derivatives U(3) , (3) ; : : : ; the same characteristic can be observed. It is this feature which is exploited in the ANM12; 18 because it needs only one matrix inversion to estimate any higher time derivatives. Neglecting the higher time derivatives than the one in (13) and taking into account (15), we have U() = U(0) + U(1) ( − t0 ); () = (0) + (1) ( − t0 ) ; U |t0 U(1) + ; |t0 (1) = 0 (16) Of course, we can keep the higher derivatives until n order in the computation. In this case, the method is called as asymptotic numerical method10 or incremental perturbation method.9 The interested reader can nd more details for large strain viscoplastic problems in Reference 12. A short description is given in Appendix II. In the present work, we follow the same idea as Mole et al.13 by using only the rst time derivative to improve the classical Newton–Raphson method. Introducing the rst time derivative of (7) at time t0 into (16)3 we obtain Z 0∗ ext [(M : D(1) ) : D 0 − p(1) tr D] d = R(1) − Rv ∗ 0 Z 1 ∗ (17) tr D(1) + p(1) p d = Rp ∗ g 0 with D∗ (1)ij ext R(1) @V(1) j 1 @V(1)i 2c(m − 1) 0 = + D ⊗ D0 ; M = cI I + 2 @xj @xi 3DA2 Z d d = ((1) t d + t(1) ) · V d ; t(1) = 0; if t d is prescribed @ 0 d = b d (tr (l)I − lt )N; t(1) Copyright ? 1999 John Wiley & Sons, Ltd. if b = b d is prescribed Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS R v∗ Z = 0 161 {b[tr(l)I − lt ] : l + [(Rc)D 0 − M : RD ] : D 0 } d (Rc) = 23 kn( + 0 ) (n−1) DA ; RD = 12 (lc + lct ); lc = ll Z p p∗ t R = − tr(l)I + l : [tr(l)I − l ] p d g 0 in which the integration and all the integral functions are estimated on the conguration at time t0 : From this set of equations we obtain a linear system of unknown elds V(1) and p(1) once V(0) and p(0) are known. Some details about the derivation of equations (17) are presented in Appendix III. Comparing system (10) with (17), it can be noted that they have the same structure with their respective unknowns (dV, dp in (10) and V(1) , p(1) in (17)). If system (10) is constructed by using the last convergent solution at time t0 , the only dierence between the operators of their unknowns concerns the eective viscoplastic modulus M c and M: This dierence comes from the work hardening of the material. If n = 0, M is identical to M c , which has been dealt with by Mole et al.13 If t = 0, M is also identical to M c : This last point is taken into account in the ANM to nd out the solution at time t0 and the corresponding time derivatives. In general, M c and M are dierent. However, M can be approximated by M c for most applications. In fact, the term involving work hardening in M c can be rewritten as −1 + 0 nDA t =n 1 + + DA t The value of this term is small with regard to |m − 1| for most metallic materials where n is around 0·2 and m is less than 0·5. Particularly, equivalent strain increment is usually controlled under 5 per cent in numerical simulations. If the term involving work hardening is neglected, it can decrease the convergent rate in the NRM. On the contrary, the substitution of M with M c in (17) gives a good approximation of V(1) and p(1) without constructing and inversing the stiness matrix once more. This will be shown by the numerical examples. After nite element discretization, system (17) becomes [K]{U(1) } = {R∗ } (18) where {U(1) } involves the rst time derivatives of velocity {V(1) } and pressure {p(1) } at time t0 : [K] stands for the stiness matrix, which will become [K c ] if M in (17) is replaced by M c : ext , Rv ∗ and Rp ∗ : This residual Finally, {R∗ } denotes the residual force array coming from R(1) force can be calculated in a similar way as {R} in (11). 5. IMPROVED NEWTON–RAPHSON METHOD Three points will be discussed in this section, the adaptive step size, the trial solution for the rst iteration at each step and the consequent computational procedure. 5.1. Adaptive step size As indicated by Chung and Wagoner,4 an adaptive step size is generally better than an optimal equal one in terms of computational eort. This is why the adaptive step size is implemented in some commercial codes (ABAQUS, etc.). The problem for this strategy concerns essentially Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 162 H.-L. CAO AND M. POTIER-FERRY their chosen criteria by which the step is controlled. Normally, the usual criteria are the maximum admissible strain increment and the maximum incremental rotation. These limitations allow to avoid divergence if the non-linearities originate from the variation of geometry or that of the equivalent strain. However, when the non-linearity comes from the variation of strain rates, the previous limitations will lose their eciency. As demonstrated in Section 4, the rst time derivatives of velocity and pressure may be computed easily once the velocity and pressure elds are known. We can use these data to control the step size apart from the classical limitations on strain increment and incremental rotation. Based on both the limitations of velocity variation and that of the velocity, the adaptive step size becomes more ecient even if the strain rate varies signicantly. In the realized nite element code, the following criteria are chosen: limitation on the increment of velocity gradient |l| and limitation on the relative variation of velocity gradient with respect to its velocity gradient |l(1) |=|l|. Of course, other limitations can also be implemented into the code, limitations based on strain increment and the variation of strain rate, for example. All these will not inuence the adaptive step size considerably. Compared with the method proposed by Mole et al.,13 the introduction of the present adaptive step size is more satisfactory, which can be observed from the numerical examples presented in Section 6. 5.2. Trial solution for the rst iteration at each step If V(0) and p(0) at time t0 are obtained, we can use the last inversed stiness matrix [K c ]−1 instead of [K]−1 to compute V(1) and p(1) for system (18). The next time step size t can be chosen according to the criteria discussed above. After the choice of step size, the trial solution for the rst iteration at time t0 + t can be predicted as follows: V = V(0) + V(1) t; p = p(0) + p(1) t (19) As a consequence, the trial conguration x and the trial equivalent strain at the time t0 + t can be estimated as x = x0 + tV(0) + 12 t 2 V(1) ; = 0 + tDA(0) + 12 t 2 DA(1) (20) To keep the consistancy with (20), the conguration x and equivalent strain at the end of each iteration can be updated as x = x0 + 12 t(V(0) + V); = 0 + 12 t(DA(0) + DA ) (21) where V is the newly obtained solution at the time t0 + t, and DA is the equivalent strain rate corresponding to V and x: Furthermore, the viscoplastic modulus M c are modied to take into account the updated scheme (21), which is rewritten as nDA t 2c D0 ⊗ D0 (m − 1) + (22) M c = cI I + 2( + ) 3DA2 It is this modulus which is used to construct the stiness matrix [K c ] in the improved Newton– Raphson method. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 163 5.3. Computational procedure With the consideration of the two techniques above, the new computational procedure is conceived as follows for the cases where the boundary conditions are supposed to be unchanged. For the rst time step, the classical NRM is employed to nd the solution. Some measures proposed by Wang1 and others2; 4; 5 can be used to increase the convergence rate for this rst time step. In this study, we focus on the improved techniques only for the numerical solutions after the rst step. (0) Given a prescribed time step t, the initial time t0 , the initial reference conguration, the initial equivalent plastic strain eld and the trial solution, nd out the solution {U } at time t (1) = t0 + t by using the classical NRM. Loop of the deformation process for l = 2 to lmax (1) Construct the residual force array {R}∗ of system (18). (2) Compute {U(1) } by using the last inversed stiness matrix [K c ]−1 at the time step t (l−1) : (3) Choose the time step size t a according to some criteria (the two criteria described above, for example). (4) With the chosen step size, estimate the trial solution by (19), the trial conguration and the trial equivalent strain by (20) at time t (l−1) + t a : Loop of the equilibrium at each step for n = 1 to n max (5) Estimate the residual forces {R} of (11). (6) Check if the norm of the residual forces {R} is less than the admissible tolerance. If yes, pass to (12); if not, continue. (7) Form the tangent stiness matrix [K c ] of (11) with consideration of the modied viscoplastic modulus M c in (22). (8) Solve system (11) for the correction elds {dU }. (9) Compute the new solution {U } (n) = {U } (n−1) + ÿ{dU } (23) where ÿ is a relaxation coecient varying from 0 to 1. (10) Update the reference conguration and the equivalent strain eld according to (21) with the determined time step size t a : (11) end of the equilibrium loop. (12) Update the time: t (l) = t (l−1) + t a (13) Verify if the time is superior or equal to the nal time of the process. If yes, pass to (15); If not, continue. (14) End of the deformation process loop. (15) End of the process. In the above procedure lmax denotes the maximum number of time steps and n max is the maximum number of iterations. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 164 H.-L. CAO AND M. POTIER-FERRY As compared to the classical NRM with adaptive step size presented in Section 3, the proposed computational procedure has three major features. Firstly, the adaptive step size is determined according to the criteria based not only on the current velocity and pressure but also on their rst time derivatives. Secondly, the trial solution for the rst iteration at each step is calculated according to (19) in place of the solution at the preceding time step. Finally, the eective viscoplastic modulus is modied to take into account the updated scheme expressed in (22). Between the proposed computational procedure and that described in Section 3, the dierent points are concerned with the steps about the choice of the trial solution and that of the step size. The rest of the procedure is the same. This characteristic allows one to insert easily the present scheme in an existing viscoplastic nite element code. The scheme described by Mole et al.,13 which is restricted to strain rate sensitive material without work hardening, adopts the equal step sizes and considers the trial solution by (19) as the real solution for the time t0 + t: Their scheme needs only one iteration per step after the rst step of the deformation process. To ensure a good solution, however, the step size should be suciently small. If not, the obtained solution may become completely erroneous. Comparing their scheme with the proposed one, the latter takes into account the work hardening eect and adopts an improved adaptive step strategy. In addition, the former can be considered as a particular case of the latter. It is sucient to relax all limitations on the step size and put n max equal to 1. In comparison with the ANM (a short description is given in Appendix II), the present scheme needs only the rst time derivatives while the ANM requires generally until 10th or higher time derivatives. The proposed scheme thus reduces the memory storage. Since the present scheme needs only the rst time derivatives, it is not necessary to regularize strongly the used formulations (presenting a sharp variation) as recommended in Reference 17 for the ANM. Moreover, it is easier to implement the present scheme into an existing code than the ANM. 6. NUMERICAL EXAMPLES In this section, several numerical examples are presented to assess the performance of the improving techniques proposed above. For each example, ve solution procedures are employed: the classical NRM with equal time steps, the classical NRM with adaptive time steps, the one-iteration scheme using equal time steps (it becomes the scheme proposed by the Mole et al.13 in the case of strain rate sensitive material), the ANM until the 15th time derivative, and the present proposed scheme. The numerical results obtained by these ve dierent procedures are compared in terms of their computational requirements (the number of time steps and that of matrix inversions) and their nal kinematic variables (velocity and equivalent strain). Three kinds of examples are given here. The rst is concerned with the necking of a circular bar, which has been considered by Mole et al.13 The second is plane strain bending under hydrostatic pressure. The nal example deals with an axisymmetric hydrostatic bulging process, simulated by Kim and Yang14 and realized experimentally by Iseki et al.15 In the simulations, a 4-node isoparametric element with bilinear velocity interpolation and constant pressure is used. A relative residual norm of 10−2 is generally adopted. For the simulations except that using the ANM, ÿ is taken equal to 0·75 (described in (23)) to facilitate the comparison. As information, all the simulations have been carried out on a HP9000=K200 server. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 165 Figure 1. (a) The geometry of the circular bar (mm); (b) the initial mesh; (c) the deformed mesh corresponding to the nal vertical displacement of the upper head 7·1 mm 6.1. Necking of a circular bar In this example, the geometry, the initial nite element mesh, and the deformed mesh for a half of a cross-section of the specimen are shown in Figure 1. The specimen is assumed to be pulled by an uniform velocity of 1000·0 mm=s at the upper head and is xed in the radial direction. One hundred elements are used in the computation. The rst simulation is carried out with k = 1·0 Pa s−0·1 , m = 0·1, and n = 0 as the material constants, which have been considered by Mole et al.13 For Á, it is taken equal to 10−5 =s. The rst step size in the simulation is set to 10−4 s. This step is solved by using the classical NRM. After the rst step, we apply the ve computational procedures for the rest of the deformation process. In Table I are registered the equivalent strain at the centre of the specimen, the radial velocity of the external node of the central section, and the computational requirements except the rst step. All the results correspond to the nal vertical displacement of the upper head 7·1 mm. The reference results are obtained by using the classical NRM with equal time steps t = 10−5 s. From Table I, we observe that the NRM with adaptive steps needs less matrix inversions and time steps than the NRM with equal steps to reach the same accuracy in terms of : This remark conrms the observation given by Chung and Wagoner.4 The most ecient case of Mole’s scheme and that of the proposed scheme require generally less matrix inversions (2 times less) than the optimal NRM. Moreover, the present scheme obtains a very good solution (inferior to 1·0 per cent Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 166 H.-L. CAO AND M. POTIER-FERRY Table I. The equivalent strain at the centre, the radial velocity of the external node at the central section, and the computational requirements except the rst step (time steps, matrix inversions simplied as MI), obtained when the tensile process is stopped at 7·1 mm vertical displacement of the upper head. Here Vref = −985·7 mm=s and ref = 1·4180 Steps MI V=Vref =ref NRM equal step t (10−3 s) 0·1 0·25 0·5 71 28 14 138 56 28 1·0149 1·0016 1·0058 1·0045 1·0133 1·0333 NRM adaptive step |l|max and tmax = 0·75 (10−3 s) 0·05 0·1 43 21 55 37 1·0321 1·0610 0·9922 0·9925 Mole’s scheme t (10−3 s) 0·1 0·25 0·5 71 29 15 71 29 15 1·0000 1·0007 1·0021 0·9950 0·9941 0·9937 Present scheme |l|max and (|l(1) |=|l|)max 0·01 0·1 0·1 0·1 0·15 0·1 223 223 1·0000 0·9950 23 23 1·0001 0·9945 15 15 1·0005 0·9937 10 10 0·9991 0·9946 ANM (time derivatives) 15 dierence) with respect to the reference. Finally, the ANM, as expected, gives an accurate solution with less matrix inversions and time steps with respect to all other schemes. Another simulation is carried out with a material having both, work hardening eect and strain rate sensitivity. The material constants are: k = 1·0 Pa s−0·01 , m = 0·01, n = 0·20, and = 0·001. In Table II are registered the same kind of data as in Table I. Some similar remarks can be observed as in the preceding example. On the contrary, The one-iteration scheme, whenever the matrix [K c ] or [K] in (18) is used, gives similar results. It generally requires less matrix inversions and time steps than the optimal NRM. However, the obtained velocity is completely erroneous with respect to the reference in the cases of t = 0·125 × 10−3 s and t = 0·25 × 10−3 s. This gives us some sign of numerical instability of the scheme when the time step is a little bigger. 6.2. Plane strain bending In this example, a metallic blank of 24 mm in width and 0·314 mm in thickness is subjected to a time-dependent hydrostatic pressure. The material constants are taken as: k = 180·0 MPa s−0·2 , m = 0·0, n = 0·29, = 0·000769, and Á = 10−5 =s. Because of the symmetry, only a half of the blank is used in the computation. Concerning the boundary conditions, the blank is clamped at its border and xed at the centre in the horizontal direction. On the upper surface, a hydrostatic pressure is applied with a rate of 10−3 MPa=s. The blank is discretized by means of two elements across the thickness and 20 elements in horizontal direction. In Figure 2 the deformed mesh is displayed corresponding to t = 1353·3 s after the start of the deformation process. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 167 IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS Table II. The equivalent strain at the centre, the radial velocity of the external node at the central section, and the computational requirements except the rst step (time steps, matrix inversions simplied as MI), obtained when the tensile process is stopped at 7·1 mm vertical displacement of the upper head. Here Vref = −1134·1 mm=s and ref = 1·3480 Steps MI V=Vref =ref NRM equal step t (10−3 s) 0·1 0·125 0·25 70 56 100 1·0411 104 1·0503 No convergence 0·9874 0·9887 NRM adaptive step |l|max and tmax = 0·25 (10−3 s) 0·025 0·04 0·05 79 49 113 1·0259 84 1·0361 No convergence 0·9874 0·9867 One-iteration scheme with [K c ] t (10−3 s) 0·1 0·125 0·25 0·35 70 56 28 70 0·9986 56 0·7555 28 0·4646 Erroneous results 0·9941 0·9941 0·9503 One-iteration scheme with [K] t (10−3 s) 0·125 0·25 0·35 56 28 56 0·7417 28 0·3154 Erroneous results 0·9941 0·9573 Present scheme |l|max and (|l(1) |=|l|)max 0·01 0·1 0·05 0·1 0·075 0·1 205 206 0·9987 0·9941 42 46 0·9993 0·9934 ANM (time derivatives) 15 No convergence 39 40 0·9989 0·9941 Figure 2. The deformed mesh corresponding to time t = 1353·3 s in the plane strain bending process Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 168 H.-L. CAO AND M. POTIER-FERRY Figure 3. The variation of vertical velocity at the centre as a function of time Table III. The equivalent strain at the centre and the computational requirements (time steps, matrix inversions simplied as MI), obtained when the plane strain bending process is stopped at t = 1353·3 s. Here ref = 0·2562 Steps NRM equal step t (s) 1·0 5·0 10·0 12·5 MI =ref 1354 272 137 3865 0·2564 1230 0·2569 852 0·2568 No convergence 1·0008 1·0027 1·0023 NRM adaptive step |l|max and tmax = 25·0 s 0·01 0·02 0·03 0·04 225 112 74 1097 0·2566 656 0·2569 796 0·2571 No convergence 1·0016 1·0027 1·0035 One-iteration scheme with [K c ] t (s) 1·25 2·5 5·0 7·5 1059 530 266 1161 0·2559 632 0·2560 368 0·2561 Erroneous results 0·9988 0·9992 0·9996 One-iteration scheme with [K] t (s) 0·5 1·0 1·25 2646 2768 0·2559 Erroneous results Erroneous results 0·9987 Present scheme |l|max and (|l(1) |=|l|)max 0·01 0·1 0·02 0·1 ANM (time derivatives) 15 275 340 0·2561 0·9996 162 308 0·2562 1·0000 243 311 0·2561 0·9996 For the rst time step, t is set to 2·5 s to be able to attain rapidly the convergence (three iterations). With this rst time step size, the one-iteration scheme has not obtained the right solution for t superior to 0·5 s. In order to be able to increase the time step size, a bigger time step size t = 30·0 s is used at the rst step for the simulations using the one-iteration scheme. The reason Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 169 Figure 4. The deformed mesh corresponding to time t = 2090·0 s in the axisymmetric hydrostatic bulging process for this problem can be explained by the sharp time variation of the central vertical velocity at the beginning of the process shown in Figure 3. From this gure, it is observed that it is very dicult to estimate the real solution without any iteration near to t = 15·0 s if the step is too big. Table III shows the equivalent strain at the centre and the computational requirements, corresponding to the time t = 1353·3 s. It is recalled that the simulations using the one-iteration scheme adopts a time increment of 30·0 s for the rst step instead of 2·5 s as used in other simulations. The reference results are obtained by using the classical NRM with equal steps of t = 5·0 × 10·0−2 s. From this table, we observe the same phenomena as in the preceding examples. We emphasize here that the one-iteration scheme using [K c ] shows much better performance than that of using [K]: Moreover, the present scheme illustrates well its eciency with respect to others even compared to the ANM. 6.3. Axisymmetric hydrostatic bulging For this example, a metallic blank of 24 mm in diameter and 0·314 mm in thickness is subjected to a time-dependent hydrostatic pressure. The material constants are the same as in the preceding example of the plane strain bending. This example has been simulated by Kim and Yang14 using the classical NRM and experimentally realized by Iseki et al.15 In the simulations, the similar boundary conditions and the same mesh topology are applied as for the plane strain bending. In Figure 4 the deformed mesh is displayed corresponding to t = 2090 s after the start of the deformation process. For the rst time step, t is set to 2·5 s. Table IV shows the equivalent strain at the centre and the computational requirements, corresponding to t = 2090·0 s. The reference results are obtained by using the classical NRM with equal steps of t = 0·1 s. We note here that the most ecient case of the proposed scheme requires 2 times less matrix inversions and almost 3 times less time steps than the optimal case of the oneiteration scheme. With respect to the best of the NRM, the gain is 1·5 times in terms of matrix inversions. Regarding the ANM, it still stays the most satisfactory. Figure 5 depicts the variation of the time steps during the deformation process, obtained by the proposed scheme with |l|max = 0·01 and the maximum of (|l(1) |=|l|) equals to 0·1. It is remarked that the step varies signicantly during the process. In Figure 6 the variation of polar thickness Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 170 H.-L. CAO AND M. POTIER-FERRY Table IV. The equivalent strain at the centre and the computational requirements (time steps, matrix inversions simplied as MI), obtained when the axisymmetric hydrostatic bulging process is stopped at t = 2090·0 s. Here ref = 0·5278 Steps NRM equal step t (s) 6·25 12·5 25·0 37·5 336 169 85 NRM adaptive step |l|max and tmax = 25·0 s 0·01 0·02 0·03 163 90 85 One-iteration scheme with [K c ] t (s) 1·0 2·5 5·0 7·5 One-iteration scheme with [K] t (s) 5·0 6·25 7·5 Present scheme |l|max and (|l(1) |=|l|)max 0·01 0·1 0·015 0·1 0·02 0·1 ANM (time derivatives) 15 MI 1313 0·5258 776 0·5235 734 0·5161 No convergence 729 566 788 =ref 0·9962 0·9919 0·9778 0·5237 0·5181 0·5159 0·9922 0·9816 0·9775 2092 838 420 2093 0·5279 839 0·5279 421 0·5278 Erroneous results 1·0081 1·0002 1·0000 420 421 0·5277 Erroneous results Erroneous results 0·9998 206 247 0·5278 1·0000 157 224 0·5278 1·0000 137 421 0·5265 0·9975 96 97 0·5278 1·0000 Figure 5. The variation of time step size as a function of time strain versus hydrostatic pressure is shown, corresponding to, respectively, the present analysis, the viscoplastic nite element analysis using the NRM14 and the experimental data published in Reference 15. It is seen that the present numerical results agree well with the experimental data15 and the published computational results.14 Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 171 IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS Figure 6. The variation of polar thickness strain as a function of the hydrostatic pressure Table V. The equivalent strain at the centre and the computational requirements (time steps, matrix inversions simplied as MI and CPU time), obtained when the axisymmetric hydrostatic bulging process is stopped at t = 2000·0 s. Here the blank is discretized by two elements across the thickness and 100 elements in the radial direction Steps NRM equal step t (s) 2·0 3·0 6·25 NRM adaptive step |l|max and tmax = 25·0 s 0·0025 0·01 One-iteration scheme with [K c ] t (s) Present scheme |l|max and (|l(1) |=|l|)max ANM (time derivatives) MI CPU time (h) 1000 3055 0·4191 28·19 No convergence at t = 1931·5 s No convergence at t = 1558·75 s 695 2288 0·4196 21·56 No convergence at t = 1617·8 s 1·0 1·5 2·5 1999 2001 0·4193 22·79 Erroneous results from t = 1990·0 s Erroneous results from t = 1895·0 s 0·01 0·1 598 617 0·4193 6·81 365 652 0·4193 8·42 15 In order to check if the proposed scheme still keeps its eciency for problems where the number of degrees of freedom is relatively important, the preceding axisymmetric hydrostatic bulging process is simulated with a ner mesh (two elements across the thickness and 100 elements in the radial direction). In Table V the computational requirements (matrix inversions, time steps and CPU time) are registered, corresponding to t = 2000·0 s. It is observed clearly that the proposed scheme shows its superiority with respect to any other scheme in terms of matrix inversions, time steps and also CPU time. In addition, the one-iteration scheme also gives good results with respect to the classical NRM using equal steps, 1·5 times less in terms of matrix inversions and 1·24 times less in terms of CPU time. This result does not agree well with the observation given by Mole et al.13 where Mole’s scheme required 3 times more CPU time than the classical NRM. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 172 H.-L. CAO AND M. POTIER-FERRY 7. CONCLUSIONS In the present work, an improved Newton–Raphson scheme is presented for numerical analysis of large strain viscoplastic problems. The improvements, based on the perturbation technique, are concerned with the adaptive step strategy and with the rst trial solution at each step. From all the performed numerical simulations, some remarks can be made here. (1) The proposed scheme generally obtains an accurate solution with less computational requirements with respect to the classical NRM and also with respect to the one-iteration scheme (Mole’s scheme in the cases of strain rate sensitive materials). For a typical axisymmetric hydrostatic bulging process, the proposed scheme requires 5 times less matrix inversions, 1·7 times less time steps, and 4 times less CPU time with respect to the classical NRM using equal steps; (2) it is always easier for users to prescribe some limitations with physical signication, which are independent of the treated problems, than to look for the relative optimal time step used in the NRM by doing numerous tests; (3) although the ANM seems in general better than the proposed scheme in this work, it needs also more memory storage and more modications to be inserted in an existing viscoplastic nite element code. It is useful to note here that the proposed improvements are given in the framework of large strain viscoplastic problems but the basic idea is amenable to other non-linear problems. Additionally, as a restriction of the proposed improving techniques, a direct solver should be used to solve the linear system formed by nite element discretization. If an iterative solver is used, the proposed scheme will lose its advantages. APPENDIX I. AN USUAL COMPUTATIONAL PROCEDURE USED IN THE CLASSICAL NRM In a viscoplastic nite element code, the computational procedure can be summarized as follows if the boundary conditions are supposed to be unchanged. Known parameters are: the prescribed time step t, the initial reference conguration and the initial equivalent plastic strain eld. Loop of the deformation process for l = 1 to lmax (0) The trial solution is given and t (0) = t0 Loop of the equilibrium at each step for n = 1 to n max (1) Estimate the residual forces. (2) Check if the norm of the residual forces {R} in (11) is less than the admissible tolerance. If yes, pass to (9); if not, continue. (3) Construct the tangent stiness matrix [K c ]. (4) Solve the equation system (11) to obtain the correction elds {dU }. (5) Compute the new solution: {U } (n) = {U } (n−1) + ÿ{dU } where ÿ is a relaxation coecient varying from 0 to 1. (6) Choose the time step size t c according to the used criteria (the maximum admissible equivalent strain increment, for example). Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 173 (7) Update the reference conguration and the equivalent strain eld: {X } (n) = {X } (n−1) + t c {V } n {} (n) = {} (n−1) + t c {D} n (8) End of the equilibrium loop. (9) Update the time: t (l) = t (l−1) + t c : (10) Verify if the time is superior or equal to the nal time of the process. If yes, pass to (12); if not, continue. (11) End of the deformation process loop. (12) End of the process. For step (0), the trial solution can be taken as the solution of the previous step except for the rst step, where the trial solution may be the solution of the linear Newtonian problem (m = 1·0 and n = 0·0), or that of an elastic problem. In step (5), the value ÿ has an important inuence on the convergence rate. Wang1 has proposed some variation of ÿ for several typical forming processes. If step (6) is removed from the above procedure, we obtain the usual procedure used in the NRM with equal steps. APPENDIX II. AN ASYMPTOTIC NUMERICAL METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS The ANM, based on the perturbation technique presented in Section 4, consists in expanding the unknowns (velocity V and pressure p) into power series (13) and in solving the obtained linear equations for each time derivative in a recursive manner. As all the linear systems at a time step t0 have the same stiness matrix, only one matrix inversion is necessary to estimate each time derivative. With the time derivatives available, it is trivial to compute the velocity and pressure for any time according to their respective power series. As the power series have their own convergence radius, the estimated solution can be accepted when the time increment t from t0 is inside the radius. A bisection method18 is used to determine this convergence radius or the time interval during which the relative residual norm and the volumic strain rate are inferior to the prescribed tolerances. Since the time step is estimated posteriorly, there is no divergence problem as it is often encountered in the NRM. With the valid time increment t, it is easy to compute the velocity and pressure as well as other variables at time t0 + t: Following the updated Lagrangian scheme, the conguration and any used variables at time t0 + t are considered as a new starting point on the equilibrium path. As the solution estimated by the power series at time t0 + t is not always the ‘just’ solution after the updated conguration and the change of the starting point for the expansion, an iterative Newton–Raphson loop is introduced to improve the solution but without any time increment. After this, the same procedure is repeated to compute the time derivatives at time t0 + t and so on. The computational procedure for the ANM is briey described as follows: (0) Given a prescribed time step t, the initial time t0 , the initial reference conguration, the initial equivalent plastic strain eld and the trial solution, nd out the solution {U } at time t (1) = t0 + t by using the classical NRM. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 174 H.-L. CAO AND M. POTIER-FERRY Loop of the deformation process for l = 2 to lmax Loop of the adjustment at time t (l−1) for n = 1 to n max without time increment (1) Estimate the residual forces {R} of (11) using the last obtained solution. (2) Check if the norm of the residual forces {R} is less than the admissible tolerance. If yes, pass to (7); if not, continue. Here, the criteria are stringent in terms of relative residual norm and volumic strain rate. (3) Form the tangent stiness matrix [K] of system (18). (4) Solve system (11) for the correction elds {dU }. (5) Compute the new solution {U } (n) = {U } (n−1) + ÿ{dU } where ÿ is a relaxation coecient varying from 0 to 1. (6) End of the adjustment loop. Loop of the time derivatives for r = 1 to rmax (7) Form the residual force array such as {R}∗ of system (18), but corresponding to the time derivative at order r: (8) Compute {U(r) } by using the last inversed stiness matrix [K]−1 . (9) End of the time derivative loop. (10) Determine the time interval t a , during which the criteria in terms of relative residual norm and volumic strain rate are satised by using a bisection method.18 The criteria are less stringent than that in (2). (11) Estimate the solution of {p} and {V } according to their respective power series at time t (l−1) + t a : Other variables can be evaluated according to the corresponding relations regarding to {p} and {V }: (12) Update the time: t (l) = t (l−1) + t a (13) Update the reference conguration and the corresponding variables at time t (l) . (14) Verify if the time is superior or equal to the nal time of the process. If yes, pass to (16); If not, continue. (15) End of the deformation process loop. (16) End of the process. APPENDIX III. THE SET OF EQUATIONS FOR THE FIRST TIME DERIVATIVE FIELDS V(1) AND p(1) Equations (17) with unknown elds V(1) and p(1) can be deduced by putting equal to zero the rst time derivative of equations (7) at time t0 , which comes from the condition described by equation (16)3 . Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) IMPROVED ITERATIVE METHOD FOR LARGE STRAIN VISCOPLASTIC PROBLEMS 175 From (7), we can express readly their rst time derivatives at time t0 and identify them into zero as follows: Z Z d T(1) : L d = ( (1) t d + t(1) ) · V d Z 0 @ 0 p(1) p + J tr D(1) + p d = 0 J(1) tr D + g g 0 (24) where the symbol (•)(1) denotes the rst time derivative of the corresponding variable at time t0 . d and the known elds The main task is to establish the relations between T(1) , tr D(1) , J(1) , t(1) V and p at time t0 as well as their rst time derivative elds V(1) and p(1) at the same time. From the relation D = 12 (ḞF−1 + F−t Ḟ t ), we can deduce its rst time derivative: ∗ −R D(1) = D(1) D with RD = 12 (lc + lct ) (25) in which D∗(1) and lc are dened in (17). With (25) available, it becomes trivial to have the rst time derivative of tr D and that of D 0 , ∗ − tr R ; D 0 = D 0 ∗ − R 0 (26) tr D(1) = tr D(1) (1) D D (1) From the relations b 0 = cD 0 and c = 23 k( + ) n D m−1 , we can express the rst time derivative of b 0 as 2c(m − 1·0) 0 0 0 0 = c(1) D 0 + cD(1) ; c(1) = (Rc) + D : D(1) (27) b(1) 2 3DA where (Rc) is also dened in (17). Similarly, the rst time derivative T(1) can be obtained by using the relations T = J bF−t and T = J (b 0 − pI)F−t , 0 − p(1) I)F−t T(1) = J(1) bF−t − J bl t F−t + J (b(1) (28) By introducing successively (27), (26), and (25) into (28) and taking into account the relation J(1) = J tr(l), (28) can be rewritten as 0∗ − p(1) I]F−t + J {b[tr(l)I − l t ] + [(Rc)D 0 − M : RD0 ]}F−t (29) T(1) = J [M : D(1) In this relation, all terms in the right-hand side can be calculated directly from V(1) , p(1) , V and p at time t0 . Concerning (1) , it can be considered as a prescribed function due to the given load parameter . In the case where is identied as time, (1) = 1. d , it depends on the treated boundary conditions. If t d is prescribed, For t(1) d =0 t(1) (30) d If b = b d is prescribed, the rst time derivative t(1) can be deduced from the relation t d = J b d F−t N, d = J b d (tr (l)I − l t )F−t t(1) (31) d with (30) or (31), D(1) with (25), J(1) with Finally, by substituting in (24) T(1) with (29), t(1) J tr(l), and taking into consideration J = 1, F = I for the updated Lagrangian scheme, we arrive at equations (17) after putting together all terms related to V(1) and p(1) into the left-hand side and other terms in the right-hand side. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999) 176 H.-L. CAO AND M. POTIER-FERRY ACKNOWLEDGEMENTS The authors are grateful to Professor L. Toth for his kind help to improve the text. REFERENCES 1. N. M. Wang, ‘A rigid-plastic rate-sensitive nite element method for modelling sheet metal forming processes’, in J. F. T. Pittman, O. C. Zienkiewicz, R. D. Wood and J. M. Alexander (eds.), Numerical Analysis of Forming Processes, Wiley, London, 1984, pp. 117–164. 2. S. Kobayashi, ‘Thermoviscoplastic analysis of metal forming problems by the nite element method’, in J. F. T. Pittman, O. C. Zienkiewicz, R. D. Wood and J. M. Alexander (eds.), Numerical Analysis of Forming Processes, Wiley, London, 1984, pp. 45 – 69. 3. O. C. Zienkiewicz, ‘Flow formulation for numerical solution of forming processes’, in J. F. T. Pittman, O. C. Zienkiewicz, R. D. Wood and J. M. Alexander (eds.), Numerical Analysis of Forming Processes, Wiley, London, 1984, pp. 1– 44. 4. K. Chung and R. H. Wagoner, ‘Numerical improvement of viscoplastic, non-linear, nite element analysis’, Int. J. Mech. Sci., 29, 45 –59 (1987). 5. C. Bohatier and J. L. Chenot, ‘Finite element formulation for non-steady-state viscoplastic deformation’, Int. J. Numer. Meth. Engng., 21, 1697–1708 (1985). 6. A. C. Walker, ‘A nonlinear nite element analysis of shallow circular arches’, Int. J. Solids Struct., 5, 97–107 (1969). 7. J. Connor and N. Morin, ‘Perturbation techniques in the analysis of geometrically nonlinear shells’, Proc. 1970 IUTAM Conf. High Speed Computing of Elastic Structures, Universite de Liege, 1971, pp. 683 –705. 8. L. Azrar, B. Cochelin, N. Damil and M. Potier-Ferry, ‘An asymptotic numerical method to compute the post-buckling behaviour of elastic plates and shells’, Int. J. Numer. Meth. Engng., 36, 1251–1277 (1993). 9. Y. Yokoo, T. Nakamura and K. Uetani, ‘The incremental perturbation method for large displacement analysis of elastic–plastic structures’, Int. J. Numer. Meth. Engng., 10, 503 –525 (1976). 10. N. Damil and M. Potier-Ferry, ‘A new method to compute perturbed bifurcation: Application to buckling of imperfect elastic structures’, Int. J. Engng. Sci., 28, 704 –719 (1990). 11. B. Cochelin, N. Damil and M. Potier-Ferry, ‘The asymptotic numerical method, an ecient perturbation technique for non-linear structural mechanics’, Revue Europeenne des Elements Finis, 3, 281–297 (1994). 12. M. Potier-Ferry, H. L. Cao, J. Descamps and N. Damil, ‘An asymptotic numerical method for numerical analysis of large deformation viscoplastic problems’, Comput. Struct.; submitted. 13. N. Mole, J. L. Chenot and L. Fourment, ‘A velocity based approach including acceleration to the nite element computation of viscoplastic problems’, Int. J. Numer. Meth. Engng., 39, 3439 –3451 (1996). 14. Y. J. Kim and D. Y. Yang, ‘A rigid-plastic nite element formulation considering the eect of geometric change and its application to hydrostatic bulging’, Int. J. Mech. Sci., 27, 453 – 463 (1985). 15. H. Iseki, T. Jimma and T. Murota, ‘Finite element method of analysis of the hydrostatic bulging of a sheet metal’, Bull. JSME, 20, 285 –291 (1977). 16. G. Powell and J. Simons, ‘Improved iteration strategy for non-linear structures’, Int. J. Numer. Meth. Engng., 17, 1455 –1467 (1996). 17. M. Potier-Ferry, N. Damil, B. Braikat, J. Descamps, J. M. Cadou, H. L. Cao and A. Elhage Hussein, ‘Traitement des fortes non-linearites par la methode asymptotique numerique’, Compte Rendus de l’Academie des Sciences, 324, 171–177 (1997). 18. J. Descamps, H. L. Cao and M. Potier-Ferry, ‘An asymptotic numerical method to solve large strain viscoplastic problems’. Proc. 5th int. Conf. on Computational Plasticity. Models Software and Applications, Barcelona, 1997, pp. 393 – 400. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 155–176 (1999)

1/--страниц