INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 42, 289–311 (1998) A COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS AT FINITE STRAIN WITH REFERENCE TO THIN SHELL APPLICATIONS ∗ M. ROUAINIA AND D. PERIC Department of Civil Engineering, University of Wales Swansea, Singleton Park, Swansea SA2 8PP, U.K. ABSTRACT This work extends a previously developed methodology for computational plasticity at nite strains that is based on the exponential map and logarithmic stretches to the context of isotropic elasto-viscoplastic solids. A particular form of the strain-energy function, given in terms of its principal values is employed. It is noticeable that within the proposed framework, the small strain integration algorithms, and the corresponding consistent tangent operators, automatically extend to the nite strain regime. Central to the eort of this formulation is the derivation of the closed form of a tangent modulus obtained by linearization of incremental non-linear problem. This ensures asymptotically quadratic rates of convergence of the Newton–Raphson procedure in the implicit nite element solution. To illustrate the performance of the presented formulation, several numerical examples, involving failure by strain localization and nite deformations, are given. ? 1998 John Wiley & Sons, Ltd. KEY WORDS: viscoplasticity; nite strain; localization; consistent linearization; nite elements 1. INTRODUCTION Numerical modelling of nite strain plasticity problems has emerged in the last two decades or so as a crucial part of the industrial design with a wide variety of applications involving complex constitutive relations. From a purely theoretical standpoint, numerous technical papers have been published, tackling a broad range of theoretical issues such as the notion of unstressed intermediate conguration, material frame invariance, the role of the plastic spin, decomposition of strain rates and related aspects, etc. (see References 1–4 among others). On the computational side, eort has been directed toward the formulation of algorithms for integration of constitutive equations relying mostly on the so-called operator-split methodology.5 The purpose of this paper is the formulation and numerical implementation of the constitutive model for large deformations of elasto-viscoplastic solids, within the framework of the multiplicative decomposition.1 The elasto-viscoplastic constitutive equations are derived with reference to the intermediate conguration. In addition, the present formulation is restricted to isotropic yield criterion. The logarithmic stretches are employed as strain measures in conjunction with an exponential approximation of the viscoplastic ow rule, in the form rst proposed by Weber et al.6 As ∗ Correspondence to: D. Peric. Department of Civil Engineering, University of Wales Swansea, Singleton Park, Swansea SA2 8PP, U.K. E-mail: D.Peric@Swansea.ac.uk CCC 0029–5981/98/020289–23$17.50 ? 1998 John Wiley & Sons, Ltd. Received 26 November 1996 Revised 27 August 1997 290 M. ROUAINIA AND D. PERIC an immediate consequence, the eects of nite strains appear only at kinematic level and many properties of small strain algorithms are recovered. In particular, the standard return mapping algorithm and the corresponding consistent tangent operators of the small strain situations can be directly employed. In this work, the eort is directed at the formulation of a complete computational model for nite strain elasto-viscoplastic solids and the derivation of the associated consistent tangent modulus. This is required by the Newton–Raphson solution of the associated incremental boundary value problem. Here, the explicit expression of the tangent modulus is obtained by using the concept of the consistent linearization of the weak form of the momentum balance equations. The key result is the achievement of the asymptotically quadratic rate of convergence of the equilibrium iterations. Accounting for nite strain rate-independent elasto-plastic deformations in numerical simulations often leads, as a result, to the emergence of strain localization phenomena. These phenomena, which appear as a narrow region of intense deformations, frequently form a precursor to failure mechanism. Numerical investigations of strain localization have shown that the solutions exhibit strong mesh dependency. However, it is commonly known that the use of the rate-dependent viscoplastic solids as regularization mechanism, preserves the ellipticity of the governing equations and provides a mesh-objective solution.7 – 9 The remainder of the paper is organized as follows: Section 2 dispenses some preliminary notions of kinematics, the essential assumptions and summarizes the constitutive equations based on the Perzyna formulation.10 Issues pertaining to the algorithm for numerical integration of the constitutive equations are addressed in Section 3 along with the small strain consistent tangent operators. The closed form of the consistent tangent modulus, arising from the linearization of the principle of virtual work, is derived in Section 4. An assessment of the eectiveness of the proposed formulation is provided in Section 5. Typical problems, stretching of a perforated plate, pinched hemispherical shell, Scordelis–Lo roof subjected to its own weight and plane strain stretching of a thin sheet by a square punch, form the basis of this assessment. Finally, some conclusions are drawn in Section 6. 2. FINITE STRAIN ELASTO-VISCOPLASTIC CONSTITUTIVE MODEL 2.1. Kinematical preliminaries Denote by 0 ∈ R3 a reference conguration of the continuum body B with boundary @ 0 and I ⊂ R+ the time interval of interest. At time t0 , each material point in 0 is labelled by the co-ordinates X. The motion of B is described by e : 0 × I → R3 which maps, at the time instant t, a material point onto x = e(X; t) of the current or deformed conguration. The associated deformation gradient of the motion e is dened by: F= @e(X; t) @x = @X @X (1) The kinematics of large elastoplasticity is based on the concept of multiplicative split of the deformation gradient, F, into elastic and plastic parts:1 F = FeFp ? 1998 John Wiley & Sons, Ltd. (2) Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 291 This choice introduces the concept of a local, unstressed or relaxed intermediate conguration B t . One can obtain B t from the current conguration by a purely elastic unloading of a neighbourhood of a material particle. Under the assumption of inelastic incompressibility, i.e. det[F p ] = 1, and since det[F]¿0, it follows from (2) that det[F e ]¿0. Thus, from the polar decomposition theorem, F and F e admit the unique representation F = RU; Fe = ReU e (3) where the right total and elastic stretch tensors, U and U e , respectively, are positive-denite symmetric and the rotation tensors R and R e are proper orthogonal tensors. U and U e are related to the right total and elastic Cauchy–Green deformation tensors, C and C e by U 2 = C = FT F; (U e ) 2 = C e = (F e )T F e (4) In relation to the constitutive model which is described in the next section, the Hencky logarithmic strain E e , is used as the elastic Lagrangian strain measure E e = ln [U e ] (5) Its conjugate rotated stress measure, which will be used in the description of the isotropic elastic law, is denoted T := (R e )T c R e , where c denotes the Kirchho stress tensor. A detailed discussion of stress conjugate to logarithmic strain adopted here can be found in Reference 11. The time dierentiation of expression (1) and the use of the chain rule, yield the following additive split form of the velocity gradient, L: L = ḞF−1 = L e + L p (6) where L e and L p , the elastic and plastic velocity gradients, are, respectively, given by L e = Ḟ e (F e )−1 ; L p = F e Ḟ p (F p )−1 (F e )−1 (7) p = Ḟ p (F p )−1 is used to derive the evolution of Furthermore, the ‘modied’ plastic contribution L p F and, therefore, the intermediate conguration B t . 2.2. Isotropic elastic law Following Reference 12, we assume the existence of a quadratic strain energy function W in e (i = 1; 2; 3). In fact, it is generally the form of a scalar symmetric function of its stretches (i) accepted that a suciently general constitutive model with possibly wide range of applications, including metal plasticity, may be dened by employing the strain function W . This is given by W (1e ; 2e ; 3e ) = [ln(1e ) 2 + ln(2e ) 2 + ln(3e ) 2 ] + 12 ln(J e ) 2 (8) where and are positive material constants and (J e ) = 1e 2e 3e is the Jacobian. After applying the standard procedure, the following hyperelastic constitutive equation is obtained: T= @W = He : Ee @E e (9) where H e = 2l + (I ⊗ I) denotes the fourth-order isotropic constant elastic tensor, with I and l given in the component form as Iij = ij and Iijkl = 12 (ik jl + il jk ). ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 292 M. ROUAINIA AND D. PERIC 2.3. Viscoplastic constitutive equations This section provides a brief review of the equations governing a material whose elastoviscoplastic model is based on a Von-Mises yield criterion and power law isotropic hardening.10 The model relies on the existence of the viscoplastic potential (T; T ) where (T ) : 0 × I → R is the uniaxial yield stress conjugate to the equivalent viscoplastic strain E vp ∈ R. Since the attention is restricted to the isotropic case and a single scalar hardening variable with linear hardening law, the constitutive equations can be represented by f(T; T ) = F(T) − T p F(T) = 3J2 (T) (10) (11) T (E vp ) = 0 + H E vp (12) sym[Ḟ p (F p )−1 ] = h(T; T )i@T f(T; T ) (13) E˙ vp = T˙ H −1 = h(T; T )i (14) where H is the standard isotropic linear hardening parameter and 0 is the initial uniaxial yield stress. In equation (13), ∈ (0; ∞) is the uidity parameter, h·i is the ramp function dened as hxi = (x + |x|)=2 and @T f(T; T ) = @f(T; T )=@T denotes the gradient of the yield function. Among the various choices for the viscoplastic ow potential, (T; T ), the following form for isotropic hardening is adopted: N F(T) (T; T ) = −1 (15) T where N is material parameter representing the rate sensitivity of the material. The evolution problem described by (10)–(14) and (15), has a rm experimental basis and is widely accepted as a description of rate-dependent deformations of solids. Additional formulative details of the elasto-viscoplastic model used in the present study are given in Reference 13. It is emphasized that the evolution of the plastic deformation gradient F p ( ow rule) is dened by assuming the hypothesis of normality of the viscoplastic ow expressed by equation (13). 3. NUMERICAL INTEGRATION ALGORITHM Assuming that an incremental displacement u is given, the problem consists of formulating numerical procedures for updating the known variables {Fne ; Fnp ; Tn ; E vp n } at time tn to obtain p e ; Fn+1 ; Tn+1 ; E vp } at time t . For this purpose, the widely used algorithm based on oper{Fn+1 n+1 n+1 ator split methodology5 which involves two basic steps within a load increment, namely elastic predictor and plastic corrector, is adopted. In the rst step, a purely elastic trial stress is computed; if the yield condition is violated, the visco-plastic correction step based on the Newton–Raphson procedure (in terms of (Tn+1 ; E vp n+1 )) is performed on the xed intermediate conguration, such that the nal stress is fully consistent with the constitutive integration algorithms. It should be pointed out, however, that for a low rate-sensitive material (N ¿5) which causes the viscoplastic evolution problem to be sti, these algorithms fail to converge for moderately large load steps. In this case, a special numerical treatment is used consisting in taking the values (T̃n+1 ; E˜ vp n+1 ) of ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 293 Box 1. Algorithm for integration of elasto-viscoplastic constitutive equations (i) Initialization (0) = Fn Fn+1 (ii) Update current conguration and deformation gradient Fu = 1 + Grad u; (i) Fn+1 = Fu Fn+1 (iii) Evaluate elastic trial logarithmic strain e trial = Fn+1 (Fnp )−1 Fn+1 e trial e trial T e trial Cn+1 = (Fn+1 ) Fn+1 e trial e trial T e trial 1=2 Un+1 = [(Fn+1 ) Fn+1 ] e trial e trial e trial −1 Rn+1 = Fn+1 (Un+1 ) e trial e trial En+1 = ln[Un+1 ] (iv) Perform the small strain update : GOTO Box 2 (v) Update the rst Piola–Kirchho stress and intermediate conguration. e e trial = Rn+1 Rn+1 e e Pn+1 = Rn+1 Tn+1 (Rn+1 )T (Fn+1 )−T p Fn+1 = exp[ t n+1 @T fn+1 ]Fnp the rate-independent solution as the initial conditions to evaluate the viscoplastic response.13 This procedure is additionally complemented by a line search method (see References 14 –16). The governing equations for one load step of the Newton–Raphson iterative procedure for the elasto-viscoplastic material are summarized in Boxes 1 and 2. As a crucial point in the derivation of the algorithms, an exponential approximation is used in the discretization of the ow rule in the plastic corrector step rst suggested by Weber and Anand6 (see also References 17–19, 12). p Fn+1 = exp [t n+1 @T fn+1 ] Fnp (16) e trial e trial T Furthermore, if the material is assumed to be isotropic, with the consequence that Cn+1 = (Fn+1 ) e trial Fn+1 and @T fn+1 commute, it can be easily shown that the above approximation (16) results in e e trial En+1 = En+1 − t n+1 @T fn+1 ? 1998 John Wiley & Sons, Ltd. (17) Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 294 M. ROUAINIA AND D. PERIC Box 2. Stress updating procedure-Small strain (i) Initialize k = 0; (0) E vp n+1 = 0; (0) trial Tn+1 = Tn+1 (ii) Check yield condition and evaluate residuals (k) (k) vp (k) = f(Tn+1 ; T (E n+1 )) fn+1 e −1 (k) (k) trial − Tn+1 ) (H ) : (Tn+1 @T fn+1 (k) (k) + tn+1 Rn+1 = (k) vp −1 (E vp n+1 − E n ) (k) (k) IF k Rn+1 k ¡R and fn+1 ¡f THEN Set (•)n+1 = (•)(k) n+1 and EXIT ELSE (iii) Compute algorithmic modulus (k) n+1 −1 = vp TT n+1 E Tn+1 E vp T n+1 E vp E vp n+1 −1 where (k) (k) (k) (k) e −1 2 + t(k) TT n+1 = (H ) n+1 @TT fn+1 + t@F n+1 @T fn+1 ⊗ @T fn+1 vp (k) E Tn+1 = t@Evp (k) n+1 @T fn+1 vp (k) E T n+1 = t@F (k) n+1 @T fn+1 vp vp E E n+1 = 1 − t@Evp (k) n+1 (iv) Evaluate increment of stresses and internal variable (k+1) −1 Tn+1 (k) (k) : Rn+1 vp (k+1) = − n+1 E n+1 (v) Update stresses and internal variable (k+1) (k) (k+1) Tn+1 Tn+1 Tn+1 vp (k+1) = vp (k) + vp (k+1) E n+1 E n+1 E n+1 Set ? 1998 John Wiley & Sons, Ltd. k =k + 1 and GOTO (ii) Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 295 One of the most striking features of the adoption of the logarithmic strain measure along with the exponential approximation for the viscoplastic ow rule, is the preservation of the standard return algorithms of the small strain theories. It should be noted that, as consequence of the exponential approximation, the incompressibility viscoplastic ow for the pressure-insensitive yield criterion is also preserved. 3.1. Small strain tangent operator From the computational standpoint, an important issue is the derivation of the small strain tangent operator in a manner consistent with the constitutive integration algorithm of Box 2. Simo and Taylor have shown, in their pioneering work20 concerning plasticity, that such an operator preserves the important property of asymptotic quadratic convergence in the Newton–Raphson iterative procedure. By contrast, the continuum modulus is known to destroy this important property and frequently leads to a divergent iterative process, particularly when the load step size is larger. After rather lengthy but straightforward algebraic manipulation, the following explicit expression b of the present elasto-viscoplastic model, is obtained for the small consistent tangent operator H, (see Reference 13 for details): h i−1 b = @Tn+1 = TT − TEvp (Evp Evp )−1 Evp T (18) H n+1 n+1 n+1 n+1 e @En+1 vp vp vp vp E T TE E E where TT n+1 , n+1 , n+1 and n+1 are given in Box 2. 4. LINEARIZED WEAK FORMULATION: CONSISTENT TANGENT MODULUS The full Newton–Raphson iteration scheme is characterized by a convergence rate that is asymptotically quadratic. In order to preserve this feature, it is crucial that numerical models for the description of general non-linear behaviour of solids are consistently linearized. At time tn+1 , the weak form of the momentum balance, known as the principle of virtual work, can be written, with respect to the initial conguration, as Z Z G(en+1 ; W) = (Pn+1 : ∇0 W − bn+1 ·W) d 0 − tn+1 ·W d 0 = 0; ∀W ∈ V (19) 0 @ 0 where P = c·(F)−T is the rst Piola–Kirchho stress, fn+1 and tn+1 are, the body forces and the surface tractions specied on the portion @ 0 of the boundary, respectively. The symbol ∇0 denotes the material gradient and V is the space of the admissible virtual displacements. The linearization of (19) provides the basis of the Newton-type iterative schemes for solution of the momentum balance equations. In this context, the directional derivative (see Reference 16) of the internal virtual work plays an essential role. The virtual work of the internal forces Z Pn+1 : ∇0 W d 0 (20) G int (en+1 ; W) = 0 has its derivative in the direction of an incremental displacement u given by the formula Z DG int (en+1 ; W)·u = ∇0 W : Dn+1 : ∇0 (u) d 0 (21) 0 ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 296 M. ROUAINIA AND D. PERIC where Dn+1 is the consistent tangent modulus given by Dn+1 = @Pn+1 @Fn+1 (22) Remark 4.1. It should be observed that the nite strain constitutive model described in this section is formulated in the material framework. This is in contrast to the model described by de Souza Neto and Peric21 which is phrased in spatial co-ordinates. Although both models ultimately lead to identical results, the present model is more convenient for description of large deformations of bodies which show certain preferred directions. Typical problems from this class arise in modelling of nitely deformed shells. 4.1. Exact consistent tangent modulus In order to obtain a step by step derivation of the exact expression of the consistent tangent modulus, Dn+1 , it is convenient to rewrite the rst Piola–Kirchho stress Pn+1 in the alternative form Pn+1 = c·(Fn+1 )−T = Fn+1 ·Qn+1 (23) e trial −1 e trial −T Qn+1 = (Fnp )−1 (Un+1 ) Tn+1 (Un+1 ) (Fnp )−T (24) where Qn+1 is given by The consistent tangent modulus in component form, follows straightforwardly as Dijkl = ik Qlj + (Fn+1 )im @(Qn+1 )mj @(Fn+1 )kl (25) Since Qn+1 is a function of Tn+1 , the rotated stress tensor has to be derived with respect to the total deformation gradient Fn+1 . This is accomplished by application of the chain rule as follows: e trial @(Cn+1 ) @(Tn+1 ) b =H:L: @(Fn+1 ) @(Fn+1 ) (26) b is the small strain consistent tangent operator given by equation (18) and the fourth-order where H tensor L is dened as L := e trial @(En+1 ) e trial @(Cn+1 ) (27) e trial . AddiUsing equation (5), the explicit form for L is the derivative of the tensor logarithm of Un+1 tionally, the remaining derivative term on the right-hand side of equation (26) follows easily from e trial e trial T e trial := Fn+1 ·Fn+1 = (Fnp )−T · the algorithmic denition of the elastic trial Cauchy–Green tensor Cn+1 p T ·Fn+1 ·(Fn )−1 , i.e. Fn+1 e trial @(Cn+1 ) e trial e trial = (Fnp −1 )li (Fn+1 )kj + (Fn+1 )ki (Fnp −1 )lj (28) @(Fn+1 ) ijkl To completely determine equation (25), it still remains to obtain an explicit formula for the e trial , with respect to the total deformation gradient derivative of the right trial stretch tensor, Un+1 ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 297 Fn+1 . Applying again the chain rule, yields e trial −1 e trial @(Un+1 @(Cn+1 ) ) =M: @(Fn+1 ) @(Fn+1 ) (29) in which the components of the right-hand side term are given in (28), and M is dened as M := e trial −1 @(Un+1 ) e trial @(Cn+1 ) (30) e trial e trial 1=2 By exploiting the fact that Un+1 = [Cn+1 ] , the fourth-order tensor M is the derivation of e trial the inverse of the square root of Cn+1 . For more information regarding the computation of the tensors L and M, as case of the derivative of general tensor function, the reader is referred to References 21 and 22. Finally, with all expressions involved in equation (25) at hand, the explicit formula of the Cartesian components of the consistent tangent modulus follows as b : L + K : M) : N]ijkl Dijkl = ik Qlj + [(R : H (31) Here the fourth-order tensors R, K and N are conveniently dened as e e trial e trial −1 Kijkl = [(Rn+1 )(Tn+1 )]il (Fnp −1 )jk + (Fn+1 )ik (Tn+1 )lm [(Fnp −1 )(Un+1 ) ] jm e e trial −1 Rijkl = (Rn+1 )ik [(Fnp −1 )(Un+1 ) ] jl e trial Nijkl = (Fnp −1 )li (Fn+1 )kj + (32) e trial (Fn+1 )ki (Fnp −1 )lj Examination of the closed formula (31) clearly shows that the symmetric consistent tangent modulus D, contains both geometric and material nonlinearities contributions. The material related b The remaining terms are contributions appear only through the small strain consistent tangent H. related exclusively to the kinematics of large strains. 5. NUMERICAL EXAMPLES To demonstrate the performance of the elasto-viscoplastic model in conjunction with the numerical scheme discussed earlier, numerical examples are presented in this section. All numerical simulations are performed by using an incremental solution based on the backward Euler time-stepping scheme described in Sections 3 and 4 and the Newton–Raphson procedure. It is noted that the convergence properties of the Newton–Raphson procedure are ensured by using the consistent tangent operator, given in Sections 3.1 and 4.1. The convergence of the incremental nite element solution is established on the basis of the standard Euclidean norm of the out-of-balance forces. ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 298 M. ROUAINIA AND D. PERIC Figure 1. Perforated plate subjected to stretching: geometry, nite element mesh and material properties 5.1. Stretching of a perforated plate The problem consists of a two-dimensional thin plate (plane stress)† subjected to stretching. Due to symmetry of the problem under consideration, only one quarter of the plate with appropriate boundary conditions is modeled. The geometry and material properties are shown in Figure 1. The half-length of the perforated plate is L0 = 18·0, the half-width is B0 = 10·0; the radius of the hole is R0 = 5·0 and the thickness is t0 = 1·0. The material is characterised by an isotropic hardening viscoplastic response. The following properties are assigned to the plate: elastic modulus is E = 70, Poisson’s ratio is 0·2, yield stress is 0 = 0·243 and hardening parameter H = 0·2. The nite element idealization consists of 613 three-noded constant strain membrane elements. Loading is performed by controlling the displacement U2 at the top of the plate. Figure 2 shows plots of the load versus edge displacement for dierent values of the rate sensitivity exponent N = 1; 10 and 100. The average deformation rate investigated ranges between (V2 =L0 ) = 0·555 × 10−4 (1=s) and (V2 =L0 ) = 5·555 × 10−1 (1=s) and the rate being increased by a factor of ten for every new simulation. After the peak load, the load-displacement curve drops steeply. For high rate-sensitive material, N = 1, the plots highlight the eects of the deformation rate on the increase of the maximum value of the force. For low rate-sensitive material, N = 100; and for the whole range of the deformation rates, only a small dierence in the force–displacement is observed. These plots also show that the force–displacement obtained using a very low rate of deformation (V2 =L0 ) = 0·555 × 10−4 (1=s) is nearly equal to the rate-independent solution. Figures 3 and 4 show evolution of the viscoplastic zone in the plate together with the corresponding typical deformed meshes at three dierent load levels. It can be noticed that the localization † The constitutive equations and associated numerical scheme are developed directly in the plane stress form with 33 = 32 = 31 = 0, as described in Reference 15 ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 299 Figure 2. Stretching of a perforated plate: load–displacement curves for various rates of deformation and several values of the rate sensitivity exponent ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 300 M. ROUAINIA AND D. PERIC Figure 3. Stretching of a perforated plate: evolution of viscoplastic zone, plotted over initial conguration, at various stages of loading: (a) U2 = 1·15; (b) U2 = 3·65; (c) U2 = 6·15 Figure 4. Stretching of a perforated plate: deformed nite element meshes for dierent load steps: (a) U2 = 1·15; (b) U2 = 3·65; (c) U2 = 6·15 of plastication, which takes place around the hole, results in further narrowing of the weakest section of the plate and forms a precursor to failure. This simulation is successfully completed in 48, 82 and 216 load steps for N = 1; 10 and 100, respectively. A standard check of the residual forces, based on the Euclidean norm, was carried out with a convergence tolerance of 1 × 10−6 . The number of iterations required to meet this tolerance was about 3–6 iterations for each load step. As expected, the convergence rate of the Newton–Raphson method clearly exhibits the quadratic rates of asymptotic convergence whenever the solution is within the radius of convergence. This fact is presented in tabular form for four typical load steps and the rate sensitivity of material N = 100 in Table I. For the same load steps, this fact is further elucidated graphically in Figure 5 where the residual norm kRk k of the iteration k is plotted against the residual norm kRk+1 k of the subsequent iteration k +1. Observe that the slope ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 301 COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS Table I. Perforated plate: residual norms for four typical load steps for N = 100 and average rate of deformation V2 =L0 = 0·555 × 10−2 (1=s) Iter 1 2 3 4 5 6 U2 = 0·15 U2 = 0·65 U2 = 1·15 U2 = 3·65 2·046703E−01 0·351815E−01 0·146249E−02 0·721279E−05 0·182781E−09 — 1·458925E−01 0·566755E−01 0·462543E−02 0·141538E−03 0·154090E−06 0·184024E−10 0·732901E−01 0·370454E−01 0·267615E−02 0·623028E−05 0·320327E−09 — 0·776389E−01 0·759442E−01 0·442407E−03 0·901602E−07 — — Figure 5. Stretching of a perforated plate: norms of the residual for four typical load steps corresponding to the rate sensitivity exponent N = 10 1 : 2 shown in the gure corresponds to a quadratic convergence. As mentioned earlier, these results were achieved using the closed-form consistent tangent operator given in Sections 3.1 and 4.1. 5.2. Pinched hemispherical shell A pinched hemispherical shell loaded by two concentrated forces, 90◦ apart and characterized by an isotropic hardening, is considered. Again, due to symmetry of the problem, only one quarter of the hemispherical shell is considered in the nite element simulation. The geometry, the nite element mesh, consisting of 128 triangular Morley thin shell elements, and the material properties are shown in Figure 6. The radius of the shell is R = 10 and the thickness is t0 = 0·5. The material is modelled with the following parameters : elastic modulus is E = 10, Poisson’s ratio is = 0·2, yield stress is 0 = 0·02 and hardening parameter is assumed to be equal to H = 6·0. This example has been analysed by Simo and Kennedy23 by considering the stress resultant shell theory of elasto-plasticity. Figure 7 shows the sequence of deformed conguration corresponding to dierent values of imposed concentrated force P = 0·03 × 10−01 , 0·10 × 10−01 , 0·20 × 10−01 and 0·30 × 10−01 . It is ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 302 M. ROUAINIA AND D. PERIC Figure 6. Pinched hemispherical shell: geometry, nite element model and material parameters Figure 7. Pinched hemispherical shell: deformed nite element meshes at various steps of loading: (a) P = 0·03 × 10−01 ; (b) P = 0·10 × 10−01 ; (c) P = 0·20 × 10−01 ; (d) P = 0·30 × 10−01 worth noticing that the displacement of the node 289 (inward) in the nal conguration is almost equal to the radius of the hemispherical shell. Figure 8 shows the load–displacement curves for dierent uidity parameter ranging from 0·2 × 10−6 (1=s) to 0·2 × 10−2 (1=s). The range of the rate sensitivity exponent investigated is N = 1; 10 and 100. Once again, the inuence of the material rate sensitivity is demonstrated by increase of the force for the high rate-sensitive material represented by N = 1. Only a minor difference in the overall force–displacement curves is observed for the low rate-sensitive material characterized by N = 100 over the whole range of the uidity parameter . This simulation is successfully completed in 21 load steps. The number of iterations required to meet the convergence tolerance was about 3– 4 iterations for each load step. Results of the ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 303 Figure 8. Pinched hemispherical shell: load–displacement curves for various values of the uidity parameter and several values of the rate sensitivity exponent ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 304 M. ROUAINIA AND D. PERIC Table II. Pinched hemispherical shell: residual norms for four typical load steps for N = 10 and uidity parameter = 0·2 × 10−3 (1=s) Iter P = 0·5 × 10−2 P = 0·10 × 10−1 P = 0·24 × 10−1 P = 0·30 × 10−1 1 2 3 4 0·410084E−01 0·389172E−02 0·138051E−03 0·532599E−06 0·221076E−01 0·142509E−02 0·873790E−03 0·251861E−05 0·281587E−01 0·104769E−02 0·611148E−03 0·842952E−06 0·191048E−01 0·120899E−02 0·249249E−03 0·219909E−06 Figure 9. Scordelis–Lo roof: geometry, material parameters and loading conditions Euclidean norm of the global residual are shown in Table II and suggest that, again, a quadratic rate of convergence is attained by the Newton–Raphson iterative procedure. 5.3. Elasto-viscoplastic Scordelis–Lo roof subjected to gravity loading A Scordelis–Lo roof composed of elasto-viscoplastic material and subjected to gravity loading is next considered. Due to symmetry considerations, only one-quarter of the roof, with appropriate boundary conditions, has been discretized. The geometry and the material properties for this example are shown in Figure 9. The half-length of the roof is L = 7·6 m, radius is R = 7·6 m and thickness is t0 = 0·076 m. The selected material properties of the roof are: elastic modulus is E = 2·1 × 104 N=mm 2 , Poisson’s ration is = 0·2, yield stress is 0 = 4·2 N=mm 2 and hardening parameter H = 0·0. The gravity type loading is performed with the reference value of q0 = 4 kN=m 2 . The deformations are restricted to be symmetric along the line X1 = 0 and X2 = 0. Accordingly, the boundary conditions can be expressed as U1 = 0 2 = 0 at X1 = 0 U2 = 0 1 = 0 at X2 = 0 U1 = U3 = 0 2 = 0 at X2 = L where Ui are displacements in Xi directions and j are rotations about Xj -axis. The gravitational force q=q0 applied to the roof is plotted against the vertical displacement of point A in Figure 10, obtained with a mesh of 12 × 12 triangular Morley elements. Simulations are carried out for rate sensitive material exponent N = 1 and for three dierent values of the uidity ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 305 Figure 10. Scordelis–Lo roof: Load-displacement curves for point A parameter = 2·0 × 10−2 ; 2·0 × 10−3 and 2·0 × 10−4 . For comparison purposes, the underlying rate-independent solution obtained by Peric et al.24 is also shown. From this example, the following observations, already drawn from the previous tests, may be noted again: (a) the major inuence of the rate of deformation, which is altered by varying the uidity parameter , is demonstrated by the eect on the maximum value of the gravity force; (b) the descent of the force-displacement curve observed in the rate-independent case is completely averted for the higher rate sensitive material, N = 1. For the rate-independent solution (see Figure 5.8 in Reference 24), a typical localization behaviour, as indicated by the sharp drop in the force-displacement curve, starts through the initiation of plastic hinge around the X2 axis, followed by a second hinge in the X1 direction, leading to failure. Not surprisingly, the failure mode is quite dierent in the rate-dependent case. Typical deformed meshes of the roof in Figure 11 show that the appearance of localization and associated failure is less pronounced with the increase of the rate sensitivity. This is evident from Figure 11, where the formation of the second plastic hinge in X1 direction across the mid-section of the roof is clearly delayed in comparison to the rate-independent results. 5.4. Plane strain stretching of a thin sheet by a square punch A thin sheet stretched by a rectangular punch, with approximately plane strain conditions along the centre section perpendicular to the longer side of the punch, is next considered. The analysis is performed by employing a 3-D formulation, restricting the deformations to be symmetric along the line X1 = 0 and imposing plane strain boundary conditions in the X2 direction. The geometrical characteristics and material properties are shown in Figure 12. A mesh of 240 Morley thin shell nite elements is considered for the half of the blank, and 64 triangular ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 306 M. ROUAINIA AND D. PERIC Figure 11. Scordelis–Lo roof: Deformed conguration of the shell at four dierent load steps: (a) q=q0 = 1·248; (b) q=q0 = 1·688; (c) q=q0 = 2·008; (d) q=q0 = 2·408 Figure 12. Plane strain stretching of a thin sheet by a square punch: geometry and materials characteristics facet elements are used to discretize both the surfaces of the die and the punch. The problem is investigated for a coecient of friction f = 0·25 and uidity parameter = 0·002 s−1 . Figure 13 shows plots of the punch force versus punch displacement for dierent rate sensitivity exponent N = 1; 10 and 50. The loading program is altered by varying the punch speed between Vp = 0·004 – 4 m=s, increasing the speed by a factor of ten for every new test. As expected, the major inuence of the rate of loading, demonstrated by the eect on the maximum value of the punch force, is observed for the highly rate-sensitive material represented by N = 1. No drop in the punch force is observed in the simulation for the high rate sensitive material characterised by N = 1 and ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 307 Figure 13. Plane strain stretching of a thin sheet by a square punch: punch force versus punch displacement curves for various punch speeds and several values of the rate sensitivity exponent ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 308 M. ROUAINIA AND D. PERIC Figure 14. Plane strain stretching of a thin sheet by a square punch: Deformed nite element meshes at punch displacement Dp = 4 mm, Dp = 8 mm, Dp = 12 mm, Dp = 16 mm and Dp = 20 mm Figure 15. Plane strain stretching of a thin sheet by a square punch: distribution of the true strain plotted over the initial conguration for punch speed Vp = 0·04 m=s: (a) top layer; (b) bottom layer 10, as seen in Figure 13. For low rate-sensitive material (N = 50), the force–displacement diagram displays a drop after achieving a maximum value, and approaches the rate independent solution for slow loading conditions (Vp = 0·4 × 10−2 m=s). The deformed conguration for dierent values of the punch displacement are presented in Figure 14. The distribution of the true strain in the X1 direction is shown in Figure 15 for the punch speed Vp = 0·04 m=s and the rate sensitivity N = 1 and 10 for punch displacement Dp = 12 mm. For comparison, the experimental results of Stoughton25 and the underlying rate-independent solution are also presented. For the rate-independent solution, which is evidently in agreement with the experimental results, the strain accumulates in a narrow band reaching high levels thus reecting the tendency to localisation behaviour. However, the appearance of localisation is averted in the higher rate-sensitive material characterized by N = 1. ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 309 COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS Table III. Plane strain stretching of a thin sheet by a square punch: residual norms for four typical load steps for N = 1 and uidity parameter = 0·2 × 10−2 (1=s) Iter 1 2 3 4 5 Dp = 4 mm Dp = 8 mm Dp = 12 mm Dp = 16 mm 0·711233E + 04 0·786184E + 03 0·443103E + 03 0·311620E + 01 0·457315E − 04 0·571463E + 04 0·870780E + 03 0·222233E + 03 0·700009E + 00 0·127989E − 05 0·640930E + 04 0·905270E + 03 0·136072E + 03 0·221860E + 00 0·231715E − 06 0·709472E + 04 0·619294E + 03 0·280641E + 03 0·110341E + 01 0·508533E − 05 Figure 16. Plane strain stretching of a thin sheet by a square punch: norms of the residual for ve typical load steps In our numerical model, throughout the process of deformation, the contact region between blank and punch has been restricted to a small area around the punch radius. This creates high normal and frictional forces which suppress sliding of the sheet over the punch. This trend, which is reported by Stoughton,25 is closely related to the bending eects which introduce an additional restraint and highly inuence the overall deformation pattern. It should be noted that, quadratic convergence is observed in a typical load step. As an illustration, values of the Euclidean norm of the global residual are reported in Table III and Figure 16 within several typical punch displacements. These results provide a clear indication of the quadratic rates of asymptotic convergence near the solution point. The number of iterations needed to meet the convergence tolerance of 1 × 10−3 is about 4 – 7 per increment. ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) 310 M. ROUAINIA AND D. PERIC 6. CONCLUSIONS Within the ecient hyperelastic-based multiplicative nite strain framework, a model for large elasto-viscoplastic deformations has been proposed and applied to the large-scale problems involving failure by strain localization. A particularly simple algorithm for integration of the constitutive equations, has been obtained by employing logarithmic strain as a kinematic variable and an exponential approximation of the viscoplastic ow rule. In this case the small strain format of the stress updating procedure is preserved and the eects of nite strain appear only at the kinematic level independently of the constitutive model. On the basis of the consistent linearization of the incremental boundary value problem, the explicit expression of the tangent modulus has been derived. It relies on the closed formula for the derivative of a class of isotropic tensor functions. The performance of the described formulation and the accuracy of the numerical results have been demonstrated for four representative problems. It has been found that both, the material rate sensitivity and the strain rate have a pronounced eect on the appearance of strain localisation and on the description of the failure mechanism. The use of the consistent tangent operator given in Sections 3.1 and 4.1, has enabled quadratic rates of convergence to be obtained in the equilibrium iterations. This increases the speed of the computations by requiring only an average of three to six iterations for each load step in all the numerical examples. REFERENCES 1. E. H. Lee, ‘Elastic-plastic deformation at nite strains’, J. Appl. Mech., 36, 1– 6 (1969). 2. J. Mandel, ‘Equations constitutives et directeurs dans les milieux plastiques et viscoplastiques’, Int. J. Solids Struct., 9, 725 – 740 (1973). 3. S. Nemat-Nasser, ‘On nite elasto-plasticity’, Int. J. Solids Struct., 18(10), 857 – 872 (1982). 4. J. C. Simo, ‘A framework for nite elastoplasticity based on maximum plastic dissipation and multiplicative decomposition: Part I. Continuum formulation’, Comput. Meth. Appl. Mech. Engng., 66, 199 – 219 (1988). 5. J. C. Simo and T. J. R. Hughes, ‘General return mapping algorithms for rate independent plasticity’, in C. S. Desai et al. (eds.), Constitutive Laws for Engineering Materials: Theory and Applications, Elsevier, Amsterdam, 1987, pp. 221– 231. 6. G. Weber and L. Anand, ‘Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solid’, Comput. Meth. Appl. Mech. Engng., 79, 173 – 202 (1990). 7. B. Loret and J. H. Prevost, ‘Dynamic strain localisation in elasto(visco)-plastic solids, Part I: General formulation and one dimensional examples’, Comput. Meth. Appl. Mech. Engng., 83, 247 – 273 (1990). 8. A. Needleman, ‘Material rate-dependence and mesh sensitivity in localisation problems’, Comput. Meth. Appl. Mech. Engng., 67, 69 – 86 (1988). 9. L. J. Sluys, J. Block and R. de Borst, ‘Wave propagation and localisation in viscoplastic media’, in D. R. J. Owen et al. (eds.), Proc. 3rd Int. Conf. on Computational Plasticity, Pineridge Press, Swansea, 1992, pp. 299 – 312. 10. P. Perzyna, ‘Thermodynamic theory of viscoplasticity’, in Advances in Applied Mechanics, Vol. 11, Academic Press, New York, 1971, pp. 313 – 354. 11. A. Hoger, ‘The stress conjugate to logarithmic strain’, Int. J. Solids Struct., 23, 1645 –1656 (1987). 12. D. Peric, D. R. J. Owen and M. E. Honnor, ‘A model for nite strain elasto-plasticity based on logarithmic strains’, Comput. Meth. Appl. Mech. Engng., 94, 35 – 61 (1992). 13. D. Peric, ‘On a class of constitutive equations in viscoplasticity: formulation and computational issues’, Int. J. Numer. Meth. Engng., 36, 1365 –1393 (1993). 14. M. A. Criseld, Non-linear Finite Element Analysis of Solids and Structures. Vol. 1: Essentials, Wiley, Chichester, 1991. 15. M. Dutko, D. Peric and D. R. J. Owen, ‘Universal anisotropic yield criterion based on superquadratic functional representation: Part I: algorithmic issues and accuracy analysis’, Comput. Meth. Appl. Mech. Engng., 109, 73 – 93 (1993). 16. J. C. Marsden and T. J. R. Hughes, Mathematical Foundations of Elasticity, Prentice-Hall, Englewood Clis, N.J., 1983. 17. A. Cuitiño and M. Ortiz, ‘A material independent method for extending stress update algorithms from small strain plasticity to nite plasticity with multiplicative kinematics’, Engng. Comput., 9, 437 – 451 (1992). ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998) COMPUTATIONAL MODEL FOR ELASTO-VISCOPLASTIC SOLIDS 311 18. A. L. Eterovic and K. J. Bathe, ‘A hyperelastic based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures’, Int. J. Numer. Meth. Engng., 30, 1099 –1114 (1990). 19. J. C. Simo, ‘Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the innitesimal theory’, Comput. Meth. Appl. Mech. Engng., 99, 61–112 (1992). 20. J. C. Simo and R. L. Taylor, ‘Consistent tangent operators for rate-independent elastoplasticity’, Comput. Meth. Appl. Mech. Engng., 48, 101–118 (1985). 21. E. A. de Souza Neto and D. Peric, ‘A computational framework for a class of fully coupled models for elastoplastic damage at nite strains with reference to the linearization aspects’, Comput. Meth. Appl. Mech. Engng., 130, 179 –193 (1996). 22. D. E. Carlson and A. Hoger, ‘The derivative of a tensor-valued function of a tensor’, Quart. Appl. Math., 44(3), 409 – 423 (1986). 23. J. C. Simo and J. G. Kennedy, ‘On the stress resultant geometrically exact shell model. Part V. Nonlinear plasticity: formulation and integration algorithm’, Comput. Meth. Appl. Mech. Engng., 99, 133 –177 (1992). 24. D. Peric and D. R. J. Owen, ‘The Morley thin shell nite element for large deformations problems: Simplicity versus sophistication’, in N. Bicanic (ed.), Proc. Int. Conf. on Nonlinear Engineering Computations, Pineridge Press, Swansea, 1991, pp. 121–142. 25. T. B. Stoughton, ‘Finite element modeling of 1008 AK sheet steel stretched over a rectangular punch with bending eects’, in N. M. Wang et al. (eds.), Computer Modeling of Sheet Metal Forming Process, AIME, New York, 1985, pp. 143 –159. ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 289–311 (1998)

1/--страниц