Computational Damage Mechanics: Application to Metal Forming Simulation

3.06 Computational Damage Mechanics: Application to Metal Forming Simulation K. SAANOUNI UniversiteŁ deTechnologie deTroyes, GSM/LASMIS, France and J...

Download PDF

3MB Sizes 0 Downloads 56 Views

Report

Full Text

3.06 Computational Damage Mechanics: Application to Metal Forming Simulation K. SAANOUNI UniversiteŁ deTechnologie deTroyes, GSM/LASMIS, France and J. L. CHABOCHE UniversiteŁ deTechnologie deTroyes, GSM/LASMIS, France and ONERA, OMSE/LCME, Chatillon Cedex, France 3.06.1 INTRODUCTION

322

3.06.2 THEORETICAL BACKGROUND

323 323 324 325

3.06.2.1 Kinematical Framework 3.06.2.2 Thermodynamical Framework 3.06.2.3 On the Ductile Damage Modeling 3.06.3 CONSTITUTIVE EQUATIONS WITH DAMAGE 3.06.3.1 Main Assumptions and State Variables 3.06.3.2 State Potential and State Relations 3.06.3.3 Dissipation Potentials and Complementary Relations 3.06.3.3.1 Thermal dissipation analysis: fully coupled heat ﬂow equation 3.06.3.3.2 Mechanical dissipation: two yield surface formulation 3.06.3.3.3 Mechanical dissipation: one yield surface formulation 3.06.3.4 The Fully Isotropic Case 3.06.3.4.1 Two-surface formulation 3.06.3.4.2 One-surface formulation 3.06.3.5 On the Frictional Constitutive Equations 3.06.3.5.1 Time-independent formulation 3.06.3.5.2 Time-dependent formulation 3.06.4 COMPUTATIONAL ASPECTS

326 326 328 329 329 329 333 334 334 337 337 338 338 339

3.06.4.1 Formulation of the IBVP 3.06.4.2 Spatial Discretization of the IBVP 3.06.4.3 Incremental Solution Procedures 3.06.4.3.1 SI resolution strategy 3.06.4.3.2 Dynamic explicit resolution strategy 3.06.4.4 Local Integration of the Fully Coupled Constitutive Equations 3.06.4.4.1 Application to the isotropic single yield surface model 3.06.4.4.2 Consistent elastoplastic tangent operator 3.06.4.4.3 On the incremental objectivity 3.06.4.5 On the Numerical Aspects of Contact-friction Problems 3.06.4.6 On Some Aspects of Remeshing and Adaptive Meshing 3.06.4.6.1 Remeshing decision 3.06.4.6.2 New mesh generation

321

339 340 342 342 345 345 346 348 349 349 350 351 351

322

Computational Damage Mechanics: Application to Metal Forming Simulation 351 352

3.06.4.6.3 Field transfer and FEA restart 3.06.4.7 Special Treatment of Fully Damaged Elements 3.06.5 NUMERICAL PREDICTION OF THE DAMAGE: SOME APPLICATIONS 3.06.5.1 On the Identiﬁcation of the Material Parameters 3.06.5.2 Some Simple Numerical Results 3.06.5.3 Damage Prediction in Sheet Metal Forming 3.06.5.3.1 Sheet metal cutting 3.06.5.3.2 Hydroforming of thin tubes 3.06.5.4 Damage Prediction in Bulk Metal Forming 3.06.5.4.1 Cold extrusion of a 3D spider 3.06.5.4.2 Hot forging of 3D billet 3.06.6 CONCLUDING REMARKS AND PERSPECTIVES 3.06.6.1 Concluding Remarks 3.06.6.2 Perspectives 3.06.7 REFERENCES

3.06.1

INTRODUCTION

Continuum damage mechanics (CDM) is widely used to model the progressive material degradation which occurs during the deformation processes, prior to macroscopic crack initiation. It has already attained a high degree of maturity and is being used for many different applications in connection with numerical simulation methods. Many attempts have been made to propose different kinds of CDM models describing different kinds of damage mechanisms for various materials under various applied loading paths (see the books by Lemaitre and Chaboche (1990), Lemaitre (1992), Krajcinovic (1984), Voyiadjis and Kattan (1999), Besson et al. (2001), and Kattan and Voyiadjis (2002), see also Chapter 2.04). Particularly, in metal-forming processes the material is subjected to large inelastic deformations, leading to the formation of high-strainlocalization zones and, consequently, to the onset of internal (voids) or surface (ﬂat cracks) microdefects. Based on the ﬁnite element analysis (FEA), many theoretical and numerical models have been proposed to simulate, as accurately as possible, the material ﬂow (strain and stress distribution), the evolving contact with friction, and the large change in the part shape during the process until its end (see the proceedings of different editions of the specialized International Conferences NUMIFORM, NUMISHEET, ESAFORM). However, much less work has been dedicated to the prediction of the ductile damage development during sheet metal or bulk metal forming. In fact, whatever the damage (microdefect) theory, there are two principal methods to numerically simulate the damage occurrence (here we exclude models using damage theories to predict the forming limit curves). The ﬁrst one, called the uncoupled approach, consists in

352 353 355 355 356 357 361 361 363 365 365 369 372

carrying out a conventional FEA simulation to predict the evolution of the part shape. By postprocessing the FEA for a given deformed conﬁguration, the obtained thermomechanical ﬁelds are used together with a ductile failure criterion, to estimate the damage distribution in the part for the given conﬁguration (Gelin et al., 1985; Hartley et al., 1989; Liu and Chung, 1990; Clift, 1992). In this ‘‘uncoupled’’ approach, the effect of damage on the elastic and the inelastic behavior is neglected. Consequently, the redistribution of the mechanical ﬁelds (stress, strain, damage, etc.) generated by the damage (voids, microcracks) initiation and growth are not correctly described. The second approach is the ‘‘coupled’’ one in the sense that, at each time step the damage increment is directly incorporated (fully coupled) in the elasto-inelastic constitutive equations. This kind of approach has been employed by many authors using damage models based either on Gurson’s type theory (Gelin and Predeleanu, 1989; Onate et al., 1988; Gelin, 1990; Aravas, 1986; Bontcheva and Iankov, 1991; Boudeau and Gelin, 1994; Bennani and Oudin, 1995; Brunet et al., 1996, 2001; Picart et al., 1998), or on CDM in the Chaboche and Lemaitre sense (Lee et al., 1985; Mathur and Dawson, 1987; Zhu et al., 1992; Zhu and Cescotto, 1995; Saanouni and Franqueville, 1999; Saanouni et al., 2000, 2001a, 2001b, 2001c). These fully coupled approaches allow one to predict not only the large transformation of the processed part as large deformations, rotations, and evolving boundary conditions, but indicate where and when the damaged zones can appear inside the formed part also. This gives a powerful tool for the industrial metal forming optimization by choosing the process technological parameters in such a manner that the damage is avoided. In some other processes—such as shaving, blanking, or orthogonal cutting—ductile

Theoretical Background damage is designed into an operation by generating controlled damaged zones (or macroscopic cracks) through the intensively deformed zone. Again, the fully coupled approaches have shown their ability to numerically simulate accurately these sheet metal operations (Abdali et al., 1995; Homsi et al., 1996; Bezzina and Saanouni, 1996; Brokken et al., 1998; Faura et al., 1998; Samuel, 1998; Saanouni and Franqueville, 1999; Saanouni et al., 2000; Hambli, 2001). This chapter is devoted to the computational aspects of CDM with application to metal forming problems. The theoretical foundations of the ﬁnite elastoplastic constitutive equations coupled to the continuum damage are reviewed and effects of damage on both elastic and inelastic behaviour are discussed (Chaboche, 1988; Saanouni et al., 1994; Voyiadjis and Kattan, 1999 and references therein). Numerical formulation of the fully coupled problem is discussed in detail and its implementation into a general-purpose ﬁnite element (FE) code is presented. It is felt that the subject matter described in this chapter constitutes a representative sample of state-of-the-art methodology used in computational CDM. Various selected examples dealing with the prediction of the ductile damage in metal forming processes are presented. It is worth noting that no attempt has been made to present a careful account of the historical developments of the subject. Likewise, the list of references should, by no means, be regarded as a complete literature survey of the computational damage mechanics. Indeed, we limit our citations to the major works on CDM applied to metal forming processes. The chapter is divided into four main sections. Section 3.06.2 is devoted to the review of basic notions related to the kinematics of ﬁnite deformations and the thermodynamics of irreversible processes with state variables. It closes by giving the main approaches used in the literature to model the ductile damage and its effect on the elastic and inelastic behavior (strong coupling). Section 3.06.3 deals with the formulation of the constitutive equations with damage effect in the framework of the thermodynamics of irreversible processes with state variables. A general initial and induced anisotropy modeling approach is presented using both singlesurface and two-surface formulations. The classical frictional constitutive equations are brieﬂy reviewed. Section 3.06.4 is devoted to the numerical and computational aspects of initial and boundary value problem (IBVP). First, the

323

standard displacement based weak variational form is introduced in the framework of nonlinear solid mechanics, including material and geometrical nonlinearities. Both explicit and implicit resolution strategies are presented and discussed in order to solve the associated partial differential equations (PDEs) in the light of their application to the highly nonlinear IBVP problems. The fully implicit local numerical integration of the coupled constitutive equations (ordinary differential equations (ODEs)) is presented and discussed in detail. Special attention is given to the calculation of the material Jacobian matrix needed to construct the consistent tangent stiffness matrix in the case of iterative implicit global resolution strategy. Finally, some numerical aspects related to the contact between solids and the adaptive remeshing problems are brieﬂy discussed and the effect of the ductile damage is highlighted. Finally, Section 3.06.5 is concerned with some applications to the prediction of the ductile damage expected to occur during metal forming processes. In fact, a complete methodology is proposed to realize ‘‘virtual’’ metal forming processes in order to predict when and where ductile damage may take place inside the working part. This methodology is applied to some sheet or bulk metal forming processes with the double objectives, namely, to avoid or to enhance the ductile damage initiation and growth.

3.06.2

THEORETICAL BACKGROUND

The purpose of this section is to present a short overview of the kinematics, thermodynamics, and the most commonly used methods in ductile DM modeling under ﬁnite deformation hypothesis.

3.06.2.1

Kinematical Framework

Consider a three-dimensional (3D) deforming body with material points, each identiﬁed with its spatial position x in the Euclidean point space at time t corresponding to the position X at an initial time t ¼ 0. The wellknown notion of an intermediate local conﬁguration, relative to which the elastic response of the body should be characterized, is used (see Lee, 1971; Mandel, 1973; Sidoroff, 1973, among many others). This leads to a local multiplicative decomposition of the transformation gradient F ¼ grad(x(X, t)) into elastic (F e) and plastic (F p) parts. Accordingly, the following classical deﬁnitions on the deformed

324

Computational Damage Mechanics: Application to Metal Forming Simulation

conﬁguration are used: F ¼ F e F p and B ¼ F F T

ð1Þ

L ¼ F’ F 1 ¼ D þ W

ð2Þ

D ¼ 12ðL þ LT Þ

ð3aÞ

W ¼ 12ðL LT Þ

ð3bÞ

and

where B is the total Eulerian left Cauchy– Green deformation tensor associated with the Cauchy stress tensor r; L is the spatial velocity gradient in the current conﬁguration, D and W are, respectively, the pure strain rate and the material spin tensors. The overdot denotes the usual material time derivative and the superscript ‘‘T’’ stands for the transpose operation. To satisfy the objectivity requirement, the so-called rotated frame formulation (RFF) or rotation neutralized description is used. This allows us to express the constitutive equations in a rotated conﬁguration from the current one by a proper orthogonal rotation tensor Q deﬁned by Dienes (1979), Johnson and Bamman (1984), and Simo and Hughes (1998), among others): ’ ¼ WQ QT Q

with Qðt ¼ 0Þ ¼ I

ð4Þ

Accordingly, for any symmetric second-order tensor T; the objective rotational derivative with respect to the rotating frame is given by DQ T ¼ T’ þ T W Q W Q T DQ t

ð5Þ

from which the classical Jaumann and GreenNaghdi rotational derivatives can be easily obtained. Alternatively, the objective rotated (by the rotation Q) tensor TQ is given by T Q ¼ QT T Q

ð6Þ

Its time and rotational derivatives are related by DQ T Q T’ Q ¼ QT DQ t

ð7Þ

Consequently, the constitutive equations are formulated in the same way as under smallstrain hypothesis and their generalization to the large-strain case is simply achieved by replacing all the tensorial variables by their rotated corresponding quantities deﬁned by Equation (6). Note that this very helpful ad hoc generalization of the inﬁnitesimal rate-independent plasticity supposes hypoelastic nature

of the elastic response (see later). Despite the serious objections from a physical and fundamental perspective, this widely used hypothesis will be assumed here, according to the conceptual theoretical and computational simplicity of the obtained constitutive equations. The second basic question posed by the ﬁnite deformation hypothesis is how the Eulerian total strain rate (see Equation (3a)) can be decomposed into elastic (reversible) and plastic (irreversible) parts. For the forming of metallic materials, it is usually supposed that elastic strains are still very small compared to the plastic ones. Consequently, and as a ﬁrst approximation, the total Eulerian strain rate tensor can be additively decomposed using DE’eJe þ Dp

ð8Þ

where e’ Je is the Jaumann derivative (rotational objective derivative) of the elastic strain tensor, and Dp is the plastic strain rate tensor deﬁned by the constitutive equations as will be shown in the next section. It is worth noting that, by considering the small elastic strain, the hyperelasticity of the elastic response is more sound. Throughout this chapter, and for the sake of notational simplicity, the subscript Q associated to any rotated tensorial variable (see Equation (6)) will be omitted. Similarly, the superscript J for the Jaumann derivative of the small elastic strain (see Equation (8)) will be omitted. 3.06.2.2

Thermodynamical Framework

The basic framework of the irreversible thermodynamics (Truesdell and Noll, 1965/ 1992; Germain, 1973; Germain et al., 1983; Lemaitre and Chaboche, 1990) with state variables will be presented using the Eulerian description on the deformed conﬁguration. It is perfectly consistent with the presentation made in Chapter 2.04 in the present series. In the thermodynamics of irreversible processes, the local state method (Germain, 1973) attributes one couple of state variables to each thermomechanical phenomenon to be analyzed namely, (ai, Ai) for ‘‘i ’’ set of different phenomena. The force-like state variable Ai is associated to the strain-like variable ai. These are usually taken as tensorial variables of order zero (scalars), one (vectors), two, or four according to the physical nature of each considered phenomenon. In this macroscopic and phenomenological approach, the state variables are classiﬁed into two families, namely, the observable and the internal variables, as will be deﬁned later for the particular case of this study.

Theoretical Background A state potential is then chosen as a continuous scalar function of the state variables ai in the strain space, or Ai in the stress space. It should be a convex and/or concave function of its argument and should contain the origin. This potential, usually taken as the Helmholtz free energy in strain space or its dual by Legendre–Fenchel transformation (free enthalpy) in the stress space, allows the deﬁnition of the so-called state relations. For example, if we start from the free energy c(ai), the derivation of these state relations (i.e., the deﬁnition of the stress-like variables Ai from c(ai)) is based on the Clausius–Duhem inequality obtained by the combination of the ﬁrst and the second thermodynamic laws: ’ i Þ þ sT’ q gZ0 r : D r cða T

ð9Þ

where r is the stress tensor (conjugated to the total strain rate tensor D), r is the density of the body, s is the speciﬁc entropy, T is the absolute temperature, g ¼ grad(T) and q is the heat ﬂux vector. All these quantities are deﬁned with respect to the Eulerian deformed conﬁguration. From the fundamental inequality (Equation (9)), both the state relations Ai ¼ r

@cðai Þ @ai

ð10Þ

and the residual or dissipation Y inequality Y ¼ r : Dp

X

Ai a’ i

i

q gZ0 T

ð11Þ

are obtained after some simple algebraic manipulations. It is worth noting that the dissipation inequality (Equation (11)) is usually decomposed into two different terms, which are separately positive or null if thermomechanical uncoupling is assumed. The ﬁrst term deﬁnes the intrinsic dissipation Yin ¼ r : Dp

X

Ai a’ i Z0

ð12Þ

i

and the second the thermal dissipation q Yth ¼ gZ0 T

ð13Þ

In order to ensure the satisfaction of the Equations (12) and (13), the existence of dissipation potentials is assumed. In the stress space, e.g., these are chosen as scalar continuous and convex functions of the stress-like variables with the strain-like variable acting as simple parameters, namely, F(Ai; ai). The kinetic constitutive equations describing the evolution of the studied phenomena are

325

derived from these potentials according to the normality rule: @F ðAi Þ a’ ¼ 7l’ @Ai

ð14Þ

where l’ is the Lagrange multiplier to be derived from the consistency condition applied to a speciﬁc yield function in the case of nonassociative plasticity formulation as will be shown later. This framework will be used in Section 3.06.3 to derive the small thermoelastic and large plastic strain model coupled to ductile damage for metal forming. 3.06.2.3

On the Ductile Damage Modeling

Ductile damage grows continuously as plastic strain increases. Many investigations have shown that ductile fracture involves, schematically, three successive damage stages (see Franc-ois et al., 1993 and references therein) which are given in the following. Void nucleation. After a certain amount of plastic deformation, microvoids are nucleated mainly at the boundaries of some favorably oriented inclusions. At this stage, the effect of these nucleated microvoids on the physical properties of the material can be neglected. Void growth. The microvoids grow through the plastically deformed matrix, under the effect of the triaxial stress state generated by the applied external load. This growth induces, generally, some change in the void shape (anisotropy) as well as a certain volume change as the distance between voids becomes smaller. At this stage the material properties (thermal, elastic, and plastic) become strongly affected by the void growth (i.e., damage). Accordingly, the coupling between the thermomechanical behavior and the voids cannot be neglected to take into account the effect of microvoids growth on the material behavior. Void coalescence. When the distance between growing voids becomes sufﬁciently small, very large plastic deformations are concentrated into small ligaments. This induces microshear bands connecting some voids, which deﬁne the onset of macroscopic crack initiation (i.e., the ﬁnal fracture of the representative volume element (RVE) of the material). At this stage, the stress ﬁelds (Cauchy stress and internal stresses) inside the macroscopic hole (i.e., the completely damaged RVE) vanish. The successive fracture of many RVEs deﬁnes the propagation stage of a macroscopic ductile crack inside the body or the structure. Phenomenologically, these ductile-damage stages, deﬁning the macroscopic crack

326

Computational Damage Mechanics: Application to Metal Forming Simulation

initiation, should take into account the following phenomena: (i) continuous degradation of physical properties (elastic, plastic, thermal, etc.); (ii) strain softening effects on yielding of the RVE inducing zero stresses at the ﬁnal fracture, i.e., when all the damage tensor components reach their critical values; (iii) decrease in the thermal conductivity inside the fully damaged RVE; and (iv) induced volume variation as the voids grow. In the wide literature concerned with the subject, two classes of coupled models are often used to describe these mechanisms of ductile damage. Gurson’s models. In this approach, the distribution of the different kinds of existing voids is represented by a single scalar variable, namely, the volume fraction of voids f deﬁned by the ratio between the accumulated volume of individual voids Vv and the total volume of the RVE at the time t, Vt: f ¼

Vv Vv r Vv ¼ ¼ t ; Vt JV0 r0 V0

0:0rf r1:0

ð15Þ

with J ¼ detðFÞ ¼

r0 Vt ¼ rt V0

ð16Þ

where r0, V0 and rt, Vt are the material density and the volume of the RVE in the initial and deformed conﬁgurations, respectively. To realize the coupling between the plastic behavior and the ductile damage of the RVE, the basic idea, pioneered by Gurson (1977), of this approach consists in introducing the variable f into the plastic potential from which the plastic strain rate is derived by Equation (14). However, the evolution equation of the variable f is not obtained by the same normality rule but postulated as the sum of three rates representing the three mechanisms discussed above: f’ ¼ f’nucleation þ f’growth þ f’coalescence

ð17Þ

The term concerning the growth stage is derived from the principle of mass conservation to get (neglecting the volume change in the virgin material) f’growth ¼ ð1 f Þ trðDÞ

ð18Þ

The two other constitutive equations concerning the nucleation and the coalescence stages are postulated, mainly based on the experimental observation (see, e.g., Gurson, 1977; Tvergaard and Needleman, 1984; Beremin, 1981; Tvergaard, 1990; Benzerga et al., 1999) or based

on some analytical or numerical arguments (Rousselier, 1987; Gologanu et al., 1993, 1994; Suquet, 1993; Ponte Castaneda and Zaidman, 1994; Leblond et al., 1995; Garajeu and Suquet, 1997; Garajeu et al., 2000, among many others). This ‘‘physically’’ based approach to model the ductile damage is often termed in the literature as a micromechanical approach. Note that, in this approach, the ductile-damage-induced anisotropy (change in the void shape) is taken into account by introducing an appropriate void shape factor in the plastic potential as proposed in Benzerga et al. (1999). Alternatively, the elastic behavior of the material is still unaffected by the void growth. CDM approach. Starting from classical CDM theory (Lemaitre and Chaboche, 1990; Lemaitre, 1992; Krajinovic, 1984; Voyiadjis and Kattan, 1999; Besson et al., 2001; see also Chapter 2.04 by Chaboche), this approach describes, phenomenologically, the above discussed three stages of ductile damage by one (or more) scalar (isotropic) or tensorial (anisotropic) internal variable. The constitutive (or evolution) equations of damage variables are derived from speciﬁc damage potentials (or yield functions) by using the effective state variables. These are deﬁned from the classical state variables using one of the three following hypotheses: strain equivalence, stress equivalence, or energy equivalence. The coupled constitutive equations of the damaged RVE are generally deduced from the same state and dissipation potentials in which the state variables are replaced by the effective state variables (Lemaitre and Chaboche, 1990; Lemaitre, 1992). In this work, and to be short, we limit the present formulation to the CDM approach, which will be applied in the next section to develop the coupled elastoplastic-damage model for metal forming. 3.06.3

CONSTITUTIVE EQUATIONS WITH DAMAGE

This section is devoted to the presentation of the ﬁnite elastoplastic-damage coupled constitutive equations for metal forming. As discussed above, the framework of thermodynamics of irreversible processes with state variables will be used within the Eulerian (spatial) description. 3.06.3.1

Main Assumptions and State Variables

The framework of the fully local formulation is used. This means that the ﬁrst displacement gradient is supposed to be sufﬁcient to describe

Constitutive Equations with Damage the kinematics of the body deformation (see Section 3.06.6.2 for the nonlocal or higher gradient formulations). The elastic components of strain rate are assumed to be inﬁnitesimal (see Equation (8)), and the formulation is made in the rotated reference frame (i.e., all tensorial variables are rotated by Equation (6) without using the subscript Q for simplicity). In the present work two different formulations are used. The ﬁrst one uses a single (or uniﬁed) yield function and dissipation potential to derive all the dissipative phenomena. The second one considers two different yield functions and dissipation potentials (Voyiadjis and Deliktas, 1999; Voyiadjis and Kattan, 1999; Hammi, 2000). The ﬁrst yield function, expressed in the Cauchy stress space, is related to the plasticity phenomenon with damage effect. The second, expressed in the damage conjugate force (Y) space, is related to the continuous damage effect. Dealing with the simple ﬁrst displacement gradient theory, two pairs of observable (or external) state variables are used and are given in the following. (i) Total strain tensor associated with the Cauchy stress tensor (e, r). (ii) Straightforwardly speaking, in Eulerian description of the ﬁnite deformation the variables (B, r) should be used if the elastic deformation gradient F e is not supposed to be inﬁnitesimal. (iii) Absolute temperature associated with the speciﬁc entropy (T, s). To describe the plastic behavior with damage in the framework of two yield surfaces, seven pairs of internal variables are taken into account:

327

(i) small elastic strain representing the elastic ﬂow associated with the Cauchy stress tensor (ee, r); (ii) normalized heat ﬂux vector associated with the gradient of the absolute temperature (q/T, g); (iii) isotropic hardening variables (r, R) representing the size of the yield surface in strain space (r) and stress space (R); (iv) tensorial kinematic hardening variables (a, X) representing the displacement of the center of the yield surface in strain space (a) and stress space (X); (v) second-order tensorial damage variables (d, Y); (vi) scalar variables (b, B) representing the isotropic increase of the damage yield surface radius in the d or Y spaces, respectively (equivalent to the isotropic hardening for plasticity); and (vii) second-order tensorial variables (c, C) describing the translation of the damage yield surface (equivalent of the kinematic hardening). It is worth noting that, when using a unique yield surface formulation, the variables (c, C) are not necessary. Suppose that the current conﬁguration at time t contains some ductile-damage distribution, i.e., a given distribution of microdefects as voids and/or microcracks; the concept of the effective stress (Chaboche, 1988; Lemaitre and Chaboche, 1990) together with the hypothesis of total energy equivalence (Saanouni et al., 1994) are used to deﬁne the effective state variables by (Figure 1) r* ¼ M1 : r and e*e ¼ MT : ee

Wt Real Damaged Configuration

(e, ), (r, R) (, X), (D, Y), (, ), (, B)

Fictitious undamaged Configuration

~ ~~ R) (~ e, ), (r,

~ ~ X) (,

Figure 1 Scheme of the total energy equivalence (Wt) hypothesis.

ð19Þ

328

Computational Damage Mechanics: Application to Metal Forming Simulation X˜ ¼ N1 : X and a* ¼ NT : a R R˜ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ and 1 jjdjj

r˜ ¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 jjdjjr

ð20Þ ð21Þ

where Mijkl(d) and Nijkl(d) are fourth-order symmetric and positive deﬁnite operators describing the damage effect on Cauchy stress and kinematic stress, respectively (Lemaitre and Chaboche, 1990; Voyiadjis and Kattan, 1999; Hammi, 2000). They reduce to fourthorder unit tensor if the damage vanishes (i.e., d ¼ 0) and they decrease as the damage tends to its critical (maximum) value at the ﬁnal fracture of the RVE. The notation ||d|| in Equation (21) deﬁnes an appropriate Euclidean norm of the second-order damage tensor. In the following, these effective state variables will be used in the state and dissipation potentials (Lemaitre and Chaboche, 1990; Saanouni et al., 1994) to derive the complete set of fully coupled constitutive equations for metal forming processes. 3.06.3.2

State Potential and State Relations

The Helmholtz free energy c(ee, a, r, d, b, c, T) is taken as a state potential. It is assumed to be a convex function of all the strain-like state variables deﬁned above. It will be additively decomposed into three contributions, namely, thermoelastic/damage, plastic/damage, and purely damage terms: rc ¼ rcted ðee ; d; TÞ þ rcpd ða; r; d; TÞ þ rcd ðb; c; d; TÞ

where r ¼ rt is the material density in the current deformed and damaged conﬁguration and the variables located after the (;) act as simple parameters. In this work, these three terms of the state potential are chosen as follows:

rcpd ða; r; d; TÞ ¼

1 * 2a

rcd ðc; b; T; dÞ ¼

1 2c

: C : a* þ :L:cþ

1 r2 2Q˜

2 1 2Gb

@c * : ee ðT T0 Þk˜ r¼r e¼K @e s¼

@c 1 ˜ e Cv ¼ k : e þ ðT T0 Þ @T r T0

ð23Þ

ð27Þ

ð28Þ

@c ˜ ¼ Qr R¼r @r

ð29Þ

Y ¼ r

@c ¼ ðY e þ Y k þ Y i Þ @d

ð30aÞ

* 1 @K @ k˜ Y e ¼ ee : : ee ðT T0 Þ : ee 2 @d @d

ð30bÞ

1 @ C˜ :a Yk ¼ a : 2 @d

ð30cÞ

1 @ Q˜ 2 r 2 @d

ð30dÞ

@c ¼L:c @c

ð31Þ

@c ¼ Gb @b

ð32Þ

Yi ¼

C¼r

It is remarkable to note that in the state relations above the main physical properties of the material are affected by the damage according to: damaged elasticity moduli operator: * ¼ M : K : MT K damaged kinematic hardening moduli operator: C˜ ¼ N : C : NT

ð24Þ

damaged thermal properties: k˜ ¼ MT : k

ð25Þ

damaged isotropic hardening modulus: Q˜ ¼ ð1 jjdjjÞQ

where Kijkl is the symmetric fourth-order elasticity operator, kij is the symmetric secondorder thermal operator, Cijkl is the symmetric fourth-order kinematic hardening operator, Q is the scalar isotropic hardening modulus, Lijkl is the symmetric fourth-order ‘‘kinematic’’ damage operator, and G is the ‘‘isotropic’’ damage modulus. They are all assumed to be

ð26Þ

@c ¼ C˜ : a X¼r @a

B¼r ð22Þ

rcted ðee ; d; TÞ ¼ 12 e*e : K : e*e ðT T0 Þk : e*e Cv 12 r ðT T0 Þ2 T0

symmetric, positive deﬁnite, and temperaturedependent parameters. Note that T0 is the reference absolute temperature and Cv is the classical speciﬁc heat parameter. Following the state relations (Equation (10)), the stress-like variables are obtained:

It is also worth noting that in the virgin (undamaged) state of the RVE, all these effective physical properties reduce to those of the virgin undamaged material. While at the ﬁnal fracture of the RVE these properties decrease as indicated by the state relations (Equations (26), (28), and (29)) and a strain-independent entropy as shown by Equation (27).

Constitutive Equations with Damage Remark 1. By using the Legendre transform a dual state potential (here it will be the free enthalpy) c*{r, X, R, Y, B, C, s} can be calculated from c(ee, a, r, d, b, c, T} by n

rc ðAi Þ ¼

Sup ee ;a;r;d;b;c;T

X

! Ai : ai rcðai Þ

ð33Þ

3.06.3.3

Dissipation Potentials and Complementary Relations

Following the selected dissipative phenomena, the residual inequality (see Equation (11)) can be witten as Y ¼ r : Dp X : a’ R’r þ Y : d’ C : c’ Bb’ q g Z0 T

ð34Þ

This inequality should be identically veriﬁed for each selected dissipative mechanism. The stress-like variables (i.e., r, X, R, Y, B, C) are given by the state relations (Equations (26)– (32)), the ﬂuxes or rate variables should be deﬁned by using the generalized standard materials framework as discussed above (see Equation (14)). This will be achieved by introducing both yield functions and dissipation potentials for each class of dissipative phenomena. As a ﬁrst approximation and similar to the uncoupling hypothesis between mechanical and thermal dissipations (see Equations (12) and (13)), the total dissipation (Equation (34)) will be additively decomposed into three terms, namely, plastic dissipation, damage dissipation, and thermal dissipation each of which is supposed to be independently positive or zero: Yp ¼ r : Dp X : a’ R’rZ0

ð35aÞ

’ Yd ¼ Y : d’ C : c’ BbZ0

ð35bÞ

q Yth ¼ g Z0 T

ð35cÞ

Each of these dissipations will be analyzed to derive the ﬂux variables deﬁning the constitutive equation associated with each selected dissipative phenomenon. 3.06.3.3.1

ratic scalar function of the force gi with T acting as a parameter: jðg; TÞ ¼ 12g k g

Thermal dissipation analysis: fully coupled heat ﬂow equation

Usually the heat equation is derived from the Fourier dissipation potential, which is a quad-

ð36Þ

In this time sensitive dissipation case, the heat ﬂux vector is given by

i

This allows the direct calculation of the strainlike variables ai Afee ; a; r; d; b; c; Tg as a function of the stress-like variables AiA{r, X, R, Y, B, C, s}. This transformation is independently and successively applied to each state variable.

329

q ¼ k:g T

or q ¼ kðTÞ g

ð37Þ

which is nothing but the classical linear Fourier’s heat ﬂow equation. The generalized heat equation can be derived from the equation of energy conservation or ﬁrst principle of thermodynamics (Lemaitre and Chaboche, 1990). In the present case it can be formally written in the following form: divðqÞ ¼ rCv T’ Yin w @r e @X @R : e’ þ : a’ þ r’ T @T @T @T @Y ’ @C @B ’ :dþ : c’ þ b ð38Þ @T @T @T

where w is the quantity of internal (volumetric) source of heat, whereas Yin ¼ Yp þ Yd is the intrinsic dissipation deﬁned by Equations (35a) and (35b). Associated with Fourier’s linear heat conduction (Equation (37)), the generalized heat equation (Equation (38)) can be used to derive an appropriate variational formulation to solve the associated IBVP by the FEA (Zienkiewicz et al., 1981). It is easy to show that the thermal dissipation inequality (Equation (35c)) is unconditionally veriﬁed as long as the thermal conductivity tensor k is positive deﬁnite. This ensures the thermodynamic admissibility of the heat ﬂow constitutive equation (Equation (37)). 3.06.3.3.2

Mechanical dissipation: two yield surface formulation

Both plastic ﬂow and damage ﬂow are supposed to be rate-independent mechanisms. Using the nonassociative framework both phenomena will be modeled using a yield function together with dissipation potentials. Concerning the rate-independent plasticity coupled with damage, a yield function in the stress space fp(r, X, R; d, T) and a plastic potential Fp(r, X, R; d, T) are introduced in order to derive the constitutive equations for plasticity with damage effect: ˜ R˜ sp o0 fp ¼ jjr* Xjj s

ð39Þ

1 1 b ˜2 R Fp ¼ fp þ aX˜ : C1 : X˜ þ 2 2Q

ð40Þ

330

Computational Damage Mechanics: Application to Metal Forming Simulation

where the temperature-dependent material constants a and b are the non-linearity parameters for kinematic and isotropic hardening, respectively. The parameter sp is the initial size of the plastic yield surface. The damage ﬂow is also assumed to be a rate-independent mechanical process. In a similar way to the arguments leading to plastic dissipation, it is assumed that the damage evolution equations are derived from a damage yield function in the Y space fd(Y, C, B; d, T) and a damage potential Fd(Y, C, B; d, T): fd ¼ jjY Cjjd B sd o0

ð41Þ

1 1 g2 2 B Fd ¼ fd þ g1 C : L1 : C þ 2 2G

ð42Þ

where the temperature-dependent material parameters g1 and g2 characterize the nonlinearity of the ‘‘kinematic’’ (g1) and ‘‘isotropic’’ (g2) motions of the damage yield surface which has sd as the initialsize. The notation r* X˜ s deﬁnes, in the effective stress space, the following norm of the effective equivalent stress: qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ˜ ¼ jjr* Xj r* X˜ : H : r* X˜ s qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ ðr X Þ : H˜ : ðr X Þ

ð44Þ

It tends to zero as the damage approaches its critical value leading to a zero equivalent stress inside the fully (isotropic) damaged RVE. Remark 2. When operators Mijkl and Nijkl are different, the equivalent stress (Equation (43)) becomes ˜ ¼ jjr* Xjj s

jjY Cjjs ¼

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ r M : N1 : X : H˜ : r M : N1 : X

with H˜ being unchanged and given by Equation (44). Similarly, for the damage ﬂow, the norm jjY Cjjd deﬁnes, in the damage conjugate force space, the equivalent damage energy

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ðY CÞ : J : ðY CÞ

ð45Þ

where the fourth-order symmetric and positive deﬁnite operator Jijkl deﬁnes the anisotropy of the damage evolution. Its explicit form can be simpliﬁed as proposed in the literature (Lee et al., 1985; Chow and Lu, 1989). The evolution equations for both plasticity and damage are derived from the generalized normality rule: for plastic ﬂow with hardening and damage effect: @Fp ’ @fp ’ Dp ¼ l’ p ¼ lp ¼ lp n @r @r

r’ ¼ l’ p

ð46Þ

@Fp ¼ Dp al’ p a @X

ð47Þ

l’ p @Fp ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃð1 b˜rÞ @R 1 jjdjj

ð48Þ

a’ ¼ l’ p

for damage ﬂow: @Fd ’ @fd ’ d’ ¼ l’ d ¼ ld ¼ ld m @Y @Y c’ ¼ l’ d

ð43Þ

where, for the sake of simplicity, we have supposed that the damage effect on both Cauchy stress r and kinematic internal stress X is the same, which requires Mijkl ¼ Nijkl (see Equations (19) and (20)). The fourth-order symmetric and positive deﬁnite operator Hijkl deﬁnes the anisotropy of plastic ﬂow for the undamaged medium. The explicit choice of this operator deﬁnes a quadratic or a nonquadratic yield criterion (Hill, 1971; Mellor, 1977). For the damaged RVE this operator transforms to H˜ ¼ MT : H : M1

release rate is given by

@Fd ’ ¼ d g1 l’ d c @C

@Fd ’ b’ ¼ l’ d ¼ ld ð1 g2 bÞ @B

ð49Þ ð50Þ ð51Þ

where the plastic and the damage multipliers satisfy the plastic and the damage loading/ unloading rules (Kuhn–Tucker condition): fp r0; l’ p Z0; l’ p fp ¼ 0 fd r0; l’ d Z0; l’ d fd ¼ 0

ð52Þ

The second-order tensors n and m represent the outward unit normal to the plastic and damage yield surfaces, respectively: n¼

˜ @fp @Fp MT : H : ðr* XÞ ¼ ¼ ˜ @r @r jjr* Xjj pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s with jjnjj ¼ n : H1 : n ¼ 1

ð53Þ

@fd @Fd J : ðY CÞ ¼ ¼ @Y @Y jjY Cjjd pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ with jjmjj ¼ m : J1 : m ¼ 1

ð54Þ

m¼

Remark 3. From Equations (46), (49), (53), and (54), one can observe that the normal to fp (fd) is the same as the normal to Fp (Fd). This ensures that the obtained model is associative only in the r and Y spaces but not in the internal variable spaces. An equivalent plastic strain rate p’ can be deﬁned in the strain space using the following

Constitutive Equations with Damage norm (see Equation (46)): pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ p’ ¼ Dp : H1 : Dp ¼ l’ 2p n : H1 : n ¼ l’ p

ð55Þ

This indicates that the isotropic hardening strain r is not equal to the cumulative plastic strain p unless the hardening is linear (b ¼ 0) and damage is zero (see Equation (48)). Similar to the equivalent plastic strain rate, the equivalent damage rate q’ can be introduced: q’ ¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ d’ : J1 : d’ ¼ l’ 2 m : H1 : m ¼ l’ d d

ð56Þ

which indicates that q’ is different from b’ unless the evolution of b is linear, i.e., g2 ¼ 0, as shown by Equation (51). Now, both the plastic and the damage multipliers will be calculated from the classical consistency conditions applied to the yield functions fp and fd. Three situations have to be distinguished according to whether one of the two or both yield conditions are satisﬁed. Mixed damage and plastic yielding. fp ¼ 0; f’p ¼ 0 for plastic yielding and fd ¼ 0; f’d ¼ 0 for damage ﬂow. In this case both the plastic and the damage multipliers are positive (l’ d 40 and l’ p 40). They are obtained by solving the following system: * : D Hpp l’ p Hpd l’ d þ HpT T’ ¼ 0 n:K * @K @ k˜ ðT T0 Þ :D m : ee : @d @d Hdp l’ p Hdd l’ d þ HdT T’ ¼ 0

ð57Þ

Plastic yielding without damage evolution. fp ¼ 0; f’p ¼ 0 for plastic yielding and fd o0 (no damage yielding) or fd ¼ 0; f’d o0 (damage unloading). In this case the damage multiplier is zero and only the plastic multiplier is positive and is given by l’ d ¼ 0 1 * : D þ Hpt T’ n:K l’ p ¼ Hpp

ð58Þ

where the McCauley brackets hX i deﬁne the positive part of X. The six scalar quantities which have the meaning of tangent ‘‘hardening’’ moduli—Hpp, Hpd, HpT, Hdd, Hdp, HdT—are deﬁned by: * þ C˜ : n aX˜ : n þ bR˜ 40 ð60Þ Hpp ¼ Q þ n : K T * : @M : MT : ee ðT T0 Þk : @M Hpd ¼ n : 2K @d @d @M :X :m þ 2M1 : @d " 1 @jjdjj R˜ : 2 ð1 jjdjjÞ @d # ˜ ðr XÞ : @ H=@d : ðr XÞ :m ð61Þ ˜ jjr Xjj s

@K e 1 @C : e C˜ : : X : MT HpT ¼ n : M : @T @T @k : MT k˜ ðT T0 Þ @T @M 1 : K : ee C˜ : C : X : MT þn: 2 @T @MT @Q=@T @sp ð62Þ R˜ ðT T0 Þk : @T @T Q * 1 @2K Hdd ¼ G þ m : L : m m : ee : 2 : ee 2 @d 2 ˜ 1 1 @ C 1 þ X : C˜ : 2 : C˜ : X 2 @d 1 R2 @ 2 jjdjj @ 2 k˜ e :m ðT T Þ : e 0 ˜ jjdjjÞ @d 2 2 Qð1 @d 2 ð63Þ þ g1 C : m þ g2 B 40

Hdp ¼ m :

( * ˜ @K @ k˜ 1 @ C : ee X : C˜ : ðT T0 Þk : :n @d @d @d

˜ 1 @ C 1 : C˜ : X þ aX : C˜ : @d ) ˜ @jjdjj RðQ bRÞ þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Q˜ 1 jjdjj @d

ð64Þ

Damage yielding without plasticity ( fd ¼ 0; f’d ¼ 0 for damage yielding, fp o0 for no plastic yielding, and fp ¼ 0; f’p o0 for elastic unloading). In this case the plastic multiplier is zero and the damage multiplier is positive and is given by l’ p ¼ 0

* 1 @K @ k˜ m : ee : l’ d ¼ ðT T0 Þ :D Hdd @d @d

þ HdT T’

331

ð59Þ

* 1 e @2K 1 @ 2 C˜ 1 1 : ee þ X : C˜ : e : C˜ : X 2 @d@T 2 @d@T 1 @Q=@T R2 @jjdjj @L : L1 : C 2 Q˜ @T Q˜ @d

@ 2 k˜ @k m : ðT T0 Þ : þ : ee @d@T @d @G=@T @sd B ð65Þ G @T

HdT ¼ m

Remark 4. Some bulk metal forming is performed under high-temperature conditions (hot forging for example). Under these warm

332

Computational Damage Mechanics: Application to Metal Forming Simulation

forming conditions some (superplastic) materials exhibit a certain rate effect for the inelastic deformation (viscoplasticity). However, the damage ﬂow is still rate independent, according to the fact that the process is done rapidly without the occurrence of ‘‘creep’’ damage. Accordingly, the damage multiplier, l’ d 40 is still given by the aforementioned consistency condition applied to fd; whereas the ‘‘viscon plastic’’ multiplier is given by l’ vp ¼ fp =K where K and n are viscosity parameters. The inelastic constitutive equations (Equations (46)–(48)) are still the same with the plastic multiplier being replaced by the viscoplastic one. The continuous elastoplastic tangent operator Cep (usually fourth-order nonsymmetric tensor) and its thermal contribution C Tp (second-order symmetric tensor) are derived from the time derivative of the Cauchy stress tensor (Equation (26)). By considering the additivity of the Eulerian strain rates (Equation (8)) together with the constitutive equations (Equations (46)–(51)) and the plastic and damage multipliers (Equations (57)–(65)), a rather lengthy derivation gives s’ ¼ Cep : D þ C Tp T’

C

simultaneous plastic and damage yielding, i.e., l’ p 4l’ d 40: ep

* # K * :n n:K * ¼K Hpp

ð69Þ

HpT * @K @M þ 2K : : MT : ee K:nþ M: @T @T Hpp @k @MT : k k˜ ð70Þ ðT T0 Þ : MT þ @T @T

C Tp ¼

*

damage yielding and plastic unloading, i.e., l’ d 40 and l’ p ¼ 0:

* ˜ * þ 1 m : @ K : ee ðT T0 Þ@ k Cep ¼ K Hdd @d @d

@M @MT # 2K : : MT : ee ðT T0 Þk : :m @d @d ð71Þ

C Tp ¼

HdT @M @MT : MT : ee ðT T0 Þk : :m 2K : @d Hdd @d @K @M þ 2K : : MT : ee þ M: @T @T @k @MT ðT T0 Þ : k k˜ : MT þ @T @T ð72Þ

ð66Þ

where the tangent operators are given according to the following four different situations: *

C

ep

*

mutual plastic and damage unloading, i.e., l’ d ¼ 0 and l’ p ¼ 0: * Cep ¼ K

ð73Þ

* :n * # K * :n ˜ * n:K Hpd K * * : n Hpp ee : @ K ðT T0 Þ@ k : m ¼K # Hpd K @d @d Hpp ðHpp Hdd Hpd Hdp Þ Hpp " #

T T e * : @M=@d : M : e ðT T0 Þ@M =d : k : m * 2K @K @ k˜ * :n þ # Hpp ee : ðT T0 Þ : m Hpd K Hpp Hdd Hpd Hdp @d @d ð67Þ

*

@K @M C Tp ¼ M : þ 2K : : M T : ee @T @T @k @MT T : k k˜ ðT T0 Þ :M þ @T @T

Hpd ðHpd HdT Hdd HpT Þ * K:n Hpp ðHpp Hdd Hpd Hdp Þ @K @M þ M: þ 2K : : M T : ee @T @T Hpp HdT Hdp HpT @M : M T : ee 2K : þ @d Hpp Hdd Hpd Hdp @MT ðT T0 Þ :k :m @d @k @MT : MT þ ðT T0 Þ : k k˜ ð68Þ @T @T

The same kind of calculations can be applied to the other internal stresses (X, R, Y, C, B) in order to derive similar tangent moduli. This allows us to express the rates of the stress-like variables in terms of the total strain rate D and the temperature rate T’ (Saanouni et al., 1994).

plastic yielding and damage unloading, i.e., l’ p 40 and l’ d ¼ 0:

Remark 5. In the case of the two-surface formulation, the kinetics of the damage description is often supposed to have only an

C Tp ¼

ð74Þ

Constitutive Equations with Damage increase in its radius described by (b, B) without any translation of its center, i.e., (c ¼ 0, C ¼ 0). This has been widely used in the literature (Lee et al., 1985; Chow and Lu, 1989; Zhu and Cescotto, 1995; Voyiadjis and Kattan, 1999; Hammi, 2000, among others). To check the thermodynamic admissibility of the above developed two-surface constitutive equations, the inequalities (Equations (35a) and (35b)) should be separately satisﬁed. By using the constitutive equations (Equations (46)–(51)) these two inequalities become, after some simple algebraic calculations, b 1 Yp ¼ l’ p fp þ aX : C˜ : X þ R2 þ sp Z0 Q˜ h i g 2 1 d Y ¼ l’ d fd þ g1 C : L : C þ B2 þ sd Z0 G

ð75Þ

These two inequalities are satisﬁed as long as ˜ g2 =G; sd ; sp all these parameters—a; g1 ; b=Q; are positive or null. 3.06.3.3.3

Mechanical dissipation: one yield surface formulation

In many metal forming processes of monophase ductile metallic materials, the damage initiation cannot take place without signiﬁcant plastic deformation. Consequently, the most common coupled ductile plastic-damage models are based on a single yield function together with a single plastic potential to describe both plastic and damage ﬂows. This can be obtained from the above developed two-surface model (with c ¼ C ¼ 0, b ¼ B ¼ 0 and fd ¼ Fd ¼ 0) by using only fp(r, X, R; d, T) and a plastic potential Fpd(r, X, R, Y; d, T). The coupled plastic yield function is still given by Equation (39), while the plastic potential is modiﬁed to include a new term Fd(Y; d, T,y), a convex function of the damage force Y and containing other internal variables as parameters: 1 1 b ˜2 R þ Fd ðY; d; T; :::Þ Fpd ¼ fp þ aX˜ : C1 : X˜ þ 2 2Q ð76Þ

where the choice of the explicit form of Fd, in this anisotropic case, is not a trivial task. A possible choice is given by Hammi (2000) which allows us to recover the Lemaitre (1992) potential in the fully isotropic case as (see later) Fd ðY; d; TÞ ¼

b Y : J : Y S T M :: M 2 S

ð77Þ

The ﬂux variables—Dp ; a’ ; r’—are still given by Equations (46)–(48) with l’ p being replaced by

333

’ while the damage rate is now given by l; b J : Y @Fpd ’ T ¼ l M :: M d’ ¼ l’ S @Y

ð78Þ

where the unique plastic-damage multiplier, in the case of plastic ﬂow, is given by 1 * : D þ HT T’ l’ ¼ n:K H

ð79Þ

where the scalars H and HT are now given by * þ C˜ : n þ Q aX : n þ bR˜ H ¼n : K ( ˜ @Fd 1 ðr XÞ : @ H=@d : ðr XÞ : þ ˜ 2 @Y jjr* Xjj s T @M e :K :e þn: 2M : @d

@M @MT : C : a ðT T0 Þ :k 2MT : @d @d ) 1 @jjdjj R˜ 40 ð80Þ þ 2 ð1 jjdjjÞ @d @M @K : K þ MT : : M : ee HT ¼ n : 2MT : @T @T @M T T @C þ n : 2M : :CþM : :M :a @T @T @k @MT þ : k k˜ ðT T0 Þn : MT : @T @T R @Q @sp ð81Þ þ Q @T @T

From Equation (26), the time derivative of the Cauchy stress can be written as @M * : D l’ K * : n þ @Fd : : K : ee 2MT : s’ ¼ K @d @Y @M :k ðT T0 Þ @d @M @K þ T’ 2MT : : K þ MT : : M : ee @T @T T @k @M þ : k k˜ ð82Þ ðT T0 Þ MT : @T @T

This can be rewritten in the same form as Equation (66), from which the continuous tangent moduli can be obtained with the help of Equation (79) and assuming the mechanical and thermal loading condition: * * * * n:K # K:n K:n Cep ¼ K H H @Fd T @M : K : ee : 2M : # @d @Y T @M :k ðT T0 Þ @d

ð83Þ

334

Computational Damage Mechanics: Application to Metal Forming Simulation

and

ratio n as

HT * @Fd @M K:nþ : 2MT : CT ¼ : K : ee H @Y @d @M :k ðT T0 Þ @d @M @K : K þ MT : : M : ee þ 2MT : @T @T T @k @M : k k˜ ðT T0 Þ MT : þ @T @T

nE 3Ke 2Ge ¼ ð1 þ nÞð1 2nÞ 3 E Ge ¼ m ¼ 2ð1 þ nÞ E 3le þ 2m Ke ¼ ¼ 3ð1 2nÞ 3 le ¼

ð84Þ

It is worth noting that the same type of tangent moduli can be obtained for the internal stresses X, Y, and R. Alternatively, it is easy to check that, when damage vanishes, the classical nonassociative thermo-elastoplastic constitutive equations are recovered. Remark 6. In this formulation, the nonsymmetry of the elastoplastic tangent operator Cep is clearly observed from Equation (67) for the two-surface formulation, and from Equation (83) for the single-surface formulation. The symmetry of these operators is recovered when damage vanishes, i.e., uncoupled formulation. Remark 7. The present formulation can be made in the strain-like variable spaces instead of the stress-like variable spaces. To this end, the yield functions (Equations (39) and (41)) as well as the dissipation potentials (Equations (40), (42), (76), and (77)) should be transformed by the Legendre–Fenchel transformation in the corresponding dual spaces.

3.06.3.4

The damage effect tensors (see Equations (19) and (20)) become M¼N¼

K ¼ 2m1 þ le 1#1 ¼ 2Ge 1d þ Ke 1#1 C ¼ 231d ; H ¼ 321d ; k ¼ a1

ð85Þ

(ii) For the damage properties L ¼ 231d ;

J ¼ 321d

ð86Þ

where 1d ¼ 1 ð1=3Þ1#1 is the deviatoric fourth-order unit tensor and the isotropic elastic properties (le, Ke, Ge, m) are expressed in terms of Young’s modulus E and Poisson’s

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 D 1 and jjdjj ¼ D

ð88Þ

Accordingly, the overall constitutive equations in this fully isotropic case are summarized below for both two-surface and single-surface formulations. 3.06.3.4.1

Two-surface formulation

In this fully isotropic two-surface formulation, the state relations (Equations (26)–(36)) transform to r ¼ l* e ðee : 1Þ1 þ 2me * e ð3le þ 2mÞ*aðT T0 Þ1 3K˜ e 2G˜ e e * T0 Þ1 ðe : 1Þ1 þ 2G˜ e ee 3Ke aðT ¼ 3 ð89Þ

The Fully Isotropic Case

In many situations, it is sufﬁcient to assume the initial isotropy of the thermoelastic behavior and the full isotropy of the plastic, hardening, and damage yielding (represented by the scalars D and Y). In this case, the physical properties of the undamaged RVEs are given in the following. (i) For the thermo-elastoplastic properties

ð87Þ

1 Cv s ¼ 3Ke a* ðee : 1Þ þ ðT T0 Þ r T0

ð90Þ

˜ X ¼ 23Ca

ð91aÞ

˜ R ¼ Qr

ð91bÞ

Y ¼ Ye þ Yk þ Yi

ð92aÞ

Ye ¼ 12le ðee : 1Þ2 þ mðee : ee Þ a 32Ke pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃðT T0 Þðee : 1Þ 1D "

# jjrjj * 2s 2 sH 2 ð1 þ nÞ þ 3ð1 2nÞ ¼ jjrjj * s 2E˜ 3 þ

1 a ðT T0 Þðr* : 1Þ 2 ð1 DÞ

˜ 2 1 1 jjXjj s Yk ¼ Ca : a ¼ ; 3 2 C˜

1 1 R˜ 2 Yi ¼ Qr2 ¼ 2 2 Q˜

ð92bÞ

ð92cÞ

C ¼ 23Lc

ð92dÞ

B ¼ Gb

ð92eÞ

Constitutive Equations with Damage where all the elastic properties of the damaged material are indicated by (‘‘tilde’’) and are expressed with the ‘‘damaged’’ Young’s modulus E˜ ¼ ð1 DÞE: Similarly, the kinematic hardening modulus, the isotropic hardening modulus, and the thermal conductivity for the damaged RVEs p are C˜ ¼ ð1 DÞC; Q˜ ¼ ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð1 DÞQ; and a* ¼ 1 Da; respectively. The norm of stress-like variable for any secondorder symmetric tensor t is deﬁned in the stress space by qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ jjtjjs ¼ 32td : td

ð93Þ

where td ¼ t tH 1 is the deviatoric stress-like tensor, while tH ¼ ð1=3ÞtrðtÞ is the hydrostatic stress. Concerning the time-independent plastic ﬂow, the same yield function and dissipation potential (Equations (39) and (40)) are used. Then, the complementary relations are still given by the same equations (Equations (46)– (48)), in which the normal to the plastic yield surface n should be changed to 3 1 ðrd X d Þ n ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1 D jjr Xjjs

ð94Þ

The equivalent plastic strain rate (Equation (55)) is now given by p’ ¼

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ p 2 p 3D : D ¼

l’ p 1D

ð95Þ

while, for the time-independent damage ﬂow, the yield surface fd (Equation (41)) and the damage potential Fd (Equation (42)) need to be changed according to the fact that the damage force Y is now a scalar quantity. Three possibilities can be considered. The ﬁrst and the simplest one assumes that there is no translation of the damage yield surface in the stress (resp. strain) spaces leading to c ¼ 0, C ¼ 0. These functionals are then given by fd ¼ Y B sd o0 S Y sþ1 1 1 g2 2 B B þ Fd ¼ f sþ1 S 2G ð1 DÞ

ð96Þ

ð97Þ

where now the potential Fd is a power function of the variable Y, with S, s, and f are new material parameters depending on the temperature. The second possibility is to save the internal variables (c, C) and to express these functionals in the stress space (or similarly in the strain space): fd ¼ jjr Cjjd B sd o0

ð98Þ

335

1 1 g2 2 Fd ¼ fd þ g1 C : L1 : C þ B 2 2G

ð99Þ

In this case, the damage constitutive equations are still given by Equations (49)–(51) in which the normal to the damage yield surface is now given by m¼

3 ðrd Cd Þ 2 jjr Cjjs

ð100Þ

In this study, for the sake of simplicity in this isotropic case, only the ﬁrst possibility (i.e., Equations (96) and (97)) will be used to derive the damage constitutive equations. Accordingly, the damage constitutive equations (Equations (41) and (51)) should be replaced by @Fd ’ Y s 1 ¼ ld ¼ l’ d Yu D’ ¼ l’ d S ð1 DÞf @Y @Fd ’ ¼ ld ð1 g2 bÞ b’ ¼ l’ d @B

ð101aÞ

ð101bÞ

The plastic and the damage multipliers are calculated as the fully anisotropic case discussed above. They are given by considering the following three cases. Mixed damage and plastic yielding. In this case both the plastic and the damage multipliers are positive (l’ d 40 and l’ p 40) and obtained by solving the following system: d

3m*

ðr* d X˜ Þ : D Hpp l’ p Hpd l’ d þ HpT T’ ¼ 0 jjr Xjjs ð102aÞ 1 D r* þ að3l þ 2mÞðT T0 Þ1 : pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1D Hdp l’ p Hdd l’ d þ HdT T’ ¼ 0

ð102bÞ

Plastic yielding without damage. In this case the damage multiplier is zero and only the plastic multiplier is positive and can be written as *

l’ d ¼ 0

+ d ˜d ’lp ¼ 1 3m* ðr* X Þ : D þ HpT T’ Hpp jjr Xjjs

ð103Þ

Damage yielding without plasticity. In this case the plastic multiplier is zero and the damage multiplier is positive and is given by l’ p ¼ 0

1 1 D l’ d ¼ ½r* þ að3l þ 2mÞ1 : pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ þ HdT T’ Hdd 2 1D ð104Þ

336

Computational Damage Mechanics: Application to Metal Forming Simulation

The six tangent ‘‘hardening’’ moduli—Hpp, Hpd, HpT, Hdd, Hdp, HdT—are now deﬁned as 3 ðrd X d Þ ˜ Hpp ¼ 3m þ C þ Q a : X þ bR˜ 40 2 jjr Xjjs ð105Þ

Hpd ¼

Yu h sp i jjr Xjjs R˜ 1D 2

@m @C=@T X HpT ¼ n : 2ð1 DÞ ee @T C @Q=@T ˜ @sp R @T Q

ð106Þ

ð107Þ

where the normal n is given by Equation (94). Yu

1 að3l þ 2mÞ ð1 DÞ3=2 4 ðT T0 Þ ðee : 1Þ gB

Hdd ¼ G þ

ð108Þ

* Cep ¼ l1#1 þ 2m1 *

að3l þ 2mÞ Yu r* þ DT 1 2 ð1 DÞHpp að3l þ 2mÞ 1 # r* þ DT 2 d 2m* ð3=2Þ ðr* d X˜ Þ=jjr Xjjs Hpp Hdd Hpd Hdp " d ðr* d X˜ Þ Hpd # 3mH pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ * dd jjr Xjjs 1D !# 3 r* þ DT ð3l þ 2mÞ1 2 þ

ð111Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @m e @l e ðe : 1Þ1 1D 2 e þ C Tp ¼ 1 D @T @T @l @m @a DT a 3 þ2 þ ð3l þ 2mÞ @T @T @T 1 að3l þ 2mÞ 1 3m* Hpp d Hpp HdT Hdp HpT ðr* d X˜ Þ HpT þ Hpp Hdd Hpd Hdp jjr Xjjs að3l þ 2mÞ Yu HdT 1 ð112Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ r* þ DT 2 1 D Hdd

In the case of plastic yielding with damage unloading, i.e., l’ p 40 and l’ d ¼ 0; d 3 1 ðr* d X˜ Þ a ˜ 2 : r* þ jjXjj Hdp ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ s 2 1 D jjr Xjjs C b 2 ˜ þ R˜ R40 Q

* Cep ¼ l1#1 þ 2m1 * ð109Þ

! ! d d 1 ðr* d X˜ Þ ðr* d X˜ Þ # 3m* 3m* Hpp jjr Xjjs jjr Xjjs ð113Þ

1 @l e að3l þ 2mÞ ðe : 1Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ HdT ¼ mðee : ee Þ þ 2 @T 1D T T0 @l pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ a 2 þ 3 @T 1D @a ðee : 1Þ þ ð3l þ 2mÞ @T

@C=@T jjXjjs 2 @Q=@T R 2 þ þ 2 2 Q˜ C˜ þ

@G=@T @sd B @T G

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @m e @l e ðe : 1Þ1 1D 2 e þ C Tp ¼ 1 D @T @T @l @m @a þ2 þ ð3l þ 2mÞ DT a 3 @T @T @T d HpT ðr* d X˜ Þ að3l þ 2mÞ 1 3m* ð114Þ Hpp jjr Xjjs ð110Þ

The time derivative of the Cauchy stress (see Equation (89)) can be written in the same form as Equation (66), in which the tangent elastoplastic and thermal operators in the case of simultaneous plastic and damage yielding conditions are

While for the case of damage yielding with elastic unloading, i.e., l’ d 40 and l’ p ¼ 0; * Cep ¼ l1#1 þ 2m1 *

að3l þ 2mÞ Yu 1 r* þ DT 2 Hdd ð1 DÞ að3l þ 2mÞ 1 # r* þ DT 2

ð115Þ

Constitutive Equations with Damage

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @m e @l e C Tp ¼ 1 D 1D 2 e þ ðe : 1Þ1 @T @T @l @m @a þ2 þ ð3l þ 2mÞ DT a 3 @T @T @T að3l þ 2mÞ 1 að3l þ 2mÞ Yu HdT pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 ð116Þ r* þ DT 2 1 D Hdd

Finally, when there is a plastic unloading together with damage unloading, i.e., l’ d ¼ 0 and l’ p ¼ 0; ep

C

* ¼ l1#1 þ 2m1 *

ð117Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @m e @l e ðe : 1Þ1 1D 2 e þ C Tp ¼ 1 D @T @T @l @m @a þ2 þ ð3l þ 2mÞ DT a 3 @T @T @T að3l þ 2mÞ 1 ð118Þ

where the notation DT ¼ T – T0 has been used.

3.06.3.4.2

One-surface formulation

Similar to the fully anisotropic one-surface formulation (see, e.g., Section 3.06.3.2), it is often very helpful to use one-surface formulation in the case of the isotropic ﬂow. The simplest model can be obtained from the above developed two-surface formulation by neglecting the internal variables associated to the damage except the pair of scalar variables (D,Y) (Saanouni et al., 1994). In this case, the state relations are still unchanged and given by Equations (89)–(92c), and the complementary relations are given by Equations (46)–(48) with the normal deﬁned as Equation (92a). The damage evolution equation is given by Equation (101a) if the coupled plastic and damage potential is chosen as 3 a ˜ ˜ 1 b ˜2 S Y sþ1 1 F ¼f þ X:Xþ R þ 4C 2Q sþ1 S ð1 DÞf ð119Þ

where f is the coupled plastic-damage yield function deﬁned by Equation (39) with the same stress norm as Equation (93). The plastic multiplier is deduced from Equation (102a) in which the damage multiplier is zero. In the case of plastic yielding, this gives * + d ˜d ’l ¼ 1 3m* ðr* X Þ : D þ HT T’ H jjr Xjjs

with the plastic tangent modulus H:

ð120Þ

337

sp H ¼ 3m þ C þ Q þ Yu 2ð1 DÞ 3 ðrd X d Þ ˜ ˜ : X þ bR 40 a 2 jjr Xjjs

ð121Þ

and the thermal term HT is given by Equation (107) above. The elastoplastic tangent operators used in the time derivative of the Cauchy stress tensor (Equation (66)), in the case of this single yield surface formulation, are C

ep

d 1 ðr* d X˜ Þ * ¼ l1#1 þ 2m1 3m* * H jjr Xjjs ! d ðr* d X˜ Þ # 3m* jjr Xjjs ! d Yu ðr* d X˜ Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3m* jjr Xjjs H 1D að3l þ 2mÞ 1 # r* þ ðT T0 Þ 2

!

ð122Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @m e @l e ðe : 1Þ1 C Tp ¼ 1 D 1D 2 e þ @T @T

@l @m þ2 ðT T0 Þ a 3 @T @T @a ð3l þ 2mÞ að3l þ 2mÞ 1 þ @T d HT ðr* d X˜ Þ 3m* H jjr Xjjs að3l þ 2mÞ Yu HT 1 r* þ ðT T0 Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 1D H

ð123Þ

It is worth noting that when the damage vanishes, the classical elastoplastic constitutive equations are recovered as shown in Saanouni et al. (1994) and Hammi (2000). This closes the presentation of the fully coupled constitutive equations accounting for many thermomechanical phenomena as heat exchange, elastoplasticity, isotropic and kinematic hardening, and the ductile damage. All these phenomena can be considered separately as isotropic or anisotropic. The present formulation is valid for both ‘‘small-strain hypothesis’’ and ‘‘large plastic strain but small elastic strain hypothesis’’ as discussed in the Section 3.06.2.1. 3.06.3.5

On the Frictional Constitutive Equations

In metal forming processes, the contact and friction between the deformed part and tools play a prominent role. Three points should be considered in the numerical simulation of

338

Computational Damage Mechanics: Application to Metal Forming Simulation

contact with friction in metal forming: the contact-friction constitutive equations, the associated solution strategies, and the contact search algorithm. This section is devoted to a brief review of the principal models used to describe the contact-friction in metal forming processes. Generally, the contact-friction constitutive equations describe small-sliding or ﬁnite-sliding contact with arbitrary rotation of the surfaces of two bodies with respect to each other. They are strongly dependent on many factors affecting the forming process such as: the contact pressure, the sliding speed, the material of the contacting bodies, surface roughness, lubrication, temperature, etc. This large variety of inﬂuencing parameters render the development of these contact-friction constitutive equations a nontrivial task. The widely used approach in this ﬁeld is still the phenomenological method similar to the timeindependent and time-dependent plasticity framework (see Section 3.06.3). A friction yield function is introduced, in the stress space, to derive a relation between the normal stress rn, the tangential stress st, and the relative tangential velocity vt at the contact point between the part (deforming plastically) and the tool (supposed rigid or elastic body). Two families of contact-friction laws can be derived in this way, namely: the isotropic and the anisotropic friction models. Here, we limit ourselves to the isotropic case (see Wagoner and Chenot, 1997 chapter 9).

3.06.3.5.1

Time-independent formulation

To be short, we restrict ourselves to the classical empirical presentation of these models which can be formally written as if teq otcrit ) vt ¼ 0 ðsticking contactÞ if teq ¼ tcrit ) (dZ0; vt ¼ dt

ð124Þ

ðsliding contactÞ

where teq ¼ ðt21 þ t22 Þ1=2 is the equivalent tangential stress. The difference between the main time-independent friction models comes from the deﬁnition of the critical tangential stress. For example, the Coulomb-type friction model is deﬁned by tcrit ¼ ZP ¼ Zsn giving if teq oZP ) vt ¼ 0 ðsticking contactÞ if teq ¼ ZP ) (dZ0; vt ¼ dt ðsliding contactÞ teq oZP ) impossible

ð125Þ

where Z is the temperature-dependent friction coefﬁcient. The Tresca-type friction model is deﬁned by tcrit ¼ sy =2; while in the Von-Mises-type friction modelpthe ﬃﬃﬃ critical shear stress is given by tcrit ¼ sy = 3; where sy is the yield stress in tension. 3.06.3.5.2

Time-dependent formulation

The frequently used time-dependent friction model is the Norton–Hoff-type model deﬁned by if teq otcrit ) vt ¼ 0 ðsticking contactÞ

sy vt p if teq Ztcrit ) t tcrit ¼ kpﬃﬃﬃ 3 V0 ðsliding contactÞ

ð126Þ

To these friction constitutive equations, one might add the unilateral contact or Signorini conditions (similar to the loading–unloading conditions in plasticity) deﬁned by un r0; Fn r0; un Fn ¼ 0

ð127Þ

where un and Fn are the normal components of the displacement and the force vectors expressed in a local orthogonal axis system. In Equation (127), the ﬁrst term ðun r0Þ is called the nonpenetration condition; the second ðFn r0Þ expresses the fact that at each contact point the normal force is negative; while the third one ðun Fn ¼ 0Þ expresses two situations: (i) when there is contact ðun ¼ 0Þ; the normal force is active ðFn r0Þ; and (ii) when there is no more contact ðun o0Þ; the normal force should be zero ðFn ¼ 0Þ: Note that all these standard contact-friction models are available in a standard fashion in most metal forming simulation codes (see Abaqus Theory Manual for example). Their numerical implementation needs some regularization to avoid some numerical difﬁculties related to the convergence rate and to the CPU time needed to solve the associated IBVP. Finally, let us note that, for the sake of simplicity, it is supposed here that the damage has no direct effect on the contact-friction constitutive equations. In fact, it is supposed that the damage cannot grow on the material points lying on the contact surfaces, since the stress state in such points is compressive. Accordingly, it has been suggested (Saanouni and Franqueville, 1999; Hammi, 2000; Saanouni and Hammi 2000) that the elastic part of the damage energy release rate, in the isotropic case, deﬁned by Equation (92b) is slightly modiﬁed

Computational Aspects by replacing the hydrostatic pressure sH by its positive part /sHS. This allows us to take into account the fact that ductile damage (growth of voids) cannot develop at the points where the hydrostatic pressure is compressive and the equivalent (shear) stress is very low. This gives a nonsymmetric criterion compared to the Lemaitre’s original one, as shown in Figure 2. 3.06.4

COMPUTATIONAL ASPECTS

This section is devoted to the formulation of the necessary equations that enable engineers to solve the elastoplastic-damage-contact-friction problem, thanks to the FEA. First, the associated IBVP will be formulated based on the updated Lagrangian description. 3.06.4.1

Formulation of the IBVP

Numerical simulation of metal forming processes involves the use of the FEA to compute, for each applied external load increment, the successive new body conﬁgurations satisfying the equilibrium equations as well as the constitutive equations accounting for the above discussed coupled thermomechanical phenomena. These are: heat exchange, large plasticity, nonlinear hardening, ductile damage, and contact with friction subject to some constraints as initial conditions, ﬁxed boundary as well as the evolving boundary conditions due to evolving contact between the deformed part and tools. This deﬁnes a highly nonlinear IBVP which needs to be incrementally linearized in order to ﬁnd the roots using either an iterative or noniterative implicit or explicit scheme. Practically, solving this kind of IBVP needs both the spatial discretization of the

339

body using the FEM and the time discretization using the classical ﬁnite difference method (FDM). Considerable work has been published in this ﬁeld to solve IBVPs with various types of constitutive equations applied to different kinds of nonlinear mechanical structures with the main objective to obtain a convergent solution at a minimum cost. Here, the list of references is intentionally limited to the general books concerned with the numerical aspects of nonlinear mechanics of solids and structures (e.g., Oden, 1972; Owen and Hinton, 1980; Bathe, 1981; Hughes, 1987; Kleiber, 1989; Zienkiewicz and Taylor, 1994; Bonet and Wood, 1997; Belytschko et al., 2000 for general aspects of nonlinear FEA). Particularly, we refer to the books by Crisﬁeld (1991, 1997) and Simo and Hughes (1998), especially dedicated to the small and large plasticity IBVPs. Also some references concerned with the numerical aspects of large plasticity coupled with damage related to metal forming will be cited in this section. Before discussing the space and time discretizations of the above deﬁned IBVP, let us deﬁne precisely the problem with concern in the particular case of purely mechanical phenomena (thermal aspects will be neglected for the sake of brevity). The deforming part is supposed to occupy, at time tA½0; tmax ; the volume V limited by the boundary G. The part is supposed to be submitted to an imposed: displacement ﬁeld u on the part Gu of G, force ﬁeld F on the other part GF of G, contact force ﬁeld Fc on the third part Gc of G, and volumetric force ﬁeld fv on the whole V,

140

Equivalent stress (MPa)

120 100 80 60 40 20 0 -200

-150

-100

-50

0

50

100

150

200

Hydrostatic pressure (MPa)

Figure 2 Ductile damage criteria: classical Lemaitre’s criterion (dashed line) and modiﬁed criterion (solid line).

340

Computational Damage Mechanics: Application to Metal Forming Simulation

where Gu ,GF ,Gc ¼ G and Gu -GF -Gc ¼ +: The displacement- (or velocity-) based FEA consists in ﬁnding the kinematically admissible (KA) displacement ﬁeld (u) satisfying, at each time t, the weak form of the principle of virtual power (neglecting the damping phenomena):

Z

þ

r : dD dV þ

V Z

Z

F d’u dG þ

GF

¼

Z

f v d’u dV

ZV

F c d’u dG

Gc

r¨u d’u dV ;

8d’u KA

ð128Þ

V

where d’u is the vector of virtual velocity (displacement), dD is the associated virtual strain rate, u¨ is the vector of the acceleration, and r is the body density. Here the Cauchy stress tensor r(x,t) should satisfy the local constitutive equations developed in Section 3.06.3, and hence it is a function of the displacement u(x,t) through the velocity gradient L (see Equations (2) and (3)) and the internal variables representing plasticity, hardening, and damage phenomena as discussed above. The IBVP deﬁned by Equation (128) is nonlinear with respect to both the geometry and the material. Its numerical solution, for a given material and loading conditions, involves both time and space discretizations as will be discussed hereafter. In order to obtain the numerical solution of the IBVP, it is necessary to linearize properly Equation (128). For this, two approaches are commonly used: (i) the continuous form (Equation (128)) is ﬁrst linearized using the concept of directional derivatives to obtain a weak form to be discretized using the FEM; and (ii) Equation (128) is ﬁrst discretized using FEA and then linearized with respect to the nodal positions. The latter approach will be used herein. Finally, we note that for the stability conditions and the uniqueness of the solution of this IBVP, the reader is referred to Simo and Hughes (1998), among many others. 3.06.4.2

Spatial Discretization of the IBVP

The domain of interest is discretized by the displacement-based FEM to arrive at the discrete counterpart of the weak form deﬁned by Equation (128). Accordingly, the total volume of the body V is discretized into a ﬁnite number of Ne elements, each Shaving e e the elementary volume V e such as V E N e¼1 V where xe indicates, formally, the mean size of S e the element (mesh size) and N e¼1 designates the addition operation over the all elements. For the sake of simplicity, only the isoparametric

kinematically admissible solid and contact elements will be considered. These class of C(0) FEs are based on a continuous interpolation of the displacement (or velocity) ﬁeld and lead to discontinuous interpolations of strain and stress ﬁelds over a typical element V e. Now, attention is focused on a typical solid element V e with a size xe over which, and according to the Galerkin method, the real and virtual displacement ﬁelds are linearly interpolated according to (here and throughout the context the classical FE matrix notations will be used) fue ðx; tÞg ¼ NIe ðxÞ ueI ðx; tÞ

ð129aÞ

fdue ðx; tÞg ¼ NIe ðxÞ dueI ðx; tÞ

ð129bÞ

where NIe denotes the matrix of the interpolation functions (Lagrange polynomials) for the element node number I. The successive time derivatives of Equation (129) give the velocity ﬁelds: fu’ e ðx; tÞg ¼ NIe ðxÞ u’ eI ðx; tÞ fdu’ e ðx; tÞg ¼ NIe ðxÞ du’ eI ðx; tÞ

ð130Þ

as well as the acceleration ﬁeld: fu¨ e ðx; tÞg ¼ NIe ðxÞ u¨ eI ðx; tÞ

ð131Þ

Note that according to the supposed isoparametry of the used FEs the coordinates of any material point within the element (e) are expressed in terms of the nodal coordinates of node I by the same interpolation function as the displacement ﬁeld: fX e ðx; tÞg ¼ NIe ðxÞ XIe ðx; tÞ

ð132Þ

The approximation of the displacement ﬁeld discussed above concerns only the solid element without contact. Various different approaches exist to take into account the contact (with or without friction) in the dynamic IBVP deﬁned by Equation (128) (see the special issue of the JMTA, 1988 dedicated to the contact problems). These can be classiﬁed into two families: (i) the nonregularized methods as the mathematical programming or direct methods, and (ii) the regularized methods which introduce the contact forces in the internal force vector. The use of a speciﬁc contact element is the most used approach among the regularized methods, to which we restrict our discussion in this work. Accordingly, for contact elements (1D or 2D elements) the local displacements at each point are expressed in a local triad by e uc ðx; tÞ ¼ ½Re T NIe ðxÞ ueI ðx; tÞ

ð133Þ

Computational Aspects

341

where Re is an appropriate orthogonal rotation matrix between the global coordinate system and the local triad at the contact node under consideration (Charlier and Habraken, 1990). By using these approximations together with Equation (128), one gets, after some algebraic manipulations, the discrete form of the elementary weak functional J(e) in the following matrix form:

e e Fext J e ¼ hdu’ e i ½M e fu¨ e g þ Fint

Z

r½N T ½N dV e

ð135Þ

is the element consistent mass matrix; Z

½Be T fse g dV e þ

Z

Gec

Ve

Bec

T e sc dGec ð136Þ

is the element internal force vector including both contributions of solid element and its associated contact element if there is one; and e Fext ¼

Z Ve

½N e T fve dV e þ

Z Ge

½N e T fF e g dGe ð137Þ

is the element external force vector containing body forces and surface traction. The matrices B and Bc are the strain–displacement matrix for the solid and the contact elements, respectively. They are deﬁned by e @N and Bec ¼ ½Re T ½N e ½B ¼ @Xi e

ð138Þ

fFint g ¼

J ¼ dU’

e¼1

(

"

fFext g ¼

Ne a

#

½M U¨ þ

e¼1

)!

e Fext

e

;

(

Ne a

e Fint

ð141bÞ

Ne a

ð141cÞ

e Fext

It is worth noting that Equation (140) is quite similar to the case of small-strain analysis, except the fact that all computations are performed with respect to the current deformed conﬁguration. Note that the derivation of element mass matrix (Equation (135)) and vectors (Equations (136) and (137)) involves the integration of the shape functions and/or their derivatives over the deformed element (V e), i.e., using the real coordinate system (x, y, z). It is well known in FEA that these integrals can be evaluated more easily if the shape functions are written in terms of a local (often called natural or nondimensional), coordinate system (x, Z, z) located on or within the boundaries of the undeformed reference element (V r). This requires the transformation of the problem coordinates (x, y, z) to a local coordinates by using appropriate Jacobian matrices. For the isoparametric elements used here, this coordinate transformation uses the same shape (or interpolation) functions N as those used for the displacement ﬁeld (Equation (129)). Accordingly, the above deﬁned matrix and vectors when written on the reference element (V r) transform to ½M ¼

Ne a

Ne a e¼1

e

By calling the standard FE assembly operator, the discrete counterpart of the weak form of the equilibrium equations can be written as e

ð141aÞ

e¼1

‘

Ne a

½M e

ð134Þ

Ve

e Fint ¼

Ne a

e¼1

in which ½M e ¼

½M ¼

Z

r½N T ½N Jv dV r

ð142Þ

Vr

e Fint ¼

) e Fint

Z

½Be T fse gJv dV r Z e T e þ Bc sc Jc dGrc Vr

ð143Þ

Grc

e¼1

8dU’ KA

ð139Þ

e¼1

where U¨ and dU’ are the global acceleration and virtual velocity vectors for the overall structure nodes. Since the nodal velocity variations dU’ are arbitrary, from Equation (139) one arrives at the discrete system of nonlinear differential equations: ½M U¨ þ fFint g fFext g ¼ 0

ð140Þ

where the global matrix and force vectors for the overall structure are given by

e ¼ Fext

Z

½N e T fve Jv dV r Z þ ½N e T fF e gJs dGr Vr

ð144Þ

Gr

where the superscript (r) refers to the reference element expressed in the local coordinate system ðx; Z; BÞ; whereas Jv, Js, and Jc are the determinants of the Jacobian matrices transforming V e to V r, Ge to Gr ; and Gec to Grc ; respectively. These integrals (Equations (142)– (144)) are generally evaluated numerically using the well-known Gauss quadrature method

342

Computational Damage Mechanics: Application to Metal Forming Simulation

which involves the knowledge of the integrated quantities (as the stress ﬁeld) only at some discrete points (or Gauss integration points) deﬁned in the reference element. Remark 8. Note that according to the ﬁnite transformation hypothesis, in Equations (142)– (144), the matrices Be, Bec and the determinants Jv, Js, Jc are strongly dependent on the displacement ﬁeld (geometrical nonlinearities). Similarly, the stress vectors se (plasticity) and sec (contact) are also displacement dependent (material nonlinearities). However, the shape function matrix Ne is independent of the real coordinates x, y, z (i.e., of the displacement ﬁeld) and only depends on the local coordinates ðx; Z; BÞ: 3.06.4.3

Incremental Solution Procedures

The highly nonlinear algebraic system (Equation (140)) representing the dynamic equilibrium equations will be integrated forward in time in order to obtain a solution for each time increment. Let ½0; tmax be the total time interval on which the solution of the algebraic system will be found. This time interval is partitioned as follows: ½0; tmax ¼

N [

½tn ; tnþ1

ð145Þ

n¼1

with Dtn ¼ tnþ1 tn being the nth time increment (or time step). The displacement ﬁeld u(tn) is assumed to be known and noted un. Accordingly, the algebraic system (Equation (140)) at current conﬁguration (time tnþ1 ) is rewritten in the form ½M U¨ nþ1 þfRgnþ1 ¼ 0

ð146Þ

where fRgnþ1 ¼ fFint gnþ1 fFext gnþ1 : To solve this system over each time step, various kinds of explicit and implicit (iterative) schemes can be used. However, in metal forming where the material and geometrical nonlinearities are very strong, two solution procedures are generally used: the dynamic explicit (DE) and the static implicit (SI) methods. Hereafter we restrict our discussion to these two widely used solution procedures. For the other ones, the reader is referred to the general monographs dedicated to the nonlinear FEA cited above. 3.06.4.3.1

ðUÞhnþ1 ¼ ðUÞn þ ðDUÞhnþ1

ð148Þ

The displacement ﬁeld over any solid element (V e) is given by the interpolation deﬁned by Equation (129a), and the corresponding total strain ﬁeld is given by ðee Þhnþ1 ¼ ½Be ðU e Þhnþ1

ð149Þ

The same calculation should be made for the contact elements if there is any. (ii) Given the strain ﬁeld (Equation (149)) at each integration (or Gauss) point of the element under consideration, compute the stress ðsÞhnþ1 and all the internal variables ðaÞhnþ1 using any local integration scheme (see later). (iii) Evaluate the internal force vector ðFint Þhnþ1 and the external force vector ðFext Þhnþ1 and assemble the contribution of all elements using Equations (141b) and (141c). (iv) If Equation (147) is satisﬁed, i.e., Rhnþ1 ¼ 0; then the calculated quantities ðzÞhnþ1 are a good solution. (v) Otherwise modify ðDUÞhnþ1 ; set h-h þ 1 and go to (i). In fact, this procedure is based on the linearization of Equation (147) by assuming the following approximation in which hþ1 h ðDUÞhþ1 nþ1 ¼ ðUÞnþ1 ðUÞnþ1 : @Rhnþ1 hþ1 ðDUÞnþ1 þ?¼0 h @Unþ1 h 1 h hþ1 - ðDUÞnþ1 ¼ Knþ1 Rnþ1

Rhnþ1 þ

ð150Þ

where the global tangent stiffness matrix at the time tnþ1 and iteration (h) deﬁned by

SI resolution strategy

When the inertia effect is neglected, the system (Equation (146)) reduces to Rnþ1 ¼ fFint gnþ1 fFext gnþ1 ¼ 0

in which the internal and external force vectors are given by Equations (143) and (144). The problem to be solved over the time increment Dtn is as follows. Find the displacement increment Dun ; the updated displacement ﬁeld unþ1 ¼ un þ Dun ; the updated stress ﬁeld rnþ1 ¼ rn þ Drn ; and the updated internal variable ﬁelds anþ1 ¼ an þ Dan with aAfr; a; d; c; bg is any one of the internal variables, such that the equilibrium equations (Equation (147)) and the complete set of constitutive equations (see Section 3.06.3.2) hold. The solution to this problem can be obtained by an iterative procedure according to the full Newton method, which proceeds as follows (ðzÞhnþ1 being the value of the variable (z) at the time tnþ1 and the hth iteration). (i) Let ðDUÞhnþ1 be the incremental nodal displacement at the hth iteration giving the total nodal displacement:

ð147Þ

" ½K hnþ1 ¼

@Rhnþ1 h @Unþ1

# ð151Þ

Computational Aspects e h KIV nþ1 ¼

is constructed by assembling the overall elementary tangent stiffness matrices deﬁned for a typical element (V e) by ½K e hnþ1 ¼

@

@

e Fint

V r ½B

¼

þ

e h þ Fext nþ1

@ðU e Þhnþ1

R

e T

fse gJv dV r

h

@

nþ1

@

@ðU e Þhnþ1 h e e T fv Jv dV r V r ½N

nþ1

þ R

þ

@ðU e Þhnþ1 Gr ½N

h fF e gJs dG

e T

nþ1

@ðU e Þhnþ1

ð152Þ

Keeping in mind Remark 8, each of the four derivatives on the right-hand side of Equation (152) gives rise to various terms: h h e h e h ½K e hnþ1 ¼ KIe nþ1 þ KIIe nþ1 þ KIII þ KIV nþ1 nþ1 ð153Þ

where e h KI nþ1 ¼

h @ ½Be T

Z

nþ1

@ðU e Þhnþ1

Vr

þ

Z

½Be T

Vr

þ

Z

@ðU e Þhnþ1

½Be T fse g

Vr

fse gJv dV r

@ ðfse gÞhnþ1

Jv dV r

@ ðJv Þhnþ1 @ðU e Þhnþ1

dV r

ð154Þ

represents the contribution of the elastoplastic behavior to the structural stiffness; h T Bec e nþ1 sc Jc dGrc h e @ðU Þnþ1 Vr h Z e T @ sec nþ1 1 þ B Jc dGrc 2 V r c @ðU e Þhnþ1 Z e T e @ ðJc Þhnþ1 1 þ B sc dGrc 2 Vr c @ðU e Þhnþ1

e h 1 KII nþ1 ¼ 2

Z

@

ð155Þ

h @ fve nþ1 ½N e T Jv dV r @ðU e Þhnþ1 Vr Z @ ðJv Þhnþ1 þ ½N e T fve dV r @ðU e Þhnþ1 Vr

@ðU e Þhnþ1

Js dG

@ ðJs Þhnþ1 ½N e T fF e g dG @ðU e Þhnþ1 Vr

ð157Þ

Remark 9. For coupled thermomechanical problems, i.e., when the thermal effect is taken into account, the temperature should be a nodal degree of freedom (dof). Consequently, the derivatives in Equations (154)–(157) should be performed with respect to both the nodal displacement vector and the nodal temperature vector (see Zienkiewicz et al., 1981; van der Lugt and Hue´tink, 1986, among many others). The external virtual power has contribution from body forces (Equation (156)) and tractions (Equation (157)). The body forces (self-weight or gravity loading) are shown to be deformation independent and therefore its derivative with respect to the displacement ﬁeld is zero. Alternatively, usually if the load (or time) increment is sufﬁciently small, the actual conﬁguration can be supposed to have very small geometrical variation. Consequently and as a ﬁrst approximation, the derivatives of Be, Bec and Jv ; Js ; Jc with respect to the displacement over the time increment of concern can be neglected for simplicity. Accordingly, the corresponding terms in Equations (154)–(157) are zero. This leads to a very helpful simpliﬁcation in the calculation of the tangent stiffness matrix by Equation (153). The error resulting from this hypothesis is shown to be very small as long as the smallness of the geometrical variations over the time step holds and the geometrical difference between the reference and corresponding deformed elements is also small (Hibbitt et al., 1970; Bonet and Wood, 1997). In this work we limit our discussion to the contribution of the internal stress contributions. For a given solid element the contribution (Equation (154)) writes

deﬁnes the contact and friction contributions to the global stiffness; e h KIII nþ1 ¼

Z

@ ðfF e gÞhnþ1

is the external applied surface force (traction) contribution to the global tangent stiffness matrix.

nþ1

c

½N e T

Vr

þ

@ðU e Þhnþ1 R h T @ Gr Bec sec Jc dGrc R

343 Z

Z

KIe

h nþ1

¼

Z

½Be T

Vr

¼

Z

½Be T

@ ðfse gÞhnþ1 @ðU e Þhnþ1

Jv dV r

@ ðfse gÞhnþ1 @ ðfee gÞhnþ1

@ðee Þhnþ1 @ðU e Þhnþ1 Z h e ¼ ½Be T Cnþ1 ½B Jv dV r

Jv dV r

Vr

ð158Þ

Vr

ð156Þ

is the contribution of the external applied body forces; and

h in which the matrix Cnþ1 represents the socalled material Jacobian matrix deﬁning the differentiation of each component of the Cauchy stress tensor with respect to each

344

Computational Damage Mechanics: Application to Metal Forming Simulation

component of the total strain tensor (the superscript (e) referring to the speciﬁc element being omitted): h @ ðfsgÞhnþ1 Cnþ1 ¼ @ðeÞhnþ1

ð159Þ

In order to calculate this matrix as well as the internal force vector (ﬁrst term on the righthand side of Equation (143)), one should ﬁrst compute the updated value of the stress at the time tnþ1 as will be discussed later. Given the stress ﬁeld, two main procedures can be used to calculate the matrix (Equation (159)). Analytical procedures. The objective is to h by obtain a closed-form expression for Cnþ1 differentiating the stress tensor. Two methods are usually invoked. Continuum tangent matrix. This natural method uses the time derivative of the stress given by Equation (66) (without temperature h ¼ effect in this case) and leads to ½Cnþ1 ep ep h ½C nþ1 where C is the fourth-order tangent operator given by Equations (67), (69), (71), (83), (111), (113), or (115) according to the model used. However, it has been shown that the use of this continuum form of the tangent matrix leads to a linear asymptotic rate of convergence of the global iterations in the Newton scheme as deﬁned by Equation (150) (see Nagteggal, 1982; Simo and Taylor, 1985). Alternatively, it is not possible to derive this continuum operator for the time-dependent plasticity or viscoplasticity because of the absence of the consistency condition in these cases (see Remark 4). Consistent tangent matrix. First introduced in Simo and Taylor (1985), and widely used since then, this approach starts by differentiating the discretized expression of the stress tensor over the time interval instead of its continuous expression. Simo and Taylor (1985) showed that this consistent tangent matrix preserves the quadratic rate of convergence of the global Newton scheme. However, its closed-form expression is not easy to obtain for some complex models (see Section 3.06.4.3.2) and strongly depends on the constitutive equations including the type of the objective derivatives, and the time discretization scheme used. Numerical procedures. This requires numerical computation of the differential of each component of the stress tensor with respect to each component of the total strain tensor according to Equation (159). This numerical differentiation is mostly based on the numerical perturbation technique (Kojic and Bathe, 1987). Its main advantage is the independence from the used constitutive equations, objective

derivatives, as well as the time discretization schemes. However, its main drawback is that it needs a large CPU time compared to the analytical closed form procedures. Let us examine now the contribution in Equation (152) of the contact element deﬁned by Equation (155). A similar development as for the solid element leads to h e T @ sec nþ1 Bc Jc dGrc @ðU e Þhnþ1 Vr Z e T c h e ¼ Bc Cnþ1 Bc Jc dGrc

e h KII nþ1 ¼

Z

ð160Þ

Vr

in which the differential of the contact matrix Bec has been neglected because it is not simple to compute because of its dependency on the rotation matrix (see Equation (133)). The contact Jacobian matrix in Equation (160) is deﬁned by c

@ ðfsc gÞhnþ1 h Cnþ1 ¼ @ðec Þhnþ1

ð161Þ

This matrix is often calculated numerically by an appropriate perturbation technique or analytically by differentiating the expression relating the stress to the total strain at nodes lying on the contact boundaries. Of course, the use of these different methods leads to different numerical forms of the tangent stiffness matrix, but this should have no effect on the accuracy of the solution of the IBVP. However, this has a great inﬂuence on the convergence rate of the global Newton solution strategy, particularly for highly nonlinear IBVP problems including material and geometrical nonlinearities such as large plasticity, hardening, damage, and evolving contact with friction. Remark 10. Note that the tangent matrix as deﬁned above is generally nonsymmetric because the material Jacobian matrices themselves are nonsymmetric. Hence, a nonsymmetric equilibrium solver is generally required. Alternatively, the Newton-type implicit iterative resolution strategies are unconditionally stable and allow us to use relatively large time or load increments. Finally, it is important to mention that many other implicit iterative global resolution schemes can be derived from the full Newton method presented above. They are called the modiﬁed Newton methods with the quasiNewton method as a most widely used one, which reduces signiﬁcantly the number of iterations. These quasi-Newton methods differ mainly by the way with which the tangent

Computational Aspects matrix (Equation (159)) is approximated (for more details, see the general nonlinear FEA monographs cited above). Particularly, the ‘‘arc-length’’ methods, originally introduced by Riks (1978), are suitable to solve physical problems involving limit points (see Crisﬁeld, 1991 for a complete review on these continuation techniques). It is worth noting that the fully coupled elastoplastic-damage models developed in Section 3.06.3 exhibit stress–strain (or force–displacement) responses with a limit point because of the damage-induced softening behavior. 3.06.4.3.2

Let us go back to the dynamic discrete system given by Equation (140) in order to discuss its resolution using explicit integration schemes. These aim to obtain values for dynamic quantities at the end of a speciﬁc time interval tnþ1 ¼ tn þ Dtn based exclusively on values available at time tn. The most commonly used explicit operator for stress analysis is the central difference formula. Because it is conditionally stable, and according to the high nonlinearity of the system (Equation (140)), this explicit scheme should be used with an automatically varying time increment Dtn. Written at the nth time increment, the discretized dynamic system (Equation (14)) is rewritten as follows: ½M nþ1 U¨ nþ1 þfFint gnþ1 fFext gnþ1 ¼ 0

(3) Update displacements; ðDtn Þ2 U¨ n fU gnþ1 ¼ fU gn þDtn U’ n þ 2

nþ1

¼ ½ML 1 fFint gnþ1 fFext gnþ1

ð164Þ

(8) Compute the velocity ﬁeld: Dtn U’ nþ1 ¼ U’ n þ U¨ n þ U¨ nþ1 2

ð165Þ

(9) If totmax go to (2), if not STOP. This explicit scheme requires no iteration procedure and no tangent stiffness matrix, but is conditionally stable as mentioned above. Accordingly, the automatic time step control (step (2)) is essential to control the stability limits and accuracy of the obtained solution (Hughes, 1983). The basic stability limit is shown to be inversely proportional to the highest eigenvalue in the global system: Dtr

ð162Þ

in which the mass matrix (Equation (135) (generally symmetric, positive deﬁnite, but nondiagonal) is called consistent mass matrix in the sense that it is calculated from the same shape functions N that are used for the interpolation of the displacement ﬁeld. For several dynamic solution problems, it will be computationally advantageous to replace this consistent matrix by a diagonal one called the lumped mass matrix [ML]. This consists of placing concentrated point masses mi at nodes i in the directions of the assumed displacement dof’s. The concentrated masses refer to translational and rotational inertia of the element. For more details of mathematical methods used to lumping the mass matrix and its consequences for the FE solution accuracy, the reader is referred to Hughes (1987), among many others. By using the lumped mass matrix, the DE solution strategy proceeds as follows: (1) Lump the mass matrix and set initial conditions; (2) Estimate the critical time step Dtn;

ð163Þ

(4) Compute the stress and internal variables at the time tnþ1 by integrating the full set of coupled constitutive equations of concern deﬁned in Section 3.06.3 (see later); (5) Compute the internal forces for all elements (Equations (141b) and (143)); (6) Compute the external nodal forces (Equations (141c) and (144)); (7) Solve for accelerations: U¨

Dynamic explicit resolution strategy

345

2 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ2 1þW W

omax

ð166Þ

where omax is the highest frequency in the system and W is the fraction of critical damping in the highest mode. In practice, there is no need to solve the eigenvalue problem for the whole system in order to determine the highest eigenvalue. Instead, it can be estimated by determining the maximum element dilatational mode of the mesh; further the time step is estimated by

e L Dtrmin _ cd

ð167Þ

where _ is a mesh-dependent stability factor, Le is the smallest distance between adjacent nodes of any element (e) of the mesh, also called the characteristic element size, and cd is the current dilatational wave speed of the material. 3.06.4.4

Local Integration of the Fully Coupled Constitutive Equations

Each of the SI or dynamic explicit (DE) solution procedures outlined above contains a crucial step, namely, step (2) for SI and step (4)

346

Computational Damage Mechanics: Application to Metal Forming Simulation

for DE. This concerns the computation of the stress snþ1 ¼ sn þ Dsn at the end of the time interval tnþ1 knowing its value sn at the beginning of the same time interval tn. These steps depend crucially on the used constitutive equations developed in Section 3.06.3 and aim to integrate them numerically at each Gauss point lying within a solid element (e). In the wide literature concerned with the (local) numerical incremental integration of elastoplastic constitutive equations, one can ﬁnd various explicit and implicit iterative schemes (see the general monographs cited above). In fact, these constitutive equations deﬁne a set of ODEs which can be numerically integrated by using the generalized midpoint rule. In what follows, we shall restrict our attention to the most widely used implicit incremental integration scheme called the ‘‘return mapping algorithm’’ (RMA) (see the monographs by Crisﬁeld, 1991, 1997; Simo and Hughes, 1998 and references therein for more details). In this scheme, the incremental time integration of the time-independent constitutive equations (as those developed in Section 3.06.3) over the time step ðtn ; tnþ1 ¼ tn þ Dtn Þ is regarded as a strain-driven process in which the total strain is the basic variable. From the FEA point of view, this integration procedure is performed at each time step and for each Gauss point inside each element. In the present case of large plastic deformations, the constitutive equations being formulated in the rotated frame (see Section 3.06.2.1 and the beginning of Section 3.06.3.1), the integration procedure can be stated as follows: (1) Let x be a given integration or Gauss point located inside the element of interest; (2) Assume that values of all the state variables (including the stress ﬁeld rn and r all ‘‘i’’ internal variables symbolically denoted ðiÞ by an ) at the point x are known at the beginning of the time step tn; (3) Give the time increment Dtn and the associated total strain increment Den at the point x; (4) Rotate all the tensorial variables by the rotation Q solution of Equation (4) to change from the current deformed conﬁguration (global equilibrium coordinates) to the deformed rotated one; (5) Apply the RMA to compute the stress and internal variable increments in order to update them to time tnþ1 in a manner consistent with the constitutive model under consideration; and (6) Rotate back to the global coordinates. In what follows only the step (5) dealing with the incremental fully implicit integration will be applied to fully coupled constitutive equations

developed in Section 3.06.3. For the sake of simplicity and brevity, only the uniﬁed (single yield surface) isothermal and fully isotropic elastoplastic-damage model discussed in Section 3.06.3.4.2 will be considered (Zhu and Cescotto, 1995; Hammi, 2000 for the incremental integration of the two yield surface model). 3.06.4.4.1

Application to the isotropic single yield surface model

In this case the fully coupled constitutive equations are given in Section 3.06.3.4.2. Note that this integration procedure has been applied to similar coupled constitutive equations under small-strain hypothesis by Simo and Ju (1987), Benallal et al. (1988), and Doghri (1995). Keeping in mind the isotropic elasticity with the elasticity operator K given by Equation (85), the elastic predictor stage is deﬁned by the so-called ‘‘trial’’ state in which the trial stress, using the classical compact tensorial notation, is deﬁned by p rtrial nþ1 ¼ ð1 Dn ÞK : enþ1 en ¼ ð1 Dn Þ Ke enþ1 : 1 2Ge e&dn

ð168Þ

in which enþ1 ¼ en þ Den is completely known, as are all other variables containing the subscript (n) relative to the beginning of the time increment. The deviatoric strain tensor e& dn ¼ ednþ1 epn is also entirely known over the time step. When replaced in the single yield surface (see Equation (39)), one gets the ‘‘trial’’ yield function: ftrial ¼

jjrtrial nþ1 X n jjs Rn pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sp 1 Dn

ð169Þ

If ftrialo0, then the trial solution holds (elastic step) and a new time step is outworked p with the updated variables: rnþ1 ¼ rtrial nþ1 ; enþ1 ¼ p p p epn ; X nþ1 ¼ X pn ; Rnþ1 ¼ Rpn and Dnþ1 ¼ Dpn : If ftrial is equal to or greater than zero, the trial solution above cannot be a solution to the incremental problem. This means that over the actual time step the problem is incrementally plastic characterized by f ðrnþ1 ; X nþ1 ; Rnþ1 ; Dnþ1 Þ ¼ 0 and Dlnþ1 40: The objective then is to determine the unknowns rnþ1 ; X nþ1 ; Rnþ1 ; Dnþ1 ; and Dlnþ1 at the end of the time step by using the RMA. To accomplish this task we express the state variables at the end of the time increment as follows (see Hammi, 2000; Saanouni et al., 2001a, 2001b, 2001c): rdnþ1 ¼ ð1 Dnþ1 Þr& dn pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2Ge 1 Dnþ1 Dlnþ1 nnþ1

ð170Þ

Computational Aspects X nþ1 ¼ ð1 Dnþ1 ÞX n expðaDlnþ1 Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2C þ ð1 expðaDlnþ1 ÞÞ 1 Dnþ1 nnþ1 3a ð171Þ Rnþ1 ¼ ð1 Dnþ1 ÞQRn expðbDlnþ1 Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Q þ ½1 expðbDlnþ1 Þ 1 Dnþ1 b Dnþ1 ¼ Dn þ

Ynþ1 s Dlnþ1 S ð1 Dnþ1 Þb

347

in which the deviatoric tensor Z at the time tn has been introduced: Z n ¼ r& dn 23Can expðaDlnþ1 Þ

ð177Þ

Multiplying by rdnþ1 X nþ1 taken from Equation (175) leads to ð172Þ

h pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ i rdnþ1 X nþ1 ¼ 23 Rnþ1 þ 1 Dnþ1 sp nnþ1

ð178Þ

This allows us to write ð173Þ Z nþ1 ¼ jjZ nþ1 jjs nnþ1

where

ð179Þ

in which the norm of the tensor Z is given by r& dn

¼ 2Ge enþ1

d epn ¼

2Ge e&dn

ð174Þ

To be complete, one should add the yield function (Equation (39)) written at the end of the time increment: fnþ1 ¼

jjrnþ1 X nþ1 jjs Rnþ1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sp 1 Dnþ1

ð175Þ

In obtaining Equations (171) and (172), a closed-form asymptotic integration applied to both kinematic- and isotropic-hardening evolution equations has been used (see details in Hammi, 2000 or Saanouni et al. 2001a, 2001b, 2001c). This semi-analytical integration method has been initially proposed by Walker (1987) for viscoplasticity and used by many others since then (Chaboche and Cailletaud, 1996; Nesnas and Saanouni, 2002). In 3D case, Equations (170)–(175) lead to a nonlinear system of 15 scalar equations to be linearized and solved iteratively according to the Newton solution scheme. However, an important simpliﬁcation of the problem can be made according to an idea originated by Simo and Taylor (1985) and widely used since then for advanced elastoplastic or elastoviscoplastic constitutive equations (Hartmann and Haupt, 1993; Doghri, 1993, 1995; Chaboche and Cailletaud, 1996; Hammi, 2000; Saanouni and Hammi, 2000). This simpliﬁcation aims to exploit some special properties of the used constitutive equations to reduce the system Equations (169)–(173) from 15 to only two scalar equations in this fully isotropic case (see Hammi, 2000 for the general case of fully anisotropic two-surface formulation). Taking the difference between Equations (170) and (171), one gets pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ rdnþ1 X nþ1 ¼ ð1 Dnþ1 ÞZ n 1 Dnþ1 2C 2Ge Dlnþ1 þ 3a ð1 expðaDlnþ1 ÞÞ nnþ1

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 Q jjZnþ1 jjs ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 ð1 bRn 1 Dnþ1 Þ 1 Dnþ1 b 1 þ sp þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3Ge Dlnþ1 1 Dnþ1 C þ ð1 expðaDlnþ1 ÞÞ ð180Þ a

which is a function of the only two unknowns, namely, Dlnþ1 and Dnþ1 : Hence, with the help of Equation (179), the variable n can be replaced by the variable Z. Accordingly, we have to solve a small system of two scalar equations— g1 ðDlnþ1 ; Dnþ1 Þ and g2 ðDlnþ1 ; Dnþ1 Þ—deﬁned from Equations (173) and (180) as follows: g1 ðDlnþ1 ; Dnþ1 Þ ¼ jjZ nþ1 jjn

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 Q 1 ð1 bRn 1 Dnþ1 Þ þ sp pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 Dnþ1 b 1 C pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3Ge Dlnþ1 þ ð1 expðaDlnþ1 ÞÞ a 1 Dnþ1 ¼0 g2 ðDlnþ1 ; Dnþ1 Þ Ynþ1 s Dlnþ1 ¼0 ð181Þ ¼ Dnþ1 Dn S ð1 Dnþ1 Þf

where the damage energy release rate Ynþ1 is obtained from Equation (92a) (for isothermal case) and rewritten in the following discretized form: Ynþ1 ¼ 2Ge þ

e&dn

:

e&dn

3ðDlnþ1 Þ2 Dlnþ1 nnþ1 : e&dn þ 2ð1 Dnþ1 Þ

C an expðaDlnþ1 Þ 3

1 þ ð1 expðaDlnþ1 ÞÞnnþ1 a Q rn expðbDlnþ1 Þ þ 2

2

1 þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃð1 expðbDlnþ1 ÞÞ b 1 Dnþ1 ð176Þ

!

2 ð182Þ

in which the normal n should be replaced by Z, which is a function of Dlnþ1 and Dnþ1 only. The

348

Computational Damage Mechanics: Application to Metal Forming Simulation

system given by Equation (181) should be solved by using the iterative Newton–Raphson procedure in order to determine the unknowns Dlnþ1 and Dnþ1 (see Hammi, 2000 for detailed calculations). Remark 11. For plane-stress hypothesis the total strain component e33 should be deﬁned by an additional constraint, namely, g3 ðDl; D; e33 Þ ¼ s33 ¼ 0: This leads to a new scalar equation with the new unknown e33 to be determined together with Dl and D by solving the three equations—g1, g2, and g3 (see Hammi, 2000). When Dlnþ1 and Dnþ1 are computed as discussed above, Equations (170)–(172) are used to determine explicitly the value of rnþ1 ; X nþ1 ; Rnþ1 at the end of the time step. Hence, the integration of all constitutive equations is terminated. Recall that this local integration procedure is required with both the SI global solution scheme (step (2) in Section 3.06.4.3.2) and the DE solution strategy as discussed in step (4) of Section 3.06.4.3.2. However, with the SI scheme the material Jacobian matrix deﬁned by Equation (159) is also required. This will be discussed in the next section. 3.06.4.4.2

Consistent elastoplastic tangent operator

To complete the algorithmic procedure for the SI solution strategy discussed in Section 3.06.4.3.2, there only remains to compute the closed-form expression for the consistent tangent modulus deﬁned by Equation (159). This can be done by differentiating the expression of

dDnþ1 ¼ denþ1

* nþ1 2Ge ð1 Dnþ1 Þ Cnþ1 ¼ K

dlnþ1 dnnþ1 þ Dlnþ1 nnþ1 # denþ1 denþ1 Ke ðenþ1 : 1Þ1 þ r& dn 2Ge Dlnþ1 nnþ1 dDnþ1 ð184Þ # denþ1

where n has to be replaced by Z (Equation (179)) which is function of Dlnþ1 and Dnþ1 ; which in turn are the linearized system deﬁned by Equation (181). The partial derivatives entering Equation (184) are obtained by taking the differential of both g1 and g2 deﬁned in Equation (181): dg1 dg1 dlnþ1 dg1 dDnþ1 ¼ þ denþ1 dlnþ1 denþ1 dDnþ1 denþ1 dg1 dZ nþ1 þ ¼0 dZ nþ1 denþ1 dg2 dg2 dlnþ1 dg2 dDnþ1 ¼ þ denþ1 dlnþ1 denþ1 dDnþ1 denþ1 dg2 dYnþ1 þ ¼0 dYnþ1 denþ1

f

Y s1 S

! ðdY =dlnþ1 Þ

ððdg1 =dDnþ1 Þðdg2 =dlnþ1 Þ ðdg1 =dlnþ1 Þðdg2 =dDnþ1 ÞÞ

ð186Þ

dlnþ1 ððdg1 =dDnþ1 ÞðdDnþ1 =denþ1 Þ þ 3Ge nnþ1 Þ ¼ dg1 =dlnþ1 denþ1 ð187Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ drnþ1 d ð1 Dnþ1 Þ Ke ðe : 1Þ þ r& dn 2Ge 1 Dnþ1 Dlnþ1 nnþ1 ¼ denþ1 denþ1

in which n should be replaced by Z depending on Dlnþ1 and Dnþ1 : Because Dlnþ1 and Dnþ1 are dependent on the total deformation, the

ð185Þ

This leads, after some manipulations, to

3Ge ðdg2 =dlnþ1 Þnnþ1 þ sDlnþ1 =Sð1 DÞ

the stress at the time tnþ1 ; the deviatoric part of which is given by Equation (170) (for clarity the superindex h for the iteration number will be omitted from Equation (159)):

Cnþ1 ¼

differential in Equation (183) requires the differentiation of Dlnþ1 and Dnþ1 (solution of the system Equation (181)) with respect to the total strain. These latter can be obtained from the differentiation of g1 ðDlnþ1 ; Dnþ1 Þ and g2 ðDlnþ1 ; Dnþ1 Þ deﬁned in the system Equation (181) with respect to enþ1 : Despite the apparent complexity of the fully coupled constitutive equations, they can be manipulated to give (see Hammi, 2000)

ð183Þ

However, the derivatives of n in Equation (184) and of Y in Equation (186) with respect to e are given by

Computational Aspects

dnnþ1 denþ1

1 3 d 1 nnþ1 #nnþ1 : ¼ jjZ nþ1 jjs 2 dlnþ1 2Ge 1d þ a expðaDlnþ1 ÞX n # denþ1

ð188Þ

dYnþ1 2Ge Dlnþ1 3 d ¼ 1 nnþ1 #nnþ1 : denþ1 jjZ nþ1 jjs 2 Xn dlnþ1 e # K : enþ1 ð1 þ aDlnþ1 Þ denþ1

ð189Þ

This completes the calculation of the closedform expression of the consistent tangent operator deﬁned by Equation (184) and required to compute the tangent stiffness matrix needed in the implicit Newton–Raphson resolution strategy. Note that in the fully coupled elastoplastic-damage constitutive equations, the consistent tangent operator (Equation (184)) is quite different from its continuous counterpart expression deﬁned by Equation (122) for the fully isotropic and isothermal case.

3.06.4.4.3

On the incremental objectivity

As speciﬁed above, the frame indifference of the fully coupled constitutive models, outlined in Section 3.06.3, is ensured, thanks to the local rotated frame formulation as discussed in Section 3.06.2. The constitutive equations are transformed to a locally Cartesian rotating frame (by the orthogonal tensor Q deriving from Equation (4)) which ensures that any rotated tensor (Equation (6)) remains unaltered under superimposed rigid body motions. Hence, the time-stepping algorithms discussed above (Section 3.06.4.3) should be performed in this locally rotated conﬁguration and the discrete equations are mapped back to the Eulerian unrotated conﬁguration. This is called the incremental objectivity ﬁrst formalized by Hughes and Winget, 1980 and widely used in computational plasticity since then. In this rotated frame it has been shown (see Hughes and Winget, 1980; Simo and Vu-Quoc, 1986) that the fully implicit midpoint-rule-based algorithm (backward) discussed above is incrementally objective. However, it remains to integrate the evolution equation governing the orthogonal rotation Q (Equation (4)) over the time step ðtn ; tnþ1 Þ: Details concerning this point are given in Simo and Hughes (1998) (chapter 8) where it has been shown that for large elastic and plastic deformations, the generalized radial return algorithm (y method) is such that

349

(i) it is consistent with the ODE (Equation (4)) for all yA½0; 1; and second-order accurate for y ¼ 0:5; (ii) the exponential map of the orthogonal rotation tensor Q possesses a singularity which can be easily avoided by an appropriate parametrization; (iii) it generates an incrementally orthogonal transformation regardless of the y values; (iv) it places no restrictions on the value of the admissible rotation increments; and (v) for ﬁnite elastic deformation based on the hypoelastic behavior, the closed-form linearization of the algorithm is not easy, and the consistent tangent moduli is not available. However, this becomes possible if the hyperelastic behavior is used. Accordingly, the fully implicit (y ¼ 1) radial return mapping algorithm should be performed in the locally rotated conﬁguration at time tnþ1 ; exactly in the same manner as presented in Section 3.06.4.3 above. Keeping in mind that the elastic deformation component is assumed to be inﬁnitesimal (see Equation (8)), the time derivatives are simply related to the objective rotational derivatives by Equation (7). 3.06.4.5

On the Numerical Aspects of Contact-friction Problems

First we note that the easiest and most commonly used approach to model the contact-friction in the FEA codes is the use of speciﬁc contact elements which are 1D elements for 2D problems and 2D (triangular or quadrilateral) for 3D contacting bodies. In this section, we limit ourselves to give some indications of the numerical aspects of the contact problems. For more details, the reader is referred to the specialized literature devoted to the algorithmic aspects of contact-friction problems (Hughes et al., 1976; Hallquist et al., 1985; Chen and Tsai 1986; Chaudhary and Bathe, 1986; Belytschko and Lin, 1987; JMTA, 1988; Simo and Laursen, 1993; chapter 10 in Belytschko et al., 2000, among many others). In Section 3.06.3.5, the main contact-friction constitutive equations have been reviewed. The associated numerical aspects have been discussed in Section 3.06.4.1 where it has been shown that the contact-friction contributes to the virtual power principle a part of the internal force vector (Equation (136)). Even if this contribution is limited to nodes lying in the contact elements at which the equilibrium equations, the plastic yield equations, as well as the contact-friction equations should be identically satisﬁed, it introduces additional material (friction law) and kinematical (varying contact conditions) nonlinearities.

350

Computational Damage Mechanics: Application to Metal Forming Simulation

Particularly, the contact-friction makes an important contribution to the tangent stiffness matrix (see Equation (152)) in the iterative implicit resolution scheme and participates, via the internal forces vector, in the solution of the IBVP by both SI or DE solution strategies. Accordingly, the contact with friction adds two new numerical aspects to the discretized IBVP: (i) the treatment of the contact conditions (Equation (124)) and (ii) the contact search methodology. Concerning the contactfriction solution methods, they aim to solve the Signorini conditions (Equation (127)). Three main solution methods are generally used for contact-friction problems involving large geometry changes: (i) the Lagrange multiplier method, (ii) the penalty method, and (iii) the augmented Lagrangian method. Each has its own advantages and drawbacks. For the contact search algorithms, the main problem is to determine at a given time the nodes concerned with the contact keeping in mind that during the forming processes the contact area changes continuously. It consists in ﬁnding the gap d (distance) separating two points (nodes or Gauss points) lying in the contact area. This problem is very complex and CPU time consuming, especially for 3D computations. In general, the forming tools deﬁne ‘‘master’’ surfaces and the deforming part deﬁnes the ‘‘slave’’ surfaces. Both can be considered as elastically or inelastically deformable bodies or only the part is assumed to be deformable while the master surfaces are assumed to be rigid bodies. The boundaries of the slave surface are always discretized by FE (facets), while the master surfaces are either discretized by FE or represented by analytical surfaces as elementary geometric forms (planes, cylinders, spheres) or by using Bezier patches, B-splines, etc. The main FEA computer codes concerned with the contact-friction problems offer the possibility to solve: (i) the small-sliding problems and/or (ii) the ﬁnitesliding problems, with an arbitrary rotation (ﬁnite or small) of the contacting bodies being permitted. Detailed discussions of these aspects can be found in user’s manuals of the used FEA codes: Abaqus/std and Abaqus/Explicit (see ‘‘Abaqus theory manual,’’ chapter 5, Abaqus Users Manual, vol. 1). 3.06.4.6

On Some Aspects of Remeshing and Adaptive Meshing

Numerical errors are intrinsic to FEA of any physical and mechanical problems. In fact, the time and space discretizations of the continuous PDEs and ODEs governing the physical events lead inevitably to numerical errors. How

can these errors be deﬁned, measured, controlled, and minimized? These questions have confronted specialists in numerical simulations since the beginning of numerical method applications. These aspects have been extensively studied in the literature since the early 1980s, as indicated by a considerable number of published papers, and a growing number of conferences and special sessions devoted to these problems. Since the pioneering work of Babuska and Rheinboldt (1978), a veritable new ﬁeld in computational mathematics and mechanics has emerged offering new theories and methods of a priori and a posteriori error estimation. It is out of the scope of this chapter to give an exhaustive survey of this ﬁeld; indeed, we limit ourselves to a very brief description of the four related aspects indicated above. The reader is invited to make his own exploration of the literature by consulting the following general monographs dedicated to this ﬁeld and giving many related references (Babuska and Rheinboldt, 1982; Zienkiewicz and Zhu, 1987; Gallimard et al., 1996; Ainsworth and Oden, 1997; Ladeve`ze and Oden, 1998; Ladeve`ze and Pelle, 2001, and so on). Alternatively, the algorithmic aspects related to mesh generation (2D and 3D surface and volume) can be found in Baker (1989), George (1991), and George and Borouchaki (1997). Compared with geometrically linear structural mechanics problems, metal forming problems possess additional difﬁculties related to the very large changes in the geometry of the computational domain. In fact, the deformed parts undergo geometrical variation (large plastic deformations and rotations) and are characterized by inhomogeneous spatial distribution of thermomechanical ﬁelds with evolving localized zones (stress, plastic strain, damage, temperature, etc.). Accordingly, several remeshings have to be performed during the simulation in order to preserve the reliability of the obtained results by minimizing errors generated by either the geometrical transformations or the heterogeneous thermomechanical ﬁelds. As of early 2000s, the most published works in this ﬁeld are related to linear and nonlinear problems solved over a ﬁxed domain, i.e., a domain submitted to small geometric variations with respect to its initial conﬁguration. Concerning 2D problems, the remeshing, adaptive meshing, and related topics are solved from both theoretical and practical aspects. An increasing number of general purpose or specialized FEA codes offer very helpful remeshing and adaptive meshing facilities with the possibility of controlling the entire numerical simulation process through adaptive

Computational Aspects algogithms (see, e.g., George and Borouchaki, 1997). However, for the general 3D (surfaces and volumes) remeshing and adaptive meshing of the geometrically complex structures, fewer works have been published. Particularly, many efforts continue to be made in order to solve some open problems related to the 3D adaptive meshing in metal forming (Zienkiewicz et al., 1988; Cheng, 1988; Yang et al., 1989; Habraken and Cescotto, 1990; Coupez and Chenot, 1992; Dyduch et al., 1992; Fourment and Chenot, 1994; Dyduch et al., 1995; Espinosa et al., 1998). Note that various mesh reﬁnement methods can be used: R-adaptivity. After a large distortion of elements, only node positions are optimized in order to keep only geometrically acceptable elements according to some intrinsic geometrical criteria for each kind of elements. Consequently the number of elements (and nodes) as well as the associated connectivity are kept unchanged. H-adaptivity. The element size is modiﬁed (reﬁned or coarsened) according to some appropriate geometrical or physical criteria. This leads to a varying number of elements (and nodes) with a new connectivity. P-adaptivity. The number of elements is kept constant and the degree of the shape functions is modiﬁed giving a variation in the node number. For mesh reﬁnement the shape function degree (SFD) is augmented leading to an increase in the number of nodes. Similarly, for the mesh coarsening the SFD is decreased, thus decreasing the element node number. HP-adaptivity. It consists in combining both the H-type and the P-type methods in order to increase the mesh adaptivity efﬁciency. Whatever the remeshing and the adaptivity methods, three related aspects have to be considered: the remeshing decision, the new mesh generation and ﬁeld transfer and restart of the analysis. 3.06.4.6.1

Remeshing decision

To decide when the remeshing is required during the analysis, some appropriate criteria are needed and should be automatically executed during the analysis. These are generally based on a priori and a posteriori error estimators (see Ainsworth and Oden, 1997; Ladeve`ze and Oden, 1998) where an exhaustive summary of error estimators is given). The main goal of any error estimator is the evaluation of the absolute global error as an addition of the estimated local error for each element. These error estimators (EEs) can be formally classiﬁed into three families:

351

(i) constitutive relation-based EE, (ii) equilibrium residual-based EE, and (iii) gradientbased EE. The latter seems to be the simplest and widely used one at the engineering level. It aims to obtain a continuous (average) approximation to the solution gradient by postprocessing the gradient of the FE approximation. For metal forming these criteria can be classiﬁed into three groups: (i) Estimation of the element distortion due to the large transformation of the body. Distortion criteria are based on the large variation of the geometry (skewness, angles, elongation, etc.) of FEs with respect to their reference state (George and Borouchaki, 1997). (ii) Estimation of the element size needed to avoid inter-penetration between the deformed part and tools. This geometric estimator is based on the curvature of the tool angles inside the contact zones (Borouchaki et al., 2002). (iii) Adaptation of the element size to the local or global variation of some physical ﬁelds as temperature, displacement gradient, stress, plastic strain, etc. Some works have been devoted to the adaptive mesh reﬁnement for nonlinear problems undergoing local and nonlocal damage (Espinosa et al., 1998; Svedberg and Runesson, 2000; Rodriguez-Ferran and Huerta, 2000; Ridriguez-Ferran et al., 2001) and for metal forming problems (Borouchaki et al., 2002). 3.06.4.6.2

New mesh generation

When the remeshing is required during the FEA according to the fulﬁllment of one (or more) criterion discussed above, the new mesh should be generated. This remeshing operation uses the same mesh generator as the initial mesh of the undeformed conﬁguration. When the domain boundaries are ﬁxed (i.e., inﬁnitesimal transformation hypothesis), the mesh generation in 2D (curves, surfaces) and 3D (surfaces and volumes) is now well founded (see the book by George and Borouchaki, 1997). However, when the domain boundaries undergo large geometrical variations (as in metal forming), additional difﬁculties have to be solved. In 2D metal forming, the major difﬁculties have been satisfactorily solved. However for 3D large deformation problems further developments continue to be made in order to obtain efﬁcient fully automatic mesh generators for complicated 3D shapes with adaptive nonlinear analyses (Coupez and Chenot, 1992). 3.06.4.6.3

Field transfer and FEA restart

After generating the new mesh, the thermomechaical ﬁelds should be transferred from the

352

Computational Damage Mechanics: Application to Metal Forming Simulation

old mesh to the new one, such that the computation may continue without excessive error. To achieve this successfully, two main methods are usually used, namely, the standard interpolation and the projection techniques. However, for problems with high physical (ﬁnite elastoplasticity, damage, friction, etc.) and geometrical (evolving contact, large transformation, etc.) nonlinearities, some improvements should be made. In fact, on the old mesh the concerned ﬁelds satisfy (in the FE sense) both the global equilibrium equations (displacement and temperature ﬁelds at nodes) as well as the constitutive equations (stress and internal variable ﬁelds at integration points). Accordingly, the transfer operation when large deformations are accounted for, should address the following issues (Cheng, 1988; Dyduch et al., 1992; Espinosa et al., 1998; Villon et al., 2000): (i) the satisfaction of the static or dynamic equilibrium equations (static or dynamic admissibilty); (ii) consistency with constitutive equations (plastic admissibility); (iii) compatibility with the displacement or velocity ﬁelds on the new mesh (kinematic admissibility); (iv) compatibility with evolving boundary conditions (unilateral contact conditions) and with the rupture condition (critical value of damage); and (v) minimization of numerical ‘‘diffusion’’ of state variables. 3.06.4.7

Special Treatment of Fully Damaged Elements

According to the constitutive models developed in Section 3.06.3 the ductile damage is ‘‘strongly’’ coupled to the other thermomechanical phenomena as thermoelasticity, plasticity, and hardening. Accordingly, the damage growth induces a decrease in the stress-like variables (not only Cauchy stress but also internal stresses except Y) generated by a decrease in physical properties of the material as elastic and hardening moduli (see Equations (33)–(36)). This means that the fully damaged Gauss points are stress free and no longer contribute to the elementary tangent stiffness matrix (see Equation (159)) or to the internal forces vector (Equation (143)). Hence, this point should be excluded from the integration domain (points) of the element to save CPU time. Straightforwardly speaking, this happens only if the damage is isotropic or if all the components of the anisotropic damage tensor reach their critical values, giving zero rigidity in all the space directions. When all Gauss points

within a given element are fully damaged, the corresponding stiffness matrix (or some of its components if the damage is anisotropic) is zero. Consequently, this element (in the case of isotropic damage) has no more contribution to the global tangent stiffness matrix and should be excluded from the assemblage process. This fully damaged element can be ‘‘killed’’ from the structure as shown in Mariage and Saanouni (2001). Note that in the present local formulation, the IBVP solution is mesh dependent. Consequently, damage localizes always inside one row of elements without any possibility to cover many elements in the direction orthogonal to the damage propagation direction (see Saanouni et al., 1986, 1989; Hammi, 2000, among many others). Accordingly, there is no problem to continue the FEA when fully damaged zones propagate as long as all nodes are still connected to a partly damaged element. However, elements lying on boundaries of the deformed domain need special attention especially when they are concerned with the contact between different bodies (tools and part). Two problems arise in that case: (i) boundary nodes exclusively connected to fully damaged elements should be dropped from the stiffness matrix in order to avoid any singular behavior; and (ii) the contact algorithm is affected when a boundary node is dropped because the contact element connection should be changed. The best way to treat the fully damaged elements consists in remeshing the domain after dropping the fully damaged elements and smoothing the newly created boundaries of the deformed part. To our knowledge this very difﬁcult remeshing problem with new boundaries creation (holes, cracks) has not been yet addressed in the literature. This aspect is under progress in our team at the time of this writing. 3.06.5

NUMERICAL PREDICTION OF THE DAMAGE: SOME APPLICATIONS

The aim of this section is to give some applications of the proposed theory to the numerical prediction of the ductile damage in metal forming processes. Only the isotropic ductile damage coupled to rate-independent elastoplasticity under constant reference temperature (thermal coupling being neglected) will be considered. All the examples presented herein have been solved with Abaqus/Std or Abaqus/Explicit, thanks to the use of user routines Umat and Vumat developed in our laboratory (Hammi, 2000; Saanouni and Franqueville, 1999; Saanouni et al., 2000,

Numerical Prediction of the Damage: Some Applications 2001a, 2001b, 2001c). Our main goal is to illustrate the capability of the proposed methodology to predict when and where ductile damage zones may form inside the deformed part during the process. Moreover, one can show the possibility to act on the process parameters (tools velocity, friction, material ductility, etc.) in order to avoid damage as in forging, deep drawing, hydroforming, etc., or to enhance damage in order to simulate some sheet metal cutting operations. This will be termed as an ‘‘optimization’’ of the process plan operation. Note that, when implementing the damage effect in Abaqus/Std (Umat), Abaqus/Explicit (Vumat), a ﬂag parameter (icoupl) has been introduced allowing us to perform a fully coupled analysis when icoupl ¼ 1 and an uncoupled analysis when icoupl ¼ 0. This allows us to compare the coupling effect using exactly the same analysis computational tools. 3.06.5.1

On the Identiﬁcation of the Material Parameters

In engineering mechanics the use of any kind of constitutive equations in order to simulate numerically the behavior of a given structure by FEA, needs the knowledge of the material parameters entering these constitutive equations. This difﬁcult and as yet not satisfactorily solved problem requires ‘‘automatic’’ computation of the material parameters of concern, by comparison with available experimental database. From a theoretical point of view, this deﬁnes a mathematical optimization problem using appropriate inverse methods termed here as an identiﬁcation procedure. The importance of this identiﬁcation procedure is proportional to the increasing ‘‘sophistication’’ of the socalled ‘‘advanced’’ constitutive equations describing many coupled physical phenomena. Many papers, special sessions, and specialized conferences have been devoted to this very important problem in structural mechanics (see Bock, 1983; Chicone and Gerlach, 1987; Bui, 1992; Bui and Tanaka, 1994; Mahnken and Stein, 1994, among many others). In this work, devoted to the identiﬁcation of the isotropic, isothermal elastoplastic constitutive equations accounting for both mixed (isotropic and kinematic) hardening and ductile damage, the following simple procedure is used for a given material: (i) the used database contains only uniaxial tension tests conducted on a given specimen until the ﬁnal fracture; (ii) the friction constitutive equation is limited to Coulomb model in which the friction parameter Z (see Equation (125)) is taken to be

353

constant as follows: (a) Zo0:1 for lubricated contact, (b) Z ¼ 0:3 for metal–metal natural contact, and (c) Z40:45 for sticking contact; (iii) the elastoplastic parameters (n, E, a, C, b, Q, sp) are determined at the hardening stage when damage effect is very small and can be neglected; and (iv) the damage parameters (S, s, f) are determined using the softening stage of the stress–strain curve. This is performed by integrating the differential constitutive equations on a single material point submitted to the tensile loading path, thanks to SiDoLo identiﬁcation software (Cailletaud and Pilvin, 1994). Using material parameters determined above, the real specimen is simulated by FEA and the global force– displacement curve is compared to the experimental one. If needed, the material parameters are adjusted and new FEA simulations are performed until the experimental and numerical force–displacement curves compare well. In this work, for the sake of the industrial conﬁdentiality, no material experimental results are given. Indeed, a wide range of ‘‘virtual’’ materials will be presented by choosing appropriate values of material parameters as follows. Material A0. It represents a material class characterized by a ‘‘short’’ hardening stage, i.e., the hardening saturation is obtained for cumulative plastic strain less than 5% for a stress level of B520 MPa. Figure 3 shows the stress–plastic strain response of this material obtained with E ¼ 200 GPa, n ¼ 0.3, sp ¼ 400 MPa, b ¼ 50, Q ¼ 1,000 MPa, a ¼ 100, and C ¼ 10 GPa. Concerning the damage parameters, two different fracture behaviors representing ‘‘brittle’’ and ‘‘ductile’’ fracture will be considered deﬁning two different materials A1 and A2, given in the following. Material A1: deﬁned by s ¼ f ¼ 1 and S ¼ 0.8 MPa giving (see Figure 3) a maximum stress of smaxD499 MPa and a plastic strain at fracture (ductility) of B20%. This material will be called as a ‘‘brittle’’ material since its ductility is relatively low. Material A2: deﬁned by s ¼ f ¼ 1 and S ¼ 2.8 MPa giving (see Figure 3) a maximum stress of smaxD512 MPa and a plastic strain at fracture (ductility) of B65%. This material will be called as a ‘‘ductile’’ material since its ductility is high. Material B0. It represents a material class characterized by a ‘‘large’’ hardening stage for which the hardening saturation is reached for cumulative plastic strain of B65% at a high stress level of B1,600 MPa. Figure 4 shows the

354

Computational Damage Mechanics: Application to Metal Forming Simulation [×103] 0.5

Equivalent stress

0.4

0.3

0.2 A0 A1 A2

0.1

0.0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

Cumulative plastic deformation

Figure 3 Stress vs. plastic strain curves for material A.

[×103] 1.6

Equivalent stress

1.2

0.8

0.4

0.0 0.0

B0 B1 B2

0.1

0.2

0.3

0.4

0.5

0.6

Cumulative plastic deformation

Figure 4 Stress vs. plastic strain curves for material B.

stress–plastic strain response of this material obtained with E ¼ 200 GPa, n ¼ 0.3, sp ¼ 400 MPa, b ¼ 5, Q ¼ 1,000 MPa, a ¼ 10, and C ¼ 10 GPa. Similar to material A, the damage parameters are taken to deﬁne two different submaterials B1 and B2, given in the following. Material B1 (brittle): deﬁned by s ¼ f ¼ 1 and S ¼ 10 MPa giving (see Figure 4) a maximum stress of smaxD1,018 MPa and a plastic strain at fracture (ductility) of B18.4%;

Material B2 (ductile): deﬁned by s ¼ f ¼ 1 and S ¼ 90 MPa giving (see Figure 4) a maximum stress of smaxD1,382 MPa and a plastic strain at fracture (ductility) of B65%. These ‘‘virtual’’ materials cover a wide class of real metallic materials ranging from aluminum alloys to high-strength steels. They will be used to simulate some metal forming processes. However for some cases, other material parameters corresponding to real materials will be used without giving the corresponding experimental results.

Numerical Prediction of the Damage: Some Applications Note that the overall obtained results will be discussed from the qualitative point of view excluding the quantitative aspect including the quantitative comparison with experimental data.

3.06.5.2

Some Simple Numerical Results

The ﬁrst example concerns the well-known tensile plane-strain 1 mm 3 mm specimen shown in Figure 5 and made of material A1. This example aims to illustrate the proposed predictive methodology on a very simple example for which the plastic ﬂow localization mode is well known (shear band). Also this example is very helpful for comparison between the uncoupled and fully coupled solutions against the localization mode as well as the ﬁnal ductile fracture (i.e., macroscopic crack path). Figure 6 compares the distribution of the damage inside the specimen for both fully

u

coupled and uncoupled calculations. One can clearly note that only the fully coupled case gives a realistic localization mode, i.e., the double shear bands, whereas the uncoupled calculation (using exactly the same model) leads to an unrealistic (axisymmetric) necking mode. Moreover, in the fully damaged zone the Cauchy stress is zero in the case of fully coupled calculation (Figure 7(a)) while for the uncoupled case the stress is still maximum (!) inside the ‘‘damaged’’ zone (Figure 7(b)). From this simple example one can conclude that only the fully coupled formulation gives a realistic localization mode and hence predicts a realistic ‘‘crack’’ path. Finally, it can be concluded that an essential future of the coupling between thermomechanical behavior and damage is the fact that fully damaged points should have a decreasing stiffness as the damage increases. The above results have been obtained by the iterative implicit equilibrium solver (SI solver) available in Abaqus/Std. However, for the more complex cases, e.g., metal forming processes, the DE solver available in Abaqus/Explicit is preferred. For each case, ﬁrst the SI solver is used until the convergence rate is very low; then the DE analysis is performed in order to determine the best initial time step giving approximately the same results as the SI solver. 3.06.5.3

t Displacement control

Figure 5 FE representation of the plane-strain tensile specimen.

355

Damage Prediction in Sheet Metal Forming

Numerical simulation in sheet metal forming is increasingly used in the industrial world in order to improve the associated processes (see the periodically organized international conference NUMISHEET dedicated to the virtual sheet metal forming). In this section, the utility of the proposed coupled methodology is applied to some sheet metal forming processes in order

Figure 6 Predicted localization modes for plane-strain tensile specimen: (a) fully coupled calculation and (b) uncoupled calculation.

356

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 7 Stress distribution for plane-strain tensile specimen: (a) fully coupled calculation and (b) uncoupled calculation.

to improve them with respect to ductile damage occurrence. Two cases will be considered: (i) damage enhancement in order to simulate sheet metal cutting, and (ii) damage minimization in order to obtain an ‘‘undamaged’’ part. 3.06.5.3.1

Sheet metal cutting

In sheet metal cutting processes as blanking, slitting, shearing, or punching the ductile damage is designed into an operation by generating controlled macroscopic cracks through the intensively shear deformed zones. To show the capability of the proposed fully coupled approach to predict correctly these types of processes, a simple example of orthogonal cutting of a 3.0 mm thickness and 75.0 mm diameter (axisymmetric shape) sheet is worked out. The punch–die clearance is taken constant and equal to 0.3 mm. Moreover, both the punch and the die (tools) are supposed as rigid bodies. Because the punch–die clearance is very small (0.3 mm) compared to the blank diameter (75 mm), the 2D plane-strain hypothesis is adopted with a symmetry around the central axis. Linear quadrangular and triangular FEs (CPE3 and CPE4 for plane strain of ABAQUS library) are used with a regular mesh reﬁnement inside the area concerned with the cutting operation (i.e., located between the die and punch cutting angles). In order to study the effect of the tool wear, both punch angle radius and the die angle radius have been changed in order to schematize the ‘‘new tools’’ with 0.1 mm radius and the ‘‘worn tools’’ 1 mm cutting edge radius. The punch moves vertically with a constant displacement rate of 1 mm s–1, whilst the die and the blank-holder (rigid body) are completely ﬁxed in space. Despite the anisotropy which generally characterizes thin sheets, the present simula-

tions have been performed using the fully isotropic material model neglecting the kinematic hardening, deﬁned by: E ¼ 200 GPa, n ¼ 0.3, sp ¼ 400 MPa, Q ¼ 1,000 MPa, b ¼ 5, a ¼ C ¼ 0, s ¼ b ¼ 1, and S ¼ 5 MPa. This material exhibits nonlinear isotropic hardening with a maximum Cauchy stress for uncoupled case sN ¼ 600 MPa when the hardening is saturated for ep E85%: For the coupled case, and according to the damage-induced softening, the maximum stress is about smax ¼ 480 MPa for ep E16% while the ductility (plastic strain at fracture) is B32.5%. The friction coefﬁcient between the blank and tools is taken as Z ¼ 0:3 representing ‘‘dry’’ friction without lubrication. The variation of the blanking force vs. the punch displacement is represented in Figure 8 for both new and worn tools. As expected, the cutting operation is clearly easier with new tools: the maximum clearing force is B900 N for new tools and 1,090 N for worn tools. Similarly, the blanking operation is terminated after 1 mm of punch displacement with new tools and after 1.8 mm for worn tools. Note that the cutting operation is assumed to be terminated when the global punch-force is zero. However, the qualitative shape of these curves compares well with the experimentally observed one as can be found in Mielnick (1991). One can easily verify that the Cauchy stress inside the cracked zones is zero giving a physically sound crack according to the fully coupled model Figure 9 (see also Figure 10). These results, obtained with a simplest sheet metal cutting example, show the capability of the proposed approach to simulate numerically any other cutting process and allows its improvement virtually by studying the effect of any major parameter associated with this kind of processes (see Cherouat and Saanouni, 2003).

Numerical Prediction of the Damage: Some Applications

357

1.20E+03

Blanking force (N)

1.00E+03 8.00E+02 6.00E+02 4.00E+02

New tools

Worn tools

2.00E+02 0.00E+00

0

0.2

0.4

0.6 0.8 1 1.2 1.4 Punch displacement (mm)

1.6

1.8

2

Figure 8 Blanking force vs. punch displacement curves for both new and worn tools.

Figure 9 Macroscopic crack for punch displacement: (a) d ¼ 0.73 mm for new tools and (b) d ¼ 1.78 mm for worn tools.

3.06.5.3.2

Hydroforming of thin tubes

With the availability of advanced technology, hydroforming of thin parts has become an economic alternative to various kinds of classical stamping processes. Particularly, tube hydroforming is a well-known and widely used

technology for mass production, due to advancements in computer controls and highpressure hydraulic systems. The tendency of getting more and more geometrically complicated parts demands a systematic numerical simulation of the hydroforming processes. This allows the ‘‘virtual’’ modiﬁcation of the process

358

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 10 Stress distribution for punch displacement: (a) d ¼ 0.73 mm for new tools and (b) d ¼ 1.78 mm for worn tools. Note the zero stress values inside the fully damaged zones (macroscopic crack or cutting area).

conditions to ﬁnd the best process parameters for the ﬁnal product and gives an efﬁcient way to reduce the cost and time. However, the main difﬁculty in many hydroforming processes is to ﬁnd the best evolution of the imposed internal pressure and tools axial displacement, to avoid the buckling and the fracture of the tube during the process (tube explosion). The proposed modeling gives a powerful numerical tool for virtual improvement of hydroforming processes with respect to ductile damage occurrence before their physical realization (Hammi, 2000; Saanouni et al., 2001a, 2001b, 2001c). For illustration purposes, a geometrically complex tube hydroforming, represented in Figure 11 and taken from Thuvesen and Rask (1999), is worked out using the proposed methodology. The blank tube has an initial length of 245 mm, an inner diameter of 44.7 mm, and an outer diameter of 48.6 mm. The length of the hydroformed tube at the end of the process is 210 mm. The tube is made of aluminum alloy deﬁned by: E ¼ 84.0 GPa, n ¼ 0.3, Q ¼ 600.0 MPa, b ¼ 3.0, sy ¼ 120.0 MPa, a ¼ b ¼ 1.0, and S ¼ 10.0 MPa and 200 MPa. This gives a material response presented in Figure 12 where uncoupled and fully coupled simulations are compared. With S ¼ 10 MPa

the maximum stress is B215 MPa reached for 30% of plastic deformation and the ﬁnal fracture is obtained for 50% of plastic deformation (ductility). However, with S ¼ 200 MPa the material ductility is very large B280% and the maximum stress is B290 MPa attained for 95% of plastic deformation. Moreover, for the uncoupled case the isotropic hardening is saturated at 200% of plastic deformation giving a limit stress of B325 MPa. The controlled process parameters are the internal ﬂuid pressure and the axial feed of the two ends of the tube. The main purpose of the numerical simulation is to ‘‘virtually’’ ﬁnd the best evolution of the imposed parameters (pressure and axial displacement), giving a ﬁnal tube free from any damage, i.e., to avoid the explosion of the tube. An important parameter is the material ductility (That is, value of the accumulated plastic strain at fracture) which is governed in the used model by the material constant S controlling the ductility of the material (see Equation (100)). After many numerical simulations with the moderately ductile aluminum (S ¼ 10.0 MPa), the best evolution of the imposed internal pressure and displacement of the tube ends, which avoid any tube buckling, are given in Figure 13. While the axial displacement has a

Numerical Prediction of the Damage: Some Applications

359

Figure 11 Schematic representation of the process: (a) initial tube and (b) the shape of the tool which deﬁnes the ﬁnal tube shape.

350 300

σeq (MPa)

250 200 150 100

Coupled S = 200 Coupled S = 10 Uncoupled

50 0

0

0.5

1

1.5 εpeq

2

2.5

3

Figure 12 Material response: uncoupled vs. coupled simulations.

Internal pressure (MPa)

35 Pressure

30

30

Displacement

25

25

20

20

15

15

10

10

5

5

0

0

1

2

3

5 4 Time (s)

6

7

8

0

Axial displacement of tube ends (mm)

35

Figure 13 The best imposed loading paths for the tube hydroforming.

constant rate, the internal pressure is found to vary during the process. By using the material constants given above together with the loading parameters given in Figure 13, two fully coupled calculations have been worked out, namely: with moderate ductility material (S ¼ 10.0 MPa) and large ductility material (S ¼ 200.0 MPa). Figures 14(a)–(c) show the predicted damage distribution obtained with the S ¼ 10.0 MPa showing a macroscopic crack initiation after 12.0 mm displacement of the tube ends (that is, at 35% of the process, Figure 14(b). This crack, located at the spherical area of the tube, grows until the ﬁnal explusion of the tube (Figure 14(c).)

360

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 14 (a) Damage maps for moderate material ductility at: (a) 27% of the process; (b) 35% of the process; and (c) 90% of the process. (d) Final deformed conﬁguration (at fracture).

To avoid this crack initiation the process has been simulated with very large material ductility represented by S ¼ 200.0 MPa. Figures 15(a) and (b) compare numerical and experimental (Figure 15(b) from Thuvesen and Rask, 1999) ﬁnal conﬁgurations of the tube. From

these ﬁgures it is clear that the maximum value of damage is still lower than 0.04 and distributed as indicated in Figure 15(a). On the top of the spherical part the damage is less than 0.02 which is considered as sufﬁciently low.

Numerical Prediction of the Damage: Some Applications

Figure 14

3.06.5.4

Damage Prediction in Bulk Metal Forming

Now, the proposed fully coupled methodology will be applied to improve some bulk metal forming by avoiding the ductile damage occurrence in order to obtain a ﬁnal part free from any defect. Some examples are now discussed.

361

(continued)

3.06.5.4.1

Cold extrusion of a 3D spider

Let us now examine the workability of the spider taken from the mechanical system represented in Figure 16(a) (see Saanouni and Hammi, 2000). The so-called radial cold extrusion process schematized in Figure 16(b) is used. In fact, the posed problem is to

362

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 15 (a) Numerically predicted ﬁnal shape of the hydroformed tube for large ductility material Dmax ¼ 0.02. (b) Experimental ﬁnal shape of the hydroformed tube (source Thuvesen and Rask, 1999).

improve the workability of the spider with respect to ductile damage occurrence. We start from a cylindrical part (36 mm long and 14 mm diameter) made of the material A deﬁned in Figure 3 with the Coulomb friction parameter taken equal to: Z ¼ 0:05 (lubricated contact) between the die and the part; and

Z ¼ 0:5 (stick contact) between the part and both upper and lower punches. Figures 17(a) and (b) show the initial and the ﬁnal conﬁguration of the part, respectively. 13,920 eight-noded trilinear hexahedral elements have been used to discretize the overall workpiece without considering any symmetry (three planes of symmetry in this case).

Numerical Prediction of the Damage: Some Applications

Figure 16

363

Formulation of the spider forming problem.

Figure 17 FE modeled initial and ﬁnal conﬁguration of the part.

Predicted results obtained with the ‘‘small’’ ductility material A1 are presented in Figures 18–20 where we can see that four cracks (fully damaged zones) have been initiated after tool displacement of 8.4 mm corresponding to 63% of the overall process (Figure 20). When the material A2 with 63.0% ductility is used, the process is terminated with a relatively low value of the maximum damage (Dmax ¼ 0.22) as can be veriﬁed in Figures 21– 23. It is worth noting that the maximum accumulated plastic strain for the ductile material is B300% located inside the fully damaged zones (Saanouni and Hammi, 2000). 3.06.5.4.2

Hot forging of 3D billet

This example is worked out with an aluminum alloy having a ductility of 24% and deﬁned by: E ¼ 70.0 GPa, n ¼ 0.33 for elastiQ ¼ 2,500.0 MPa, city; sy ¼ 355.0 MPa, b ¼ 31.0 for plastic behavior with isotropic hardening; and s ¼ 1.0, b ¼ 48.0, S ¼ 250.0 MPa for damage (Saanouni et al., 2001a, 2001b, 2001c). The contact between the

part and dies is assumed to be fully lubricated (sliding contact) with Z ¼ 0:05: Both the lower and upper dies are assumed to be rigid bodies and meshed as shown in Figure 24. The lower die is ﬁxed while the upper one moves vertically with 1 mm s–1 velocity. After some numerically fully coupled simulations, the best process plan has been obtained. Figure 25 shows three different deformed conﬁgurations corresponding approximately to 5%, 20%, and 100% of the process plan. The corresponding distributions of the ductile damage are shown in Figures 26–28 from which it can be noted that the damage concentrates only on the head of the billet where a central hole will be machined. This damage distribution is found to be quite different in the case of uncoupled analysis as can be clearly seen in Figures 29–31. The ﬁrst damage initiation takes place within the billet hand (Figure 29) much later than in coupled case. Many other cracks initiate simultaneously in the billet head as well as along the billet hand (Figure 30) and propagate to cover the

364

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 18 Damage distribution after 1.0 mm punch displacement (brittle).

Figure 19 Damage distribution after 6.3 mm punch displacement (brittle).

major part of the billet (Figure 31). This conﬁrms that only the fully coupled procedure gives somewhat qualitatively realistic solution. All these results even if they are obtained with simple isotropic damage give realistic results from a qualitative point of view. They conﬁrm the predictive capability of the proposed fully coupled approach and indicate the possibility to use it as a virtual tool to improve any metal forming process with respect to the ductile damage occurrence. This approach

works well either to avoid the damage initiation or to enhance the damage initiation and growth. However many improvements of this kind of models still need to be developed in order to take into account many important phenomena such as the induced anisotropy (by the plastic ﬂow or the damage), thermal coupling (mainly for bulk metal forming), mesh dependency, etc. These aspects, under progress in our laboratory, will be brieﬂy discussed in Section 3.06.6 as perspectives for this work.

Concluding Remarks and Perspectives

Figure 20

365

Damage distribution after 8.4 mm punch displacement (brittle).

Figure 21 Damage distribution after 5.4 mm punch displacement (ductile).

3.06.6 3.06.6.1

CONCLUDING REMARKS AND PERSPECTIVES Concluding Remarks

For many decades, analytical (closed form solutions) methods based on slab methods, bounding methods, and slip-line ﬁeld methods have been used for predicting the metal ﬂow in metal forming operations. The trial-error-

based experimental method is generally used to validate these closed-form predictive approaches in order to improve existing metalforming processes or to develop new ones. Since the early 1980s the focus of metal forming research has changed signiﬁcantly in nature. The changing economic imperatives of the industrial world along with the increasing cost of experimentation focused work on process improvement and optimization based

366

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 22 Damage distribution after 10.0 mm punch displacement (ductile).

Figure 23 Damage distribution after 13.4 mm punch displacement (ductile).

on FEA. This task was signiﬁcantly facilitated by the considerable progress in the following ﬁelds related to metal forming. (i) The availability of more and more ‘‘accurate’’ mathematical models to describe the main physical phenomena (thermal, mechanical, metallurgical, etc.) involved during the large deformation of various materials. Particularly, the ‘‘multiphysical’’ coupling between these physical phenomena in the framework of the continuum mechanics leads to a

highly nonlinear set of 3D constitutive equations (ODE), convenient for implementation into specialized or general purpose FE codes. (ii) The disposability of very efﬁcient numerical methods to solve these highly nonlinear evolution problems (IBVP) using the FEA. A proper space and time discretization of the IBVP associated with efﬁcient numerical methods to solve the associated PDEs as well as to integrate the fully coupled ODEs. The increasing capabilities of the CAD/CAM tools to

Concluding Remarks and Perspectives

367

Figure 24 Representation of the billet forging tools.

Figure 25 Some deformed conﬁgurations of the billet.

represent interactively the complicated geometry of parts and tools associated with the initial meshing, remeshing, and adaptive meshing facilities. (iii) The spectacular increase in computer power including the parallel computing machines leads to a continuous decrease in the CPU time cost. Based on these favorable arguments more and more FE-based computer codes are available for engineers for applications to various metal forming processes. However, some tools toward constructing real ‘‘integrated virtual’’ metal forming systems are under progress in some research laboratories. The main goal is to offer integrated systems able to realize virtually any metal forming processes in order to

optimize ‘‘automatically’’ (with less and less human intervention) their plan before the physical realization. However the most known metal forming simulation codes are based on very simple constitutive equations accounting for linear thermoelasticity, Norton–Hoff viscoplasticity, and/or Prandlt–Reuss or total Hencky type plasticity with power-law isotropic hardening. This work is a modest contribution towards this goal by adding the possibility, for metal forming practitioners, to optimize virtually the process of concern with respect to ductile damage occurrence. Advanced fully coupled constitutive equations accounting for the coupling between thermal effects, elasticity, ﬁnite (visco)plasticity with isotropic and kinematic

368

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 26 Ductile damage localization inside the billet for fully coupled analysis (5%).

Figure 27

Ductile damage localization inside the billet for fully coupled analysis (20%).

nonlinear hardening and the ductile damage have been proposed. Related numerical aspects have been developed and fully discussed. In the case of isotropic damage, the proposed modeling has been shown to be useful to improve many sheet or bulk metal forming processes, either by avoiding the damage occurrence or enhancing and improving the damage initia-

tion and growth. All the presented results are obtained using the single surface and fully isotropic version of our model (see Section 3.06.3.4.2) under the isothermal condition. Research work is under progress in our laboratory aiming to implement the twosurface model including: (i) the initial and induced anisotropy (plasticity and damage),

Concluding Remarks and Perspectives

369

Figure 28 Ductile damage localization inside the billet for fully coupled analysis (100%).

Figure 29

Ductile damage localization inside the billet for uncoupled analysis (80%).

(ii) the ‘‘strong’’ thermal coupling, and (iii) damage-gradient regularization (see perspectives hereafter). 3.06.6.2

Perspectives

To improve the proposed fully coupled methodology further improvements still need to be made in order to offer user-friendly and interactive integrated virtual metal forming systems. These concern theoretical, numerical,

as well as geometrical aspects of metal forming simulation. The most important problem is related to the induced damage softening or negative hardening, which can be viewed as a competing phenomenon against the positive hardening due to the plastic ﬂow. As shown by applications above, when employed with FEA, these models lead inevitably to the localization of the plastic ﬂow and consequently to the damage even if the structure is initially homogeneous.

370

Computational Damage Mechanics: Application to Metal Forming Simulation

Figure 30 Ductile damage localization inside the billet for uncoupled analysis (95%).

Figure 31

Ductile damage localization inside the billet for uncoupled analysis (100%).

It has been shown by all investigators in this ﬁeld that this softening-induced ﬂow localization leads to the emergence of numerical instabilities and spurious mesh sensitivities (see the review given in Hammi, 2000). It is now well established that this undesirable behavior is related to an induced change in the IBVP during the loading process, i.e., the loss of ellipticity of the IBVP. Accordingly, considerable attention has been paid to this problem since the end of 1980s and many ‘‘localization limiters’’ have been proposed in order to regularize the IBVP under consideration. Particularly, for IBVP characterized by a

progressive decrease of stiffness with increasing load, basically two families of localization limiters or regularization methods have been proposed. (i) Nonlocal or integral regularization in which internal variables causing the softening behavior are regularized using an appropriate integral over a local interaction volume surrounding the material point (Gauss point) under concern (Bazant et al., 1984). In damage mechanics, this idea has been used independently by Pijaudier-Cabot and Bazant (1987) for concrete structure and by Saanouni and Lesne (1987) for elasto-viscoplastic (creep)

Concluding Remarks and Perspectives structure. It has been developed and used by many authors to solve various problems in damage mechanics (see papers presented in Mazars and Bazant, 1989; EuroMech378, 1998; EuroMech417, 2000; Pijaudier-Cabot and Benallal, 1993; Leblond et al., 1994; de Vree et al., 1995; Stro¨mberg and Ristinmaa, 1996; Nilsson, 1997; Jirasek, 1998; Mroz and Seweryn, 1998; Ganghoffer et al., 1999; Ganghoffer and de Borst, 2000; Comi, 2001). Note that all these nonlocal models can be viewed as a particular case of the straightforward formulation of nonlocal continuum mechanics framework pioneered by Eringen (Eringen, 1972, 1981). (ii) Higher-gradient regularization: In classical continuum mechanics the ﬁrst displacement gradient (deﬁning the strain tensor) is sufﬁcient to completely characterize the kinematics of local transformation. This simple continuum theory can be generalized by using additional state variables deﬁned by higher displacement gradients as proposed and developed since the 1960s (see Mindlin, 1965). This idea has been extensively and variously used in the literature (see papers presented in Mazars and Bazant, 1989; EuroMech378, 1998; EuroMech417, 2000; Zbib and Aifantis, 1989; Fleck et al., 1994; Benallal and Tvergaard, 1995; de Borst et al., 1996, 1999; Comi and Corigliano, 1996; Comi and Perego, 1996; Fremond and Nedjar, 1996; Fleck and Hutchinson, 1997; Li and Cescotto, 1997; Svedberg and Renesson, 1997, 1998; Faciu et al., 1998; Geers et al., 1998; Lorentz and Andrieux, 1999; Mahnken and Kuhl, 1999; Steinmann, 1999; Menzel and Steinmann, 2000; Bassani, 2001, among many others). In fact, the nonlocal as well as the highergradient formulations deﬁne a general framework useful to model the behavior of continuum with microstructure initially introduced by the brothers E. Cosserat and F. Cosserat (1909) and developed later by Toupin (1962), and Green and Rivlin (1964). Note that some works have tried to establish a link between these different frameworks (Borino et al., 1998; de Borst et al., 1998; Peerlings et al., 2001). The objective in metal forming is to generalize the fully coupled models presented in Section 3.06.3 by adding the ﬁrst-order damage gradient in the state or dissipation potentials. This work is under progress and should be applied to metal forming in the near future. The second challenge is the use of polycrystalline plasticity models in metal forming simulation. In fact, it is clear from metal forming, where plastic deformations are very large, that metallurgical aspects of materials (micrsotructural evolution) play a prominent role. Accordingly, it is legitimate to take into

371

account these microstructural transformations and their interaction with thermomechanical ﬁelds during the metal forming (Mathur and Dawson, 1989; Kalidindi et al., 1992; Kalidindi and Anand, 1994; Beaudoin et al., 1994; Munhoven et al., 1995; Nakamachi and Dong, 1996, 1997). Since the early 1990s, various reﬁnements were introduced to enrich continuum polycrystalline plasticity models with microstructural ingredients as crystallographic texture and the microstructure morphology and topology. They are mainly based either on polycrystal Taylor simulations or polycrystal plasticity self-consistent simulation (see the general monographs dedicated to this ﬁeld as Mura, 1987; Havner, 1992; Nemat-Nasser and Hori, 1993; Khan and Huang, 1995 (chapters 9–11); Raabe, 1998 and references therein). Undoubtedly the use of these physically based models to describe the microstructural changes in metal forming (rotation lattice, texture induced anisotropy, grain morphology, etc.) will be increasing in the next few years. The objective is to optimize metal forming processes in order to obtain a desired distribution of the texture (morphology and crystallographic orientation). Other open problems in metal forming simulation are summarized as follows. (i) The interaction between the 3D adaptive meshing using the damage-based error estimators with the possibility to drop the fully damaged elements. (ii) Automatic optimization of the forming conditions or the process plan with respect to some key parameters. These are either geometrical parameters such as the shape of tools, mechanical ﬁelds as damage values, or the number of pre-forms needed to achieve very complex processes.

ACKNOWLEDGMENTS Numerical results presented in this chapter are some among the many obtained by the Computational Nonlinear Mechanics team of the LASMIS laboratory at the University of Technology of Troyes (France). Thanks are due to Ph.D. students, especially Y. Hammi, who implemented the ﬁrst umat and vumat within Abaqus. Special thanks are due to our colleague Dr. A. Cherouat for performing the computation of several examples. The ﬁnancial support of the Conseil Regional Champagne Ardenne and the FEDER via the PMMC network, and the support of some companies as IRSID/Maizie`re-les-Metz (Dr. P. Autesserre) and the CETIM/Saint Etienne (Dr. M. Cristescu) are gratefully acknowledged.

372 3.06.7

Computational Damage Mechanics: Application to Metal Forming Simulation REFERENCES

A. K. Abdali, K. Benkrid and P. Bussy, 1995, Numerical simulation of sheet cutting. In: ‘‘Numiform’95,’’ eds. S. Shen and P. Dawson, MAB, pp. 807–813. M. Ainsworth and J. T. Oden, 1997, A posteriori error estimation in ﬁnite element analysis. Comp. Meth. Appl. Mech. Eng., 142, 1–88. N. Aravas, 1986, The analysis of void growth that leads to central burst during extrusion. J. Mech. Phys. Solids, 34, 55–79. I. Babuska and W. C. Rheinboldt, 1978, A posteriori error estimates for ﬁnite element method. Int. J. Num. Meth. Eng., 12, 1597–1615. I. Babuska and W. C. Rheinboldt, 1982, Computational error estimates and adaptive processes for some nonlinear structural problems. Comp. Meth. Appl. Mech. Eng., 34, 895–937. T. J. Baker, 1989, Developments and trends in 3-D mesh generation. Appl. Num. Math., 5, 275–304. J. L. Bassani, 2001, Incompatibility and a simple gradient theory of plasticity. J. Mech. Phys. Solids, 49 , 1983– 1996. K. J. Bathe, 1981, ‘‘Finite Element Procedures in Engineering Analysis,’’ Prentice-Hall, Englewood Cliffs, NJ. Z. P. Bazant, T. Belytschko and T. P. Chang, 1984, Continuum theory for strain softening. J. Eng. Mech. ASCE., 110, 1666–1692. A. Beaudoin, P. Dawson, K. Mathur, U. Kocks and D. Korzekwa, 1994, Application of polycrystal plasticity to sheet forming. Comp. Meth. Appl. Mech. Eng., 117, 49– 70. T. Belytschko and J. I. Lin, 1987, A 3-D impactpenetration algorithm with erosion. Comp. Struct., 25, 95–104. T. Belytschko, W. K. Liu and B. Moran, 2000, ‘‘Nonlinear Finite Elements for Continua and Structures,’’ Wiley, New York. A. Benallal, R. Billardon and I. Doghri, 1988, An integration algorithm and the corresponding consistent tangent operator for fully coupled elastoplastic and damage equations. Comm. Appl. Num. Methods, 4, 731– 740. A. Benallal and V. Tvergaard, 1995, Nonlocal continuum effects on bifurcation in the plane strain tensioncompression test. J. Mech. Phys. Solids, 43, 741–770. B. Bennani and J. Oudin, 1995, Backward extrusion of steels, effects of punch design on ﬂow mode and void volume fraction. Int. J. Machine Tools Manuf., 35, 903– 911. A. A. Benzerga, J. Besson and A. Pineau, 1999, Coalescence-controlled anisotropic ductile fracture. J. Eng. Mater. Technol., 121, 221–229. F. M. Beremin, 1981, Cavity formation from inclusions in ductile fracture of A508 steel. Metall. Trans., 12A, 723– 731. J. Besson, G. Cailletaud, J. L. Chaboche and S. Forest, 2001, ‘‘Me´canique non line´aire des mate´riaux,’’ Hermes, Paris, pp. 143–216. S. Bezzina and K. Saanouni, 1996, Theoretical and numerical modeling of sheet metal shear cutting. In: ‘‘MatTech’96,’’ ed. J. Lu, ITT International, Paris, pp. 15–24. H. G. Bock, 1983, Recent advances in parameter identiﬁcation techniques for ODE. In: ‘‘Numerical Treatment of Inverse Problems,’’ eds. P. Deuﬂhard and E. Hairer, Birkhauser, Boston. J. Bonet and R. D. Wood, 1997, ‘‘Nonlinear Continuum Mechanics for Finite Element Analysis,’’ Cambridge University Press, Cambridge. N. Bontcheva and R. Iankov, 1991, Numerical investigation of the damage process in metal forming. Eng. Fract. Mech., 40, 387–393.

H. Borouchaki, A. Cherouat, K. Saanouni and P. Laug, 2002, Remaillage en grandes deformation. Application a` la mise en forme de structures 2D. Rev. Euro. Ele´m. Finis, 11(1), 57–79. R. de Borst, A. Benallal and O. Heeres, 1996, A gradientenhanced damage approach to fracture. J. de Physique IV, 6, 491–502. R. de Borst, M. G. D. Geers, R. H. J. Peerlings and A. Benallal, 1998, Some remarks on gradient and nonlocal theories. In: ‘‘Damage Mechanics in Engineering Materials,’’ eds. G. Z. Voyiadijs, J. W. Ju and J. L. Chaboche, Elsevier, Oxford, pp. 223–236. R. de Borst, J. Pamin and M. G. D. Geers, 1999, On coupled gradient-dependent plasticity and damage theories with a view to localization analysis. Euro. J. Mech. A/Solids, 18, 939–962. N. Boudeau and J. C. Gelin, 1994, Prediction of the localized necking in 3D sheet metal forming processes from FE simulations. J. Mater. Proc. Technol., 45, 229– 235. D. Brokken, W. Brekelmans and F. Baaijens, 1998, Numerical modelling of the metal blanking process. J. Mater. Proc. Technol., 83, 192–199. M. Brunet, F. Sabourin and S. Mgnib-Touchal, 1996, The prediction of necking and failure in 3D sheet forming analysis using damage variable. J. de Physique IV, 6, 473–482. M. Brunet, F. Morestin and H. Walter, 2001, Numerical analysis of failure in sheet metal forming with experimental validation. Rev. Euro. Ele´m. Finis, 2-3-4, 275– 293. H. D. Bui, 1992, ‘‘Introduction aux proble`mes inverses en me´canique des mate´riaux,’’ Eyrolles, Paris. H. D. Bui and M. Tanaka (eds.), 1994, ‘‘Inverse Problems in Engineering Mechanics,’’ Springer, Berlin, New York. G. Cailletaud and P. Pilvin, 1994, Identiﬁcation and inverse problems related to material behaviour. In: ‘‘Proceedings of the 2nd International Symposium on Inverse Problems, Clamart, France, November 1994,’’ pp. 79–86. J. L. Chaboche, 1988, Continuum damage mechanics: Parts I and II. J. Appl. Mech., 55, 59–72. J. L. Chaboche and G. Cailletaud, 1996, Integration methods for complex plastic constitutive equations. Comp. Meth. Appl. Mech. Eng., 133, 125–155. R. Charlier and A. M. Habraken, 1990, Numerical modelisation of contact with friction phenomena by the ﬁnite element method. Comp. Geotech., 9, 59–72. A. B. Chaudhary and K. J. Bathe, 1986, A solution method for static and dynamic analysis of 3-D contact problems with friction. Comp. Struct., 24, 855–873. W. H. Chen and P. Tsai, 1986, Finite element analysis of elastodynamics sliding contact problems with friction. Comp. Struct., 22, 925–938. J. H. Cheng, 1988, Automatic adaptive remeshing for ﬁnite element simulation of forming processes. Int. J. Num. Meth. Eng., 26, 1–18. A. Cherouat and K. Saanouni, 2003, Numerical simulation of sheet metal blanking process using a coupled ﬁnite elastoplastic damage modelling. Int. J. Form. Processes, accepted. C. Chicone and J. Gerlach, 1987, Identiﬁability of distributed parameters. In: ‘‘Inverse and Ill-posed Problems,’’ eds. H. W. Engl and C. W. Groetsch, Academic Press. C. L. Chow and T. J. Lu, 1989, On evolution laws of anisotropic damage. Eng. Fract. Mech., 34(3), 679–701. S. E. Clift, 1992, Fracture in forming processes. In: ‘‘Numerical Modelling of Material Deformation Processes,’’ eds. P. Hartley, et al., Springer, Berlin, pp. 406– 418.

References C. Comi, 2001, A nonlocal model with tension and compression damage mechanisms. Euro. J. Mech. A/ Solids, 20, 1–22. C. Comi and A. Corigliano, 1996, On uniqueness of the dynamic ﬁnite-step problem in gradient-dependent softening plasticity. Int. J. Solids Struct., 33, 3881–3902. C. Comi and U. Perego, 1996, A generalized variable formulation for gradient dependent softening plasticity. Int. J. Num. Meth. Eng., 39, 3731–3755. T. Coupez and J. L. Chenot, 1992, Large deformations and automatic remeshing. In: ‘‘Computational Plasticity,’’ eds. D. R. J. Owen, et al., Pineridge Press, Swansea, pp. 1077–1088. M. A. Crisﬁeld, 1991, ‘‘Non-linear Finite Element Analysis of Solids and Structures,’’ Wiley, Chichester, vol. 1. M. A. Crisﬁeld, 1997, ‘‘Non-linear Finite Element Analysis of Solids and Structures,’’ Wiley, Chichester, vol. 2. E. Cosserat and F. Cosserat, 1909, ‘‘The´orie des corps de´formables,’’ Hermann, Paris. J. K. Dienes, 1979, On the analysis of rotation and stress rate in deforming bodies. Acta Mech., 32, 217–232. I. Doghri, 1993, Fully implicit integration and consistent tangent modulus in elasto-plasticity. Int. J. Num. Meth. Eng., 36, 3915–3932. I. Doghri, 1995, Numerical implementation and analysis of a class of metal plasticity models coupled with ductile damage. Int. J. Num. Meth. Eng., 38, 3403–3431. M. Dyduch, S. Cescotto and A. M. Habraken, 1995, Error estimates and indicators for adaptive analysis of bulk forming. In: ‘‘Computational Plasticity. Fundamentals and Applications,’’ eds. D. R. J. Owen and E. Onate, Pineridge Press, Swansea, pp. 1355–1367. M. Dyduch, A. M. Habraken and S. Cescotto, 1992, Automatic adaptive remeshing for numerical simulations of metal forming. Comp. Meth. Appl. Mech. Eng., 101, 283–298. A. C. Eringen, 1972, On nonlocal elasticity. Int. J. Eng. Sci., 10, 233–248. A. C. Eringen, 1981, On nonlocal plasticity. Int. J. Eng. Sci., 19, 1461–1474. H. D. Espinosa, P. D. Zavattieri and G. L. Emore, 1998, Adaptive FEM computation of geometric and material nonlinearities with application to brittle failure. Mech. Mater., 29, 275–305. EuroMech Colloquium 378, 1998, ‘‘Nonlocal Aspects in Solid Mechanics,’’ Mulhouse, France (abstracts). EuroMech Colloquium 417, 2001, ‘‘NUMEDAM’00: Numerical Modelling in Damage Mechanics, Troyes, France, 2000,’’ ed. K. Saanouni, Special issue of Rev. Euro. Ele´m. Finis, vol. 10(2-3-4). C. Faciu, A. Molinari, M. Dablij and A. Zeghloul, 1998, A new rate-type gradient-dependent viscoplastic approach for stop-and-go strain band propagation. J. Phys. IV, 8, 143–150. F. Faura, A. Garcia and M. Estrems, 1998, Finite element analysis of optimum clearance in the blanking process. J. Mater. Proc. Technol., 80–81, 121–125. N. A. Fleck and J. W. Hutchinson, 1997, Strain gradient plasticity. Adv. Appl. Mech., 33, 296–361. N. A. Fleck, G. M. Muller, M. F. Ashby and J. W. Hutchinson, 1994, Strain gradient plasticity: theory and experiment. Acta Metall. Mater., 42, 475–487. L. Fourment and J. L. Chenot, 1994, Adaptive remeshing and error control for forming processes. Rev. Euro. Ele´m. Finis, 3(2), 247–279. D. Franc-ois, A. Pineau and A. Zaoui, 1993, ‘‘Comportement me´canique des mate´riaux. Volume 2: Endommagement, Me´canique de la rupture, Me´canique du contact,’’ Hermes, Paris. M. Fremond and B. Nedjar, 1996, Damage, gradient of damage and principle of virtual power. Int. J. Solids Struct., 33(8), 1083–1103.

373

L. Gallimard, P. Ladeve`ze and J. P. Pelle, 1996, Error estimation and adaptivity in elestoplasticity. Int. J. Num. Meth. Eng., 39, 189–217. J. F. Ganghoffer and R. de Borst, 2000, A new framework in nonlocal mechanics. Int. J. Eng. Sci., 38, 453–486. J. F. Ganghoffer, L. J. Sluys and R. de Borst, 1999, A reappraisal of nonlocal mechanics. Eur. J. Mech. A Solids, 18, 17–46. M. Garajeu, J. C. Michel and P. Suquet, 2000, A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids. Comp. Meth. Appl. Mech. Eng., 183, 223–246. M. Garajeu and P. Suquet, 1997, Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles. J. Mech. Phys. Solids, 45, 873–902. M. G. D. Geers, R. de Borst, W. A. M. Brekelmans and R. H. J. Peerlings, 1998, Strain-based transient-gradient damage model for failure analysis. Comp. Meth. Appl. Mech. Eng., 160, 133–153. J. C. Gelin, 1990, Finite element analysis of ductile fracture and defects formations in cold and hot forging. Ann. CIRP, 39, 215–218. J. C. Gelin, J. Oudin and Y. Ravalard, 1985, An improved FEM for the analysis of the damage and ductile fracture in cold metal forming processes. Ann. CIRP, 34(1), 209– 213. J. C. Gelin and M. Predeleanu, 1989, Finite elastoplasticity including damage—application to metal forming processes. In: ‘‘Numerical Methods in Industrial Forming Processes,’’ eds. E. G. Thompson, et al., A. A. Balkama, pp. 425–430. P. L. George, 1991, ‘‘Automatic Mesh Generation. Application to Finite Element Methods,’’ Wiley, Chichester. P. L. George and H. Borouchaki, 1997, ‘‘Triangulation de Delauney et maillage. Application a` la me´thode des e´le´ments ﬁnis,’’ Hermes, Paris. P. Germain, 1973, ‘‘Cours de Me´canique des milieux continus,’’ Masson, Paris. P. Germain, Q. S. Nguyen and P. Suquet, 1983, Continuum thermodynamics. J. Appl. Mech., 50, 1010–1020. M. Gologanu, J. B. Leblond and J. Devaux, 1993, Approximate models for ductile metals containing non spherical voids—case of axisymmetric prolate ellipsoidal cavities. J. Mech. Phys. Solids, 41, 1723–1724. M. Gologanu, J. B. Leblond and J. Devaux, 1994, Approximate models for ductile metals containing nonspherical voids. Case of axisymmetric oblate ellipsoidal cavities. J. Eng. Mater. Technol., 116, 290–297. A. E. Green and R. S. Rivlin, 1964, Multipolar continuum mechanics. Arch. Rat. Mech. Anal., 17, 113–147. A. L. Gurson, 1977, Continuum theory of ductile rupture by void nucleation and growth: Part I. Yield criteria and ﬂow rules for porous ductile media. J. Eng. Mater. Technol., 99, 2–15. A. M. Habraken and S. Cescotto, 1990, An automatic remeshing technique for ﬁnite element simulation of forming processes. Int. J. Num. Meth. Eng., 30(8), 1503– 1525. J. O. Hallquist, G. L. Goudreau and D. J. Benson, 1985, Sliding interfaces with contact-impact in large-scale Lagrangian computations. Comp. Meth. Appl. Mech. Eng., 51, 107–137. R. Hambli, 2001, Finite element model fracture prediction during sheet-metal blanking processes. Eng. Fract. Mech., 68, 365–378. Y. Hammi, 2000, Simulation nume´rique de l’endommagement dans les proce´de´s de misen en forme. Ph.D. thesis, University of Technology of Troyes, France. P. Hartley, S. E. Clift, J. Salimi-Namin, C. E. N. Sturgess and I. Pillinger, 1989, The prediction of ductile fracture initiation in metalforming using a ﬁnite element method and various fracture criteria. Res. Mech., 28, 269–293.

374

Computational Damage Mechanics: Application to Metal Forming Simulation

S. Hartmann and P. Haupt, 1993, Stress computation and consistent tangent operator using nonlinear kinematic hardening models. Int. J. Num. Meth. Eng., 36, 3801– 3814. K. Havner, 1992, ‘‘Finite Plastic Deformation of Crystalline Solids,’’ 1st edn., Cambridge University Press, Cambridge. H. D. Hibbitt, R. V. Marcal and J. R. Rice, 1970, A ﬁnite element formulation for problems of large strain and large displacement. Int. J. Solids Struct., 6, 1069–1086. R. Hill, 1971, ‘‘Mathematical Theory of Plasticity,’’ Oxford University Press, London. M. Homsi, L. Moranc-ay and J. M. Roelandt, 1996, Remeshing processes applied to the simulation of metal cutting. Rev. Euro. Ele´m. Finis, 5, 297–321. T. J. R. Hughes, 1983, Analysis on transient algorithms with particular emphasis on stability behavior. In: ‘‘Computational Methods for Transient Analysis,’’ eds. T. Belytschko and T. J. R. Hughes, North-Holland, Amsterdam, pp. 67–155. T. J. R. Hughes, 1987, ‘‘The Finite Element Method,’’ Prentice-Hall, Englewood Cliffs, NJ. T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. C. Curnier and W. Kanoknukuchai, 1976, A ﬁnite element method for a class of contact-impact problems. Comp. Meth. Appl. Eng., 8, 249–276. T. J. R. Hughes and J. Winget, 1980, Finite rotation effects in numerical integration of rate constitutive equations arising in large deformation analysis. Int. J. Num. Meth. Eng., 15, 1862–1867. M. Jirasek, 1998, Nonlocal models for damage and fracture: comparison of approaches. Int. J. Solids Struct., 35(31–32), 4133–4145. G. Johnson and D. J. Bamman, 1984, A discussion of stress rates in ﬁnite deformation problems. Int. J. Solids Struct., 20, 725–737. JMTA, 1988, Special issue, Journal de Me´canique The´orique et Applique´e, No. 1, vol. 7. S. Kalidindi and L. Anand, 1994., Macroscopic shape change and evolution of crystallographic texture in pretextured FCC metals. J. Mech. Phys. Solids, 42, 459–490. S. Kalidindi, C. Bronkhorst and L. Anand, 1992, Crystallographic texture evolution in bulk deformation processing of FCC metals. J. Mech. Phys. Solids, 40, 537–569. P. I. Kattan and G. Z. Voyiadjis, 2002, ‘‘Damage Mechanics with Finite Elements Practical Applications with Computer Tools,’’ Springer, Berlin, New York. A. S. Khan and S. Huang, 1995, ‘‘Continuum Theory of Plasticity,’’ Wiley, New York. M. Kleiber, 1989, ‘‘Incremental Finite Element Modelling in Non-linear Solid Mechanics,’’ Ellis Horwood, Chichester. M. Kojic and K. J. Bathe, 1987, The effective-stressfunction algorithm for thermo-elasto-plasticity and creep. Int. J. Num. Meth. Eng., 24, 1509–1532. D. Krajcinovic, 1984, Continuum damage Mechanics. Appl. Mech. Rev., 37(1), 1–6. P. Ladeve`ze and J. T. Oden, 1998, ‘‘Advances in Adaptive Computational Methods in Mechanics,’’ Elsevier, Amsterdam. P. Ladeve`ze and J. P. Pelle, 2001, ‘‘La ma#ıtrise du calcul en me´canique line´aire et non lineaire,’’ Hermes-Lavoisier, Paris. J. B. Leblond, G. Perrin and J. Devaux, 1994, Bifurcation effects in ductile metals with local damage. ASME J. Appl. Mech., 61, 236–242. J. B. Leblond, G. Perrin and J. Devaux, 1995, An improved Gurson-type model for hardenable ductile metals. Euro. J. Mech., A/Solids, 14(4), 499–527. E. H. Lee, 1971, Elastic plastic deformation at ﬁnite strains. J. Appl. Mech., 36, 276–279. H. Lee, K. E. Peng and J. Wang, 1985, An anisotropic damage criterion for deformation instability and its

application to forming limit analysis of metal plates. Eng. Fract. Mech., 21(5), 1031–1054. J. Lemaitre, 1992, ‘‘A course on Damage Mechanics,’’ Springer, Berlin. J. Lemaitre and J. L. Chaboche, 1990, ‘‘Mechanics of Solid Materials,’’ Cambridge University Press, Cambridge. X. Li and S. Cescotto, 1997, A mixed element method in gradient plasticity for pressure dependent materials and modelling of strain localization. Comp. Meth. Appl. Mech. Eng., 144, 287–305. T. S. Liu and N. L. Chung, 1990, Extrusion analysis and workability prediction using ﬁnite element method. Comp. Struct., 36, 369–377. E. Lorentz and S. Andrieux, 1999, A variational formulation for nonlocal damage models. Int. J. Plasticity, 15, 119–138. R. Mahnken and E. Kuhl, 1999, Parameter identiﬁcation of gradient enhanced damage models with the ﬁnite element method. Euro. J. Mech. A/Solids, 18, 819–835. R. Mahnken and E. Stein, 1994, The identiﬁcation of parameter for visco-plastic models via FEM and gradient-methods. Modell. Simulation Mater. Sci. Eng., 2, 597–616. J. Mandel, 1973, Equations constitutives et directeurs dans les milieuc plastiques et viscoplastiques. Int. J. Solids Struct., 9, 725–740. J. F. Mariage and K. Saanouni, 2001, ‘‘Simulation nume´rique de l’endommagement en Forgeage estampage,’’ Report UTT/LASMIS/01, November 2001. K. K. Mathur and P. Dawson, 1987, Damage evolution modeling in bulk forming processes. In: ‘‘Computational Methods for Predicting Material Processing Defects,’’ ed. M. Predeleanu, Elsevier, Amsterdam, pp. 251–262. K. K. Mathur and P. Dawson, 1989, On modeling the development of crystallographic texture in bulk forming processes. Int. J. Plasticity, 5, 67–94. J. Mazars and Z. P. Bazant (eds.), 1989, ‘‘Cracking and Damage. Strain Localization and Size Effect,’’ Elsevier, NY. P. B. Mellor, 1977, Forming of anisotropic sheet metal. In: ‘‘Engineering Plasticity: Theory of Metal Forming Processes,’’ ed. H. Lippman, Springer, New York, vol. 1, pp. 67, 146. A. Menzel and P. Steinmann, 2000, On the continuum formulation of higher gradient plasticity for single and polycrystals. J. Mech. Phys. Solids, 48, 1777–1796. E. D. Mielnick, 1991, ‘‘Metalworking Science and Engineering,’’ McGraw-Hill, NY. R. D. Mindlin, 1965, Second gradient of strain and surface tension in linear elasticity. Int. J. Solids Struct., 1, 417– 437. Z. Mroz and A. Seweryn, 1998, Nonlocal failure and damage evolution rule: application to a dilatant crack model. J. Phys. IV, 8, 257–268. S. Munhoven, A. M. Habraken, J. Winters, R. Schouwenaars and P. Van Houte, 1995, Application of an anisotropic yield locus based on texture to a deep drawing simulation. In: ‘‘NUMIFORM: Simulation of Material Processing: Theory, Methods and Applications, June 1995,’’ eds. S. Shen and P. Dawson, A. A. Balkama, pp. 767–772. T. Mura, 1987, ‘‘Micromechanics of Defects in Solids,’’ 2nd edn., Martinus Nijhoff Publishers, Dordrecht. J. C. Nagteggal, 1982, On the implementation of inelastic constitutive equations with special reference to large deformation problems. Comp. Meth. Appl. Mech. Eng., 33, 469–484. E. Nakamachi and X. Dong, 1996, Elastic/crystalline viscoplastic FEA of dynamic deformation of sheet metal. Int. J. Computer-aided Eng. Software, 13, 308–326. E. Nakamachi and X. Dong, 1997, Study of texture effect on sheet failure in a limit dome height test by using

References elastic/crystalline viscoplastic FEA. J. Appl. Mech., 64, 519–524. S. Nemat-Nasser and M. Hori, 1993, ‘‘Micromechanics: Overall Properties of Heterogeneous Materials,’’ NorthHolland, Amsterdam. K. Nesnas and K. Saanouni, 2002, An integral formulation of coupled damage and viscoplastic constitutive equations. Formulation and computational issues. Int. J. Damage Mech., 11, 367–397. C. Nilsson, 1997, Nonlocal strain softening bar revisited. Int. J. Solids Struct., 34, 4399–4419. J. T. Oden, 1972, ‘‘Finite Elements of Nonlinear Continua,’’ McGraw-Hill, New York. E. Onate, M. Kleiber and C. A. de Saracibar, 1988, Plastic and viscoplastic ﬂow of void containing metal—applications to axisymmetric sheet forming problem. Int. J. Num. Meth. Eng., 25, 227–251. D. R. J. Owen and E. Hinton, 1980, ‘‘Finite Elements in Plasticity, Theory and Practice,’’ Pineridge Press, Swansea. R. H. J. Peerlings, M. G. D. Geers, R. de Borst and W. A. M. Brekelmans, 2001, A critical comparison of nonlocal and gradient-enhanced softening continua. Int. J. Solids Struct., 38, 7723–7746. P. Picart, O. Ghouati and J. C. Gelin, 1998, Optimization of metal forming process parameters with damage minimization. J. Mater. Proc. Technol., 80–81, 597–601. G. Pijaudier-Cabot and Z. P. Bazant, 1987, Nonlocal damage theory. ASCE J. Eng. Mech., 113, 1512–1533. G. Pijaudier-Cabot and A. Benallal, 1993, Strain localization and bifurcation in nonlocal continuum. Int. J. Solids Strcut., 30, 1761–1775. C. Polizzotto, G. Borino and P. Fuschi, 1998, A thermodynamically consistent formulation of nonlocal and gradient plasticity. Mech. Res. Com., 25, 75–82. P. Ponte Castaneda and M. Zaidman, 1994, Constitutive models for porous materials with evolving microstructure. J. Mech. Phys. Solids, 42, 1459–1497. D. Raabe, 1998, ‘‘Computational Materials Science,’’ Wiley-VCH, Weinheim. E. Riks, 1978, The application of Newton’s method to the problem of elastic stability. J. Appl. Mech., 39, 1060– 1066. A. Rodriguez-Ferran, I. Erbos and A. Huerta, 2001, Adaptive analysis based on error estimation for nonlocaldamage models. Rev. Euro. Ele´m. Finis, 2-3-4, 193– 207. A. Rodriguez-Ferran and A. Huerta, 2000, Error estimation and adaptivity for non-local damage models. Int. J. Solids Struct., 37, 7501–7528. G. Rousselier, 1987, Ductile fracture models and their potential in local approach of fracture. Nucl. Eng. Des., 105, 97–111. K. Saanouni, J. L. Chaboche and C. Bathias, 1986, On the creep crack growth prediction by a local approach. Eng. Fract. Mech., 25(5/6), 677–691. K. Saanouni, J. L. Chaboche and P. M. Lesne, 1989, On the creep crack-growth prediction by a non local damage formulation. Euro. J. Mech. A/Solids, 8(6), 437–459. K. Saanouni, A. Cherouat and Y. Hammi, 2001a, Numerical aspects of ﬁnite elastoplasticity with isotropic ductile damage for metal forming. Rev. Euro. Ele´m. Finis, 2-3-4, 327–351. K. Saanouni, A. Cherouat and Y. Hammi, 2001b, Optimization of hydroforming processes with respect to fracture. In: ‘‘4th Esaform Conference on Material Forming, Lie`ge, Belgium, April 23–25, 2001,’’ ed. A. M. Habraken, vol. 1, pp. 361–364. K. Saanouni, A. Cherouat, J. F. Mariage and K. Nesnas, 2001c, On the numerical simulation of 3D forging of some complex parts considering the damage occurrence. In: ‘‘4th Esaform Conference on Material Forming,

375

Lie`ge, Belgium, April 23–25, 2001,’’ ed. A. M. Habraken, vol. 2, pp. 593–597. K. Saanouni, C. Forster and F. Ben Hatira, 1994, On the anelastic ﬂow with damage. Int. J. Damage Mech., 3, 140–169. K. Saanouni and Y. Franqueville, 1999, Numerical prediction of damage during metal forming processes. Numisheet 99, Besanc-on, September 1999, pp. 13–17. K. Saanouni and Y. Hammi, 2000, Numerical simulation of damage in metal forming processes. In: ‘‘Continuous Damage and Fracture,’’ ed. A. Benallal, Elsevier, Paris, pp. 353–363. K. Saanouni and P. M. Lesne, 1987, Non-local damage model for creep crack-growth prediction. In: ‘‘Proceedings of the MECAMAT International Seminar on High Temerature Fracture Mechanisms and Mechanics, Dourdan, France, October 1989,’’ eds. P. Bensoussan and J. P. Mascarell, EGF Publication 6, MEP, London, pp. 449–466. K. Saanouni, K. Nesnas and Y. Hammi, 2000, Damage modeling in metal forming processes. Int. J. Damage Mech., 9, 196–240. M. Samuel, 1998, FEM simulations and experimental analysis of parameters of inﬂuence in the blanking process. J. Mater. Proc. Technol., 84, 97–106. F. Sidoroff, 1973, The geometrical concept of intermediate conﬁguration and elastic–plastic ﬁnite strain. Arch. Mech., 25, 299–308. J. C. Simo and T. J. R. Hughes, 1998, ‘‘Computational Inelasticity,’’ Springer, New York. J. C. Simo and J. W. Ju, 1987, Strain and stress-based continuum damage models: II. Computational aspects. Int. J. Solids Struct., 23, 841–869. J. C. Simo and T. A. Laursen, 1993, An augmented Lagrangian treatment of contact problems involving friction. Comp. Struct., 42(1), 97–112. J. C. Simo and R. L. Taylor, 1985, Consistent tangent operators for rate independent elasto-plasticity. Comp. Meth. Appl. Mech. Eng., 48, 101–118. J. C. Simo and L. Vu-Quoc, 1986, A 3-D ﬁnite strain rod model: Part ii. Computational aspects. Comp. Meth. Appl. Mech. Eng., 58, 79–116. P. Steinmann, 1999, Formulation and computation of geometrically nonlinear gradient damage. Int. J. Num. Meth. Eng., 46, 757–779. L. Stro¨mberg and M. Ristinmaa, 1996, FE-formulation of a nonlocal plasticity theory. Comp. Meth. Appl. Mech. Eng., 136, 127–144. P. Suquet, 1993, On bounds for the overall potential of power law materials containing voids with arbitrary shape. Mech. Res. Commun., 19, 51–58. T. Svedberg and K. Runesson, 1997, A thermodynamically consistent theory of gradient-regularized plasticity coupled to damage. Int. J. Plasticity, 13, 669–696. T. Svedberg and K. Runesson, 1998, An algorithm for gradient-regularized plasticity coupled to damage based on a dual mixed FE-formulation. Comp. Meth. Appl. Mech. Eng., 161, 49–65. T. Svedberg and K. Runesson, 2000, An adaptive ﬁnite element algorithm for gradient theory of plasticity with coupling to damage. Int. J. Solids Struct., 37, 7481–7499. D. Thuvesen and I. Rask, 1999, Simulation of product engineering process—tube hydroforming. In: ‘‘Numisheet’99, September 13–17, 1999, Besanc-on, France.’’ R. A. Toupin, 1962, Elastic materials with couple-stresses. Arch. Rat. Mech. Anal., 11, 385–414. C. Truesdell and W. Noll, 1992, ‘‘The Non-linear Field Theories of Mechanics,’’ 2nd edn., Springer, Berlin (1st edn., 1965). V. Tvergaard, 1990, Material failure by void growth to coalescence. Adv. Appl. Mech., 27, 83–151.

376

Computational Damage Mechanics: Application to Metal Forming Simulation

V. Tvergaard and A. Needleman, 1984, Analysis of the cup-cone fracture in a round tensile bar. Acta Metall., 32, 157–169. J. van der Lugt and J. Hue´tink, 1986, Thermal mechanically coupled ﬁnite element analysis in metalforming processes. Comp. Meth. Appl. Mech. Eng., 54, 145–160. P. Villon, H. Borouchaki and K. Saanouni, 2000, Constrained interpolation based on plasticity criterion. In: ‘‘Proceedings of ECCOMAS, Barcelona, September 11–14, 2000,’’ CD-ROM paper no.: pdf.308. G. Z. Voyiadjis and B. Deliktas, 2000, A coupled anisotropic damage model for the inelastic response of composite materials. Comp. Meth. Appl. Mech. Eng., 183, 159–199. G. Z. Voyiadjis and P. I. Kattan, 1999, ‘‘Advances in Damage Mechanics: Metals and Metal Matrix Composites,’’ Elsevier, Amsterdam. J. H. P. de Vree, W. A. M. Brekelmans and M. A. J. van Gils, 1995, Comparison of nonlocal approaches in continuum damage mechanics. Comp. Struct., 55(4), 581–588. R. H. Wagoner and J. L. Chenot, 1997, ‘‘Fundamentals of Metal Forming,’’ Wiley, USA. K. P. Walker, 1987, A uniformly valid asymptotic integration algorithm for uniﬁed viscoplastic constitutive models. In: ‘‘Advances in Inelastic Analysis,’’ eds.

S. Nakazawa et al., ASME, AMD vol. 88 and PED vol. 28, pp. 13–27. H. T. Yang, M. Heinstein and J. M. Shih, 1989, Adaptive 2D ﬁnite element simulation of metal forming processes. Int. J. Num. Meth. Eng., 28, 1409–1428. H. Zbib and E. C. Aifantis, 1989, A gradient-dependent ﬂow theory of plasticity: application to metal and soil instabilities. Appl. Mech. Rev., 42, 295–304. Y. Y. Zhu and S. Cescotto, 1995, A fully coupled elastovisco-plastic damage theory for anisotropic materials. Int. J. Solids Struct., 32(11), 1607–1641. Y. Y. Zhu, S. Cescotto and A. M. Habraken, 1992, A fully coupled elastoplastic damage modeling and fracture criteria in metal forming processes. J. Mater. Proc. Technol., 32, 197–204. O. C. Zienkiewicz, Y. C. Liu and G. C. Huang, 1988, Error estimation and adaptivity in ﬂow formulation for forming problems. Int. J. Num. Meth. Eng., 25, 23–42. O. C. Zienkiewicz, E. Onate and J. C. Heinrich, 1981, A general formulation for coupled thermal ﬂow of metals using ﬁnite elements. Int. J. Num. Meth. Eng., 17, 1497– 1514. O. C. Zienkiewicz, R. L. Taylor, 1994, ‘‘The Finite Element Method,’’ 4th edn. McGraw-Hill, New York, vol. 1 and 2. O. C. Zienkiewicz and J. Z. Zhu, 1987, A simple arror estimator and adaptive procedure for practical engineering analysis. Int. J. Num. Meth. Eng., 24, 337–357.

Copyright 2003, Elsevier Ltd. All Rights Reserved. No part of this publication may be reproduced, stored in any retrieval system or transmitted in any form or by any means: electronic, electrostatic, magnetic tape, mechanical, photocopying, recording or otherwise, without permission in writing from the publishers.

Comprehensive Structural Integrity ISBN (set): 0-08-043749-4 Volume 3; (ISBN: 0-08-044158-0); pp. 321–376

Computational Damage Mechanics: Application to Metal Forming Simulation

Computational Damage Mechanics: Application to Metal Forming Simulation

Recommend Documents